甲基化 + RNA 双组学：MethylTree 谱系分析

作者: SeekGene

时长: 61 分钟

字数: 13.4k 字

更新: 2026-07-15

阅读: 0 次

1. 模块简介

本模块基于 MethylTree 开发，用于从单细胞 DNA 甲基化数据中提取细胞间的谱系相关信号，并结合 RNA 注释信息分析细胞群体内部的克隆结构、分化关系和状态异质性。

传统单细胞 RNA 分析主要反映细胞当前的转录状态，但难以直接区分这些细胞是否来自相同祖先细胞或是否具有共同的克隆来源。MethylTree 利用 DNA 甲基化中逐渐积累的 epimutation 信号作为一种天然的谱系记录，用于推断细胞之间的 lineage relationship。与依赖人工遗传标记或体细胞突变的方法不同，这类方法更适合在人类样本或缺乏预先标记体系的数据中开展非侵入式谱系分析。

本教程按照实际运行顺序组织：

构建 MethylTree 输入文件：将单细胞 *_allc.gz 文件和细胞注释表整理为 data.tsv.gz 和 sample_sheet.tsv.gz；
运行 MethylTree 主程序：完成相似性计算、谱系树构建、热图展示和 clone inference；
查看和解释输出结果。

在生物学应用上，MethylTree 可以帮助回答以下问题：

细胞之间是否存在共同谱系来源：通过计算细胞间甲基化相似性，识别具有相近谱系历史的细胞群体。
细胞群体内部是否存在克隆或亚克隆结构：根据甲基化相似性和 clone inference 结果，发现潜在的 methyl clone 或谱系相关细胞模块。
分化关系或细胞命运选择是否具有谱系基础：在发育、分化或肿瘤演化场景中，分析不同细胞状态之间是否存在共同祖先、分支关系或克隆扩增现象。

使用前提示
MethylTree 的分析效果高度依赖单细胞甲基化数据质量。根据软件作者建议，实践中单细胞 genome coverage 最好不低于约 2%；覆盖度过低时，构建出的 cell × genomic region 矩阵会包含大量缺失值，细胞间相似性计算不稳定，进而影响相似性热图、谱系树和 clone inference 的可靠性。
本教程默认将 CpG 甲基化信号按 500 bp genomic bin 聚合，也可根据数据情况调整为 1 kb。500 bp 和 1 kb 是 MethylTree 分析中较常用、较合适的窗口大小，可以在降低数据稀疏性的同时尽量保留局部 epimutation 信号。过大的窗口虽然可以进一步减少缺失值，但可能会平均掉局部甲基化差异，不利于谱系信号识别。

2. 构建 MethylTree 输入文件

本节将单细胞甲基化文件和细胞注释表整理为 MethylTree 主程序可直接读取的输入文件。完成后会生成：

data.tsv.gz：甲基化矩阵输入文件；
sample_sheet.tsv.gz：细胞注释与分组信息文件。

后续 MethylTree 谱系树构建和可视化均基于这两个文件进行。

2.1 准备输入数据

开始分析前，请准备以下文件：

单细胞甲基化文件：*_allc.gz，每个细胞对应一个 allc 文件，通常由标准分析结果中的 allcools.tar.gz 解压获得。
细胞注释表：meta.tsv，至少包含一列可与 allc 文件名对应的 barcode。如果需要展示细胞类型、clone 或 cluster 信息，可同时包含 CellAnnotation、clone_id、resolution.0.5_d30 等字段。

每个 allc 文件名去掉 _allc.gz 后，应能与 meta.tsv 中的 barcode 对应。例如，AAAGAAGAAGGATTGTT_allc.gz 对应的细胞 ID 为 AAAGAAGAAGGATTGTT。

如果 RNA barcode 带有样本后缀，例如 AAAGAAGAAGGATTGTT_7，而 allc 文件名中不包含 _7，可在后续参数区设置 remove_barcode_suffix = True。脚本会去掉最后一个下划线后的后缀，用于匹配 allc 文件。

2.2 allc 文件格式

本教程默认输入的 allc 文件为标准格式，每一行对应一个甲基化位点：

text

chrom    pos    strand    context    mc    cov    methylated

各字段含义如下：

字段	说明
`chrom`	染色体名称
`pos`	位点坐标
`strand`	链方向
`context`	甲基化上下文，例如 `CG`、`CHG`、`CHH`
`mc`	该位点支持甲基化的 UMI/read 数
`cov`	该位点的总 UMI/read 覆盖数
`methylated`	是否显著甲基化的标记；如果未进行显著性检验，通常为 `1`

后续计算主要使用 chrom、pos、context、mc 和 cov 字段；strand 和 methylated 保留在原始文件中，不参与 bin-level 甲基化比例计算。

本教程仅使用 CpG 甲基化信号，即保留 context 以 CG 开头的位点。

脚本会按照设定的 bin_size 将 CpG 位点聚合到固定长度的 genomic bin 中，并计算每个细胞在每个 bin 内的甲基化比例：

text

value = 100 × sum(mc) / sum(cov)

最终 data.tsv.gz 中的 value 表示某个细胞在某个 genomic bin 内的 CpG methylation rate，通常范围为 0–100。

2.3 导入依赖包

python

# ============================================================
# 2.3 导入依赖包
# ============================================================

import os
import re
import gzip
import math
import json
import shutil
import tarfile
import time
from pathlib import Path
from collections import Counter, defaultdict
from concurrent.futures import ProcessPoolExecutor, as_completed

import numpy as np
import pandas as pd

print("Python executable:", os.sys.executable)

output

Python executable: python3

2.4 设置输入构建参数

本步骤用于配置 MethylTree 输入构建阶段的核心参数，包括输入路径、barcode 匹配方式、注释列选择以及细胞和 region 过滤阈值。运行前请根据实际数据情况修改相应参数。

python

# ============================================================
# 2.4 设置输入构建参数
# ============================================================

# 样本名称，会用于输出目录命名
sample_name = "WTJW969"

# 输入文件
# meta_file：细胞注释表
# allc_dir：每个细胞 *_allc.gz 文件所在目录，脚本会递归查找
meta_file = "../input/meta.tsv"
allc_dir = "../input/allcools"

# 输出目录
# 本节会在该目录下生成 data.tsv.gz、sample_sheet.tsv.gz 和 QC 文件
outdir = "../output/01_methyltree_input"


# ============================================================
# 文档展示辅助函数
# ============================================================

def compact_path(path, max_len=100):
    """将过长路径压缩为适合文档展示的单行文本。"""
    s = str(path)
    if len(s) <= max_len:
        return s
    keep = max_len - 3
    head = max(30, keep // 2)
    tail = keep - head
    return s[:head] + "..." + s[-tail:]

def print_path(label, path, max_len=100):
    print(f"{label}: {compact_path(path, max_len=max_len)}")

def print_compact_log(log_text, max_len=100):
    """打印运行日志，同时压缩过长的路径类输出。"""
    for raw_line in str(log_text).replace("\r", "\n").splitlines():
        line = raw_line.rstrip()
        if not line:
            continue
        stripped = line.strip()
        if len(stripped) > max_len and ("/" in stripped or "\\" in stripped):
            print(compact_path(stripped, max_len=max_len))
        else:
            print(line)


# ============================================================
# metadata 列名配置
# ============================================================

# meta.tsv 中用于匹配 allc 文件名的 barcode 列
barcode_col = "barcode"

# 如果 meta.tsv 中 barcode 带有 _1、_2 等后缀，而 allc 文件名不带后缀，则设为 True
remove_barcode_suffix = False


# ============================================================
# 细胞筛选配置
# ============================================================

# None 表示分析全部细胞
# 字典表示按指定列筛选细胞，例如只分析 T_cells 中的 cluster 0
subset_filters = {
    "CellAnnotation": ["T_cells"],
    "resolution.0.5_d30": [0]
}


# ============================================================
# MethylTree 分组标签配置
# ============================================================

# clone_key_col：用于生成 sample_sheet.tsv.gz 中的 Clone_key
# 如果使用 RNA cluster，通常配合 clone_key_prefix = "c" 生成 c0、c1、c2 等标签
# 如果使用真实 clone 或细胞类型，clone_key_prefix 可设为 None
clone_key_col = "resolution.0.5_d30"
clone_key_prefix = "c"

# celltype_key_col：用于生成 Celltype_key
celltype_key_col = "CellAnnotation"

# cluster_key_col：用于生成 Cluster_key；如不需要可设为 None
cluster_key_col = "resolution.0.5_d30"


# ============================================================
# bin 聚合与过滤参数
# ============================================================

# bin_size：CpG 甲基化信号聚合窗口大小，常用 500 或 1000
bin_size = 500

# min_cell_frac_per_region：region 覆盖细胞比例过滤阈值
# 0.05 表示一个 region 至少需要在 5% 输入细胞中有覆盖，才会被保留
min_cell_frac_per_region = 0.05

# min_regions_per_cell：细胞过滤阈值
# 一个细胞在过滤后的 region 中至少覆盖该数量的 region，才会进入最终输入
min_regions_per_cell = 500

# 并行线程数
threads = 32

2.5 检查输入文件并匹配 barcode

本步骤用于检查输入文件是否完整，并将 meta.tsv 中的细胞 barcode 与对应的 *_allc.gz 文件进行匹配。脚本会扫描 allc 文件目录、读取细胞注释表、生成初始 sample_sheet，并过滤未匹配到 allc 文件的细胞。

如果缺失 allc 文件的细胞较多，建议重点检查 barcode_col、remove_barcode_suffix 和 allc_dir 参数是否设置正确。

python

# ============================================================
# 2.5.1 路径标准化与输出目录创建
# ============================================================

def clean_label(x):
    """将任意标签转换成适合用于目录名的字符串。"""
    x = str(x)
    x = x.replace("/", "-").replace("\\", "-")
    x = re.sub(r"\s+", "_", x)
    return x


def frac_to_cov_label(frac):
    """例如 0.01 -> cov1pct, 0.005 -> cov0p5pct。"""
    return f"cov{frac * 100:g}pct".replace(".", "p")


# 路径标准化
meta_file = Path(meta_file).expanduser()
allc_dir = Path(allc_dir).expanduser()
outdir = Path(outdir).expanduser()


# 根据细胞筛选条件生成输出标签
if subset_filters is None:
    subset_label = "all_cells"
else:
    label_parts = []
    for key, values in subset_filters.items():
        value_label = "_".join([clean_label(x) for x in values])
        label_parts.append(f"{clean_label(key)}_{value_label}")
    subset_label = "__".join(label_parts)


# 根据 min_cell_frac_per_region 自动生成 cov 标签
# 例如 0.01 -> cov1pct, 0.05 -> cov5pct, 0.005 -> cov0p5pct
cov_label = frac_to_cov_label(min_cell_frac_per_region)

out_label = f"{subset_label}_bin{bin_size}_{cov_label}_cell{min_regions_per_cell}"


# result_dir：本次 MethylTree 输入文件的最终输出目录
result_dir = outdir / f"{sample_name}_{out_label}"

# tmp_dir：逐细胞 bin-level 中间文件目录
# 如果重新运行且参数不变，已有单细胞中间文件会被直接复用
# 如需完全重跑，可手动删除该目录
tmp_dir = result_dir / f"intermediate_{bin_size}bp_by_cell"

result_dir.mkdir(parents=True, exist_ok=True)
tmp_dir.mkdir(parents=True, exist_ok=True)


# 参数检查
print("sample_name:", sample_name)
print("meta_file:", meta_file)
print("allc_dir:", allc_dir)
print("outdir:", outdir)
print("result_dir:", result_dir)
print_path("tmp_dir", tmp_dir)

print("\n[metadata]")
print("barcode_col:", barcode_col)
print("remove_barcode_suffix:", remove_barcode_suffix)

print("\n[cell selection]")
print("subset_filters:", subset_filters)
print("subset_label:", subset_label)

print("\n[grouping keys]")
print("clone_key_col:", clone_key_col)
print("clone_key_prefix:", clone_key_prefix)
print("celltype_key_col:", celltype_key_col)
print("cluster_key_col:", cluster_key_col)

print("\n[binning and filtering]")
print("bin_size:", bin_size)
print("min_cell_frac_per_region:", min_cell_frac_per_region)
print("cov_label:", cov_label)
print("min_regions_per_cell:", min_regions_per_cell)
print("threads:", threads)

output

sample_name: WTJW969
meta_file: ../input/meta.tsv
allc_dir: ../input/allcools
outdir: ../output/01_methyltree_input
result_dir: ../output/01_methyltree_input/WTJW969_Cel...n_T_cells__resolution.0.5_d30_0_bin500_cov5pct_cell500
tmp_dir: ../output/01_methyltree_input/WTJW969_CellAn...30_0_bin500_cov5pct_cell500/intermediate_500bp_by_cell

[metadata]
barcode_col: barcode
remove_barcode_suffix: False

[cell selection]
subset_filters: {'CellAnnotation': ['T_cells'], 'resolution.0.5_d30': [0]}
subset_label: CellAnnotation_T_cells__resolution.0.5_d30_0

[grouping keys]
clone_key_col: resolution.0.5_d30
clone_key_prefix: c
celltype_key_col: CellAnnotation
cluster_key_col: resolution.0.5_d30

[binning and filtering]
bin_size: 500
min_cell_frac_per_region: 0.05
cov_label: cov5pct
min_regions_per_cell: 500
threads: 32

python

# ============================================================
# 2.5.2 检查 allc 文件
# ============================================================

allc_files = sorted(allc_dir.rglob("*_allc.gz")) if allc_dir.exists() else []

print("allc_dir:", allc_dir)
print("allc files:", len(allc_files))

if len(allc_files) == 0:
    raise FileNotFoundError(
        f"No *_allc.gz found in allc_dir: {allc_dir}\n"
        "Please check whether allc_dir points to the directory containing per-cell *_allc.gz files."
    )

print("\nPreview allc files:")
for p in allc_files[:5]:
    print(f"- {p.name}: {compact_path(p)}")

# 从 allc 文件名中提取 cell ID，后续会用于和 meta.tsv 中的 barcode 匹配
allc_cell_ids = [Path(p).name.replace("_allc.gz", "") for p in allc_files]

print("\nunique allc cell IDs:", len(set(allc_cell_ids)))
print("duplicated allc cell IDs:", len(allc_cell_ids) - len(set(allc_cell_ids)))

print("\nPreview allc cell IDs:")
for x in allc_cell_ids[:5]:
    print(x)

output

allc_dir: ../input/allcools
allc files: 2196

Preview allc files:
../input/allcools/WTJW969_forward_AAAG_1_merged_fr_bam_allcools/AAAGAAGAAGTGAATTG_allc.gz
../input/allcools/WTJW969_forward_AAAG_1_merged_fr_bam_allcools/AAAGAAGAGAAGGTTGG_allc.gz
../input/allcools/WTJW969_forward_AAAG_1_merged_fr_bam_allcools/AAAGAAGAGAAGTAGTA_allc.gz
../input/allcools/WTJW969_forward_AAAG_1_merged_fr_bam_allcools/AAAGAAGATTATTGTGT_allc.gz
../input/allcools/WTJW969_forward_AAAG_1_merged_fr_bam_allcools/AAAGAAGATTTAGGTGG_allc.gz

unique allc cell IDs: 2196
duplicated allc cell IDs: 0

Preview allc cell IDs:
AAAGAAGAAGTGAATTG
AAAGAAGAGAAGGTTGG
AAAGAAGAGAAGTAGTA
AAAGAAGATTATTGTGT
AAAGAAGATTTAGGTGG

python

# ============================================================
# 2.5.3 读取 metadata 并生成初始 sample_sheet
# ============================================================

meta = pd.read_csv(meta_file, sep="\t")

print("metadata shape:", meta.shape)
print("metadata columns:")
print(meta.columns.tolist())

# 检查必要列
required_cols = [barcode_col, clone_key_col, celltype_key_col]
if cluster_key_col is not None:
    required_cols.append(cluster_key_col)

# subset_filters 中用到的列也需要检查
if subset_filters is not None:
    required_cols.extend(list(subset_filters.keys()))

required_cols = list(dict.fromkeys(required_cols))

missing_cols = [col for col in required_cols if col not in meta.columns]
if len(missing_cols) > 0:
    raise ValueError(f"metadata missing required columns: {missing_cols}")

meta[barcode_col] = meta[barcode_col].astype(str)

print("\nTotal cells in metadata:", meta.shape[0])

if celltype_key_col in meta.columns:
    print("\nCell type counts before filtering:")
    print(meta[celltype_key_col].value_counts(dropna=False))

if clone_key_col in meta.columns:
    print("\nclone_key_col counts before filtering:")
    print(meta[clone_key_col].value_counts(dropna=False).sort_index())


# ============================================================
# 根据 subset_filters 筛选细胞
# ============================================================

sub_meta = meta.copy()

if subset_filters is not None:
    print("\nApplying subset_filters:")
    print(subset_filters)

    for key, values in subset_filters.items():
        values = [str(x) for x in values]
        before_n = sub_meta.shape[0]

        sub_meta = sub_meta[
            sub_meta[key].astype(str).isin(values)
        ].copy()

        after_n = sub_meta.shape[0]
        print(f"{key} in {values}: {before_n} -> {after_n}")

if sub_meta.shape[0] == 0:
    raise ValueError("No cells selected. Check subset_filters.")

print("\nSelected cells after filtering:", sub_meta.shape[0])


# ============================================================
# 生成用于匹配 allc 文件的 barcode
# ============================================================

# gex_cb：保留原始 barcode
sub_meta["gex_cb"] = sub_meta[barcode_col].astype(str)

# m_cb：用于匹配 allc 文件名
if remove_barcode_suffix:
    sub_meta["m_cb"] = (
        sub_meta[barcode_col]
        .astype(str)
        .str.replace(r"_[0-9]+$", "", regex=True)
    )
else:
    sub_meta["m_cb"] = sub_meta[barcode_col].astype(str)


# ============================================================
# 生成 MethylTree 注释列
# ============================================================

# Clone_key
if clone_key_prefix is None:
    sub_meta["Clone_key"] = sub_meta[clone_key_col].astype(str)
else:
    raw_clone = sub_meta[clone_key_col].astype(str)

    # 如果原始值已经带有前缀，例如 c0/c1，则不重复添加
    if raw_clone.str.startswith(clone_key_prefix).all():
        sub_meta["Clone_key"] = raw_clone
    else:
        sub_meta["Clone_key"] = clone_key_prefix + raw_clone

# Celltype_key
sub_meta["Celltype_key"] = sub_meta[celltype_key_col].astype(str)

# Cluster_key
if cluster_key_col is None:
    sub_meta["Cluster_key"] = sub_meta["Clone_key"].astype(str)
else:
    if clone_key_prefix is not None and cluster_key_col == clone_key_col:
        sub_meta["Cluster_key"] = sub_meta["Clone_key"].astype(str)
    else:
        sub_meta["Cluster_key"] = sub_meta[cluster_key_col].astype(str)


# ============================================================
# 生成初始 sample_sheet
# ============================================================

sample_sheet = pd.DataFrame({
    "sample": sub_meta["m_cb"].astype(str),
    "HQ": True,
    "Clone_key": sub_meta["Clone_key"].astype(str),
    "Celltype_key": sub_meta["Celltype_key"].astype(str),
    "CellAnnotation": sub_meta[celltype_key_col].astype(str),
    "Cluster_key": sub_meta["Cluster_key"].astype(str),
    "orig.ident": sub_meta["orig.ident"].astype(str) if "orig.ident" in sub_meta.columns else sample_name,
    "Sample": sub_meta["Sample"].astype(str) if "Sample" in sub_meta.columns else sample_name,
    "raw_Sample": sub_meta["raw_Sample"].astype(str) if "raw_Sample" in sub_meta.columns else sample_name,
    "gex_cb": sub_meta["gex_cb"].astype(str),
    "m_cb": sub_meta["m_cb"].astype(str),
})

sample_raw_out = result_dir / "sample_sheet.raw.tsv.gz"
sample_sheet.to_csv(sample_raw_out, sep="\t", index=False, compression="gzip")


# ============================================================
# 输出检查
# ============================================================

print("\nSaved raw sample_sheet:")
print(sample_raw_out)

print("\nSelected cells:", sample_sheet.shape[0])

print("\nCelltype_key counts:")
print(sample_sheet["Celltype_key"].value_counts(dropna=False))

print("\nClone_key counts:")
print(sample_sheet["Clone_key"].value_counts(dropna=False).sort_index())

print("\nCluster_key counts:")
print(sample_sheet["Cluster_key"].value_counts(dropna=False).sort_index())

print("\nPreview sample_sheet:")
print(sample_sheet.head().to_string(index=False))

output

metadata shape: (2196, 9)
metadata columns:
['barcode', 'orig.ident', 'nCount_RNA', 'nFeature_RNA... 'raw_Sample', 'resolution.0.5_d30', 'CellAnnotation']

Total cells in metadata: 2196

Cell type counts before filtering:
CellAnnotation
T_cells 1416
Monocyte 420
B_cell 201
NK_cell 159
Name: count, dtype: int64

clone_key_col counts before filtering:
resolution.0.5_d30
0 613
1 453
2 317
3 310
4 201
5 159
6 63
7 40
8 40
Name: count, dtype: int64

Applying subset_filters:
{'CellAnnotation': ['T_cells'], 'resolution.0.5_d30': [0]}
CellAnnotation in ['T_cells']: 2196 -> 1416
resolution.0.5_d30 in ['0']: 1416 -> 613

Selected cells after filtering: 613

Saved raw sample_sheet:
../output/01_methyltree_input/WTJW969_CellAnnotation_...5_d30_0_bin500_cov5pct_cell500/sample_sheet.raw.tsv.gz

Selected cells: 613

Celltype_key counts:
Celltype_key
T_cells 613
Name: count, dtype: int64

Clone_key counts:
Clone_key
c0 613
Name: count, dtype: int64

Cluster_key counts:
Cluster_key
c0 613
Name: count, dtype: int64

Preview sample_sheet:

sample HQ Clone_key Celltype_key CellAnnotation Cluster_key \\
0 GGAAGTGTTTGGAGTTG True c0 T_cells T_cells c0
1 GTATAGTAGATGGGAGT True c0 T_cells T_cells c0
2 TGAGTGGATTATGATAG True c0 T_cells T_cells c0
4 AGGAAATAGAGGAAATA True c0 T_cells T_cells c0
6 AGTAAGATGAATTAGTG True c0 T_cells T_cells c0

orig.ident Sample raw_Sample gex_cb \\
0 WTJW969 sample_WTJW969_E sample_WTJW969_E GGAAGTGTTTGGAGTTG
1 WTJW969 sample_WTJW969_E sample_WTJW969_E GTATAGTAGATGGGAGT
2 WTJW969 sample_WTJW969_E sample_WTJW969_E TGAGTGGATTATGATAG
4 WTJW969 sample_WTJW969_E sample_WTJW969_E AGGAAATAGAGGAAATA
6 WTJW969 sample_WTJW969_E sample_WTJW969_E AGTAAGATGAATTAGTG

m_cb
0 GGAAGTGTTTGGAGTTG
1 GTATAGTAGATGGGAGT
2 TGAGTGGATTATGATAG
4 AGGAAATAGAGGAAATA
6 AGTAAGATGAATTAGTG

python

# ============================================================
# 2.5.4 扫描 allc 文件并匹配 sample_sheet
# ============================================================

def get_cell_id_from_allc_path(p):
    """从 *_allc.gz 文件名中提取 cell ID。"""
    name = Path(p).name
    return name.replace("_allc.gz", "")


# 扫描 allc 文件
allc_files = sorted(allc_dir.rglob("*_allc.gz"))

allc_df = pd.DataFrame({
    "sample": [get_cell_id_from_allc_path(p) for p in allc_files],
    "allc_path": [str(p) for p in allc_files],
})

# 检查重复 allc cell ID
dup_n = allc_df["sample"].duplicated().sum()
if dup_n > 0:
    print(f"[WARNING] duplicated allc cell IDs: {dup_n}")
    print(allc_df[allc_df["sample"].duplicated(keep=False)].head())

# 如果同一个 cell ID 有多个 allc 文件，默认保留第一个
allc_df = allc_df.drop_duplicates(subset=["sample"], keep="first").copy()

# 与 sample_sheet 匹配
matched_sample_sheet = sample_sheet.merge(allc_df, on="sample", how="left")

missing_mask = matched_sample_sheet["allc_path"].isna()
missing_sample_sheet = matched_sample_sheet[missing_mask].copy()

matched_sample_sheet = matched_sample_sheet[~missing_mask].copy()
matched_sample_sheet["sample"] = matched_sample_sheet["sample"].astype(str)

if matched_sample_sheet.shape[0] == 0:
    raise ValueError(
        "No cells matched to allc files. "
        "Check barcode_col, remove_barcode_suffix, and *_allc.gz filenames."
    )

# 保存匹配后的 sample_sheet
sample_out = result_dir / "sample_sheet.tsv.gz"
matched_sample_sheet.to_csv(sample_out, sep="\t", index=False, compression="gzip")


# ============================================================
# 输出匹配结果
# ============================================================

print("\nallc files:", len(allc_files))
print("unique allc samples:", allc_df["sample"].nunique())
print("sample_sheet cells before allc matching:", sample_sheet.shape[0])
print("matched cells:", matched_sample_sheet.shape[0])
print("missing allc:", missing_sample_sheet.shape[0])

if missing_sample_sheet.shape[0] > 0:
    print("\nPreview missing cells:")
    print(missing_sample_sheet[["sample", "gex_cb", "m_cb", "Clone_key", "Celltype_key"]].head(10).to_string(index=False))

print("\nMatched Celltype_key counts:")
print(matched_sample_sheet["Celltype_key"].value_counts(dropna=False))

print("\nMatched Clone_key counts:")
print(matched_sample_sheet["Clone_key"].value_counts(dropna=False).sort_index())

if "Cluster_key" in matched_sample_sheet.columns:
    print("\nMatched Cluster_key counts:")
    print(matched_sample_sheet["Cluster_key"].value_counts(dropna=False).sort_index())

print("\nSaved matched sample_sheet:")
print(compact_path(sample_out))

output

allc files: 2196
unique allc samples: 2196
sample_sheet cells before allc matching: 613
matched cells: 613
missing allc: 0

Matched Celltype_key counts:
Celltype_key
T_cells 613
Name: count, dtype: int64

Matched Clone_key counts:
Clone_key
c0 613
Name: count, dtype: int64

Matched Cluster_key counts:
Cluster_key
c0 613
Name: count, dtype: int64

Saved matched sample_sheet:
../output/01_methyltree_input/WTJW969_CellAnnotation_...n.0.5_d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz

2.6 生成 MethylTree 输入文件

本步骤用于将单细胞 *_allc.gz 文件转换为 MethylTree 输入格式，包括 CpG 位点筛选、bin-level 甲基化比例计算，以及低覆盖 region 和 cell 过滤。

已生成的中间文件会被自动复用，以减少重复计算。

python

# ============================================================
# 2.6.1 定义 allc 转 bin-level methylation 的函数
# ============================================================

def allc_to_bins_one_cell(allc_path, cell_id, bin_size):
    """
    Read one per-cell *_allc.gz file and aggregate CpG methylation
    UMI counts into fixed-size genomic bins.
    """
    bin_mc = defaultdict(int)
    bin_cov = defaultdict(int)

    with gzip.open(allc_path, "rt") as f:
        for line in f:
            parts = line.rstrip("\n").split("\t")
            if len(parts) < 6:
                continue

            chrom = parts[0]

            try:
                pos = int(parts[1])
            except Exception:
                continue

            context = parts[3]

            # Only keep CpG context
            if not context.startswith("CG"):
                continue

            try:
                mc = int(float(parts[4]))
                cov = int(float(parts[5]))
            except Exception:
                continue

            if cov <= 0:
                continue

            # allc position is 1-based; convert to 0-based genomic bin
            start = ((pos - 1) // bin_size) * bin_size
            end = start + bin_size
            region_id = f"{chrom}:{start}-{end}"

            bin_mc[region_id] += mc
            bin_cov[region_id] += cov

    rows = []
    for region_id, cov in bin_cov.items():
        mc = bin_mc[region_id]
        value = 100.0 * mc / cov
        rows.append((cell_id, region_id, value))

    return pd.DataFrame(
        rows,
        columns=["cell_id", "genomic_region_id", "value"]
    )


def process_one_cell_to_bin(args):
    """
    Process one cell and write/read its intermediate bin-level methylation file.
    Existing intermediate files are reused to avoid repeated parsing of allc files.
    """
    cell_id, allc_path, bin_size, out_cell_file = args

    if os.path.exists(out_cell_file) and os.path.getsize(out_cell_file) > 0:
        df_cell = pd.read_csv(out_cell_file, sep="\t", compression="gzip")
    else:
        df_cell = allc_to_bins_one_cell(allc_path, cell_id, bin_size)
        df_cell.to_csv(out_cell_file, sep="\t", index=False, compression="gzip")

    regions = df_cell["genomic_region_id"].astype(str).unique().tolist()

    return {
        "cell_id": cell_id,
        "out_cell_file": out_cell_file,
        "raw_region_n": len(regions),
        "regions": regions,
    }

python

# ============================================================
# 2.6.2 逐细胞生成 bin-level 中间文件
# ============================================================

start_total = time.time()

tasks = []
for row in matched_sample_sheet.itertuples(index=False):
    cell_id = str(getattr(row, "sample"))
    allc_path = getattr(row, "allc_path")
    out_cell_file = tmp_dir / f"{cell_id}.{bin_size}bp.tsv.gz"
    tasks.append((cell_id, allc_path, bin_size, str(out_cell_file)))

if len(tasks) == 0:
    raise ValueError("No cells to process. Check matched_sample_sheet.")

n_cells_input = len(tasks)
n_workers = min(threads, n_cells_input)

print("input cells:", n_cells_input)
print("threads:", threads)
print("workers used:", n_workers)
print_path("tmp_dir", tmp_dir)

region_cell_count = Counter()
cell_region_count_raw = {}
per_cell_files = []

start = time.time()
done = 0

with ProcessPoolExecutor(max_workers=n_workers) as ex:
    future_to_cell = {
        ex.submit(process_one_cell_to_bin, t): t[0]
        for t in tasks
    }

    for fut in as_completed(future_to_cell):
        cell_id_for_error = future_to_cell[fut]

        try:
            res = fut.result()
        except Exception as e:
            raise RuntimeError(
                f"Failed to process cell: {cell_id_for_error}"
            ) from e

        done += 1

        cell_id = res["cell_id"]
        out_cell_file = res["out_cell_file"]
        regions = res["regions"]

        per_cell_files.append(out_cell_file)
        cell_region_count_raw[cell_id] = res["raw_region_n"]
        region_cell_count.update(regions)

        if done % 50 == 0 or done == n_cells_input:
            elapsed = time.time() - start
            rate = done / elapsed if elapsed > 0 else 0
            remain = (n_cells_input - done) / rate if rate > 0 else np.nan

            print(
                f"[{bin_size}bp] binned {done}/{n_cells_input} cells | "
                f"elapsed {elapsed/60:.1f} min | "
                f"speed {rate:.2f} cells/s | "
                f"ETA {remain/60:.1f} min",
                flush=True
            )

per_cell_files = sorted(per_cell_files)

print("\nFinished per-cell binning.")
print("per-cell intermediate files:", len(per_cell_files))
print("raw regions before filtering:", len(region_cell_count))

print("\nPreview raw region counts per cell:")
print(
    pd.Series(cell_region_count_raw, name="raw_region_n")
    .sort_values(ascending=False)
    .head()
    .to_string()
)

output

python

# ============================================================
# 2.6.3 筛选 genomic regions
# ============================================================

min_cells_per_region = max(
    1,
    int(math.ceil(n_cells_input * min_cell_frac_per_region))
)

kept_regions = {
    region_id
    for region_id, observed_cell_n in region_cell_count.items()
    if observed_cell_n >= min_cells_per_region
}

raw_region_n = len(region_cell_count)
kept_region_n = len(kept_regions)
kept_region_frac = kept_region_n / raw_region_n if raw_region_n > 0 else 0

print("\n[Region filtering]")
print("input cells:", n_cells_input)
print("min_cell_frac_per_region:", min_cell_frac_per_region)
print("min_cells_per_region:", min_cells_per_region)
print("raw regions:", raw_region_n)
print("kept regions:", kept_region_n)
print("kept region fraction:", f"{kept_region_frac:.2%}")

if kept_region_n == 0:
    raise ValueError(
        "No genomic regions passed filtering. "
        "Consider lowering min_cell_frac_per_region."
    )

kept_regions_file = result_dir / f"selected_regions_bin{bin_size}.tsv.gz"

region_qc = pd.DataFrame({
    "genomic_region_id": list(region_cell_count.keys()),
    "observed_cell_n": [region_cell_count[r] for r in region_cell_count.keys()],
})

region_qc["kept"] = region_qc["genomic_region_id"].isin(kept_regions)

region_qc = region_qc.sort_values(
    ["kept", "observed_cell_n"],
    ascending=[False, False]
)

region_qc.to_csv(
    kept_regions_file,
    sep="\t",
    index=False,
    compression="gzip"
)

print("\nsaved:", kept_regions_file)

print("\nPreview selected region QC:")
print(region_qc.head().to_string(index=False))

output

[Region filtering]
input cells: 613
min_cell_frac_per_region: 0.05
min_cells_per_region: 31
raw regions: 5401639
kept regions: 3855725
kept region fraction: 71.38%

saved: ../output/01_methyltree_input/WTJW969_CellAnno..._bin500_cov5pct_cell500/selected_regions_bin500.tsv.gz

Preview selected region QC:

genomic_region_id observed_cell_n kept
29 chr16:46399000-46399500 609 True
30 chr16:46399500-46400000 609 True
78795 chr16:46390000-46390500 609 True
78800 chr16:46394500-46395000 609 True
86112 chr17:26603500-26604000 609 True

python

# ============================================================
# 2.6.4 合并中间文件并过滤低覆盖细胞
# ============================================================

# 实际数据文件名带有 bin 信息；后面会额外创建标准入口 data.tsv.gz
bin_data_out = result_dir / f"bin{bin_size}_data.tsv.gz"

if bin_data_out.exists():
    bin_data_out.unlink()

first = True
kept_cell_ids = []
cell_region_count_kept = {}

combine_start = time.time()
n_files = len(per_cell_files)

print("\n[Combine per-cell files]")
print("per-cell files:", n_files)
print("kept regions:", len(kept_regions))
print("min_regions_per_cell:", min_regions_per_cell)
print_path("output", bin_data_out)

if n_files == 0:
    raise ValueError("No per-cell intermediate files found. Check previous binning step.")

for i, f in enumerate(per_cell_files, start=1):
    df_cell = pd.read_csv(f, sep="\t", compression="gzip")

    # 只保留通过 region 过滤的 genomic regions
    df_cell = df_cell[
        df_cell["genomic_region_id"].astype(str).isin(kept_regions)
    ].copy()

    # 如果该细胞过滤后没有 region，则从文件名中恢复 cell_id
    cell_id = (
        df_cell["cell_id"].iloc[0]
        if df_cell.shape[0] > 0
        else Path(f).name.split(f".{bin_size}bp.tsv.gz")[0]
    )

    cell_region_count_kept[cell_id] = df_cell.shape[0]

    # 通过 cell-level region 数过滤后，写入最终数据文件
    if df_cell.shape[0] >= min_regions_per_cell:
        kept_cell_ids.append(cell_id)

        df_cell.to_csv(
            bin_data_out,
            sep="\t",
            index=False,
            compression="gzip",
            mode="wt" if first else "at",
            header=first,
        )
        first = False

    if i % 100 == 0 or i == n_files:
        elapsed = time.time() - combine_start
        rate = i / elapsed if elapsed > 0 else 0
        remain = (n_files - i) / rate if rate > 0 else np.nan

        print(
            f"[{bin_size}bp] combined {i}/{n_files} cells | "
            f"elapsed {elapsed/60:.1f} min | ETA {remain/60:.1f} min",
            flush=True
        )

kept_cell_ids = set(map(str, kept_cell_ids))

print("\nFinished combining per-cell files.")
print("cells before min_regions_per_cell:", n_files)
print("kept cells after min_regions_per_cell:", len(kept_cell_ids))

if len(kept_cell_ids) == 0:
    raise ValueError(
        "No cells passed min_regions_per_cell. "
        "Consider lowering min_regions_per_cell or min_cell_frac_per_region."
    )

if not bin_data_out.exists() or bin_data_out.stat().st_size == 0:
    raise RuntimeError(f"Output data file was not created or is empty: {bin_data_out}")

print_path("saved", bin_data_out)

print("\nPreview kept region counts per cell:")
print(
    pd.Series(cell_region_count_kept, name="kept_region_n")
    .sort_values(ascending=False)
    .head()
    .to_string()
)

output

[Combine per-cell files]
per-cell files: 613
kept regions: 3855725
min_regions_per_cell: 500
output: ../output/01_methyltree_input/WTJW969_CellAnn...on.0.5_d30_0_bin500_cov5pct_cell500/bin500_data.tsv.gz
[500bp] combined 100/613 cells | elapsed 11.0 min | ETA 56.2 min
[500bp] combined 200/613 cells | elapsed 21.8 min | ETA 45.1 min
[500bp] combined 300/613 cells | elapsed 32.7 min | ETA 34.2 min
[500bp] combined 400/613 cells | elapsed 43.4 min | ETA 23.1 min
[500bp] combined 500/613 cells | elapsed 54.2 min | ETA 12.2 min
[500bp] combined 600/613 cells | elapsed 65.0 min | ETA 1.4 min
[500bp] combined 613/613 cells | elapsed 66.4 min | ETA 0.0 min

Finished combining per-cell files.
cells before min_regions_per_cell: 613
kept cells after min_regions_per_cell: 609
saved: ../output/01_methyltree_input/WTJW969_CellAnno...on.0.5_d30_0_bin500_cov5pct_cell500/bin500_data.tsv.gz

Preview kept region counts per cell:

TGTGGAAGTTATAGGTA 1685227
TGAAGGGTGAGGAATAA 1622389
AGTAATGTAGGTTGGGA 979046
AATTATGGTTATGGATA 961325
TGAGTGGATTATGATAG 958608
Name: kept_region_n, dtype: int64

python

# ============================================================
# 2.6.5 生成最终 sample_sheet、QC 文件和 summary
# ============================================================

# 只保留通过 min_regions_per_cell 过滤的细胞
sample_final = matched_sample_sheet[
    matched_sample_sheet["sample"].astype(str).isin(kept_cell_ids)
].copy()

if sample_final.shape[0] == 0:
    raise ValueError(
        "Final sample_sheet is empty. "
        "Check min_regions_per_cell and previous filtering steps."
    )

sample_final["Subset_key"] = out_label

sample_final_out = result_dir / "sample_sheet.tsv.gz"
sample_final.to_csv(
    sample_final_out,
    sep="\t",
    index=False,
    compression="gzip"
)


# ============================================================
# cell-level QC
# ============================================================

cell_qc = pd.DataFrame({
    "sample": list(cell_region_count_raw.keys()),
    "raw_region_n": [cell_region_count_raw[x] for x in cell_region_count_raw.keys()],
    "kept_region_n": [cell_region_count_kept.get(x, 0) for x in cell_region_count_raw.keys()],
    "pass_min_regions_per_cell": [str(x) in kept_cell_ids for x in cell_region_count_raw.keys()],
})

# 加入主要注释列，方便后续查看哪些分组被过滤
anno_cols = [
    "sample",
    "Clone_key",
    "Celltype_key",
    "Cluster_key",
    "CellAnnotation",
    "raw_Sample",
]

anno_cols = [x for x in anno_cols if x in matched_sample_sheet.columns]

cell_qc = cell_qc.merge(
    matched_sample_sheet[anno_cols].drop_duplicates(subset=["sample"]),
    on="sample",
    how="left"
)

cell_qc = cell_qc.sort_values(
    ["pass_min_regions_per_cell", "kept_region_n"],
    ascending=[False, False]
)

cell_qc_out = result_dir / "cell_region_qc.tsv"
cell_qc.to_csv(cell_qc_out, sep="\t", index=False)


# ============================================================
# 标准入口文件 data.tsv.gz
# ============================================================

# 实际数据文件为 bin{bin_size}_data.tsv.gz
# 这里额外创建标准入口 data.tsv.gz，后续主程序统一读取该文件
standard_data_out = result_dir / "data.tsv.gz"

if standard_data_out.is_symlink() or standard_data_out.exists():
    standard_data_out.unlink()

try:
    os.symlink(bin_data_out, standard_data_out)
    link_note = f"created symlink: {compact_path(standard_data_out)} -> {compact_path(os.readlink(standard_data_out))}"
except OSError:
    # 某些文件系统不支持软链接时，退化为复制
    shutil.copy2(bin_data_out, standard_data_out)
    link_note = f"copied data file: {bin_data_out} -> {standard_data_out}"


# ============================================================
# summary
# ============================================================

raw_region_n = int(len(region_cell_count))
kept_region_n = int(len(kept_regions))
kept_cells_n = int(len(kept_cell_ids))

summary = {
    "sample": sample_name,
    "subset": out_label,
    "bin_size": int(bin_size),
    "input_cells": int(n_cells_input),
    "kept_cells": kept_cells_n,
    "kept_cell_frac": float(kept_cells_n / n_cells_input) if n_cells_input > 0 else 0.0,
    "raw_region_n": raw_region_n,
    "kept_region_n": kept_region_n,
    "kept_region_frac": float(kept_region_n / raw_region_n) if raw_region_n > 0 else 0.0,
    "min_cell_frac_per_region": float(min_cell_frac_per_region),
    "min_cells_per_region": int(min_cells_per_region),
    "min_regions_per_cell": int(min_regions_per_cell),
    "threads": int(threads),
    "data_out": str(bin_data_out),
    "standard_data_out": str(standard_data_out),
    "sample_out": str(sample_final_out),
    "cell_qc_out": str(cell_qc_out),
    "selected_regions_out": str(kept_regions_file),
}

for col in ["Clone_key", "Celltype_key", "CellAnnotation", "Cluster_key", "raw_Sample"]:
    if col in sample_final.columns:
        summary[f"final_{col}_counts"] = (
            sample_final[col]
            .value_counts(dropna=False)
            .astype(int)
            .to_dict()
        )

summary_out = result_dir / "build_summary.json"

with open(summary_out, "w") as f:
    json.dump(summary, f, indent=2, ensure_ascii=False)


# ============================================================
# 输出检查
# ============================================================

print("\n[Final outputs]")
print("MethylTree required files:")
print_path("data.tsv.gz", standard_data_out)
print_path("sample_sheet.tsv.gz", sample_final_out)

print("\nRecommended QC files:")
print_path("cell_region_qc.tsv", cell_qc_out)
print_path("build_summary.json", summary_out)
print_path(f"selected_regions_bin{bin_size}.tsv.gz", kept_regions_file)

print("\nData link/copy:")
print(link_note)

print("\nFinal cells:", sample_final.shape[0])
print("Final regions:", kept_region_n)

print("\nFinal Celltype_key counts:")
print(sample_final["Celltype_key"].value_counts(dropna=False))

print("\nFinal Clone_key counts:")
print(sample_final["Clone_key"].value_counts(dropna=False).sort_index())

if "Cluster_key" in sample_final.columns:
    print("\nFinal Cluster_key counts:")
    print(sample_final["Cluster_key"].value_counts(dropna=False).sort_index())

print(f"\nTOTAL elapsed: {(time.time() - start_total) / 60:.1f} min")

output

[Final outputs]
MethylTree required files:
data.tsv.gz: ../output/01_methyltree_input/WTJW969_Ce...esolution.0.5_d30_0_bin500_cov5pct_cell500/data.tsv.gz
sample_sheet.tsv.gz: ../output/01_methyltree_input/WT...n.0.5_d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz

Recommended QC files:
cell_region_qc.tsv: ../output/01_methyltree_input/WTJ...on.0.5_d30_0_bin500_cov5pct_cell500/cell_region_qc.tsv
build_summary.json: ../output/01_methyltree_input/WTJ...on.0.5_d30_0_bin500_cov5pct_cell500/build_summary.json
selected_regions_bin500.tsv.gz: ../output/01_methyltr..._bin500_cov5pct_cell500/selected_regions_bin500.tsv.gz

Data link/copy:
created symlink: ../output/01_methyltree_input/WTJW96...on.0.5_d30_0_bin500_cov5pct_cell500/bin500_data.tsv.gz

Final cells: 609
Final regions: 3855725

Final Celltype_key counts:
Celltype_key
T_cells 609
Name: count, dtype: int64

Final Clone_key counts:
Clone_key
c0 609
Name: count, dtype: int64

Final Cluster_key counts:
Cluster_key
c0 609
Name: count, dtype: int64

TOTAL elapsed: 69.8 min

python

# ============================================================
# 2.6.6 最终检查：确认 MethylTree 输入已准备完成
# ============================================================

print("\nFinal output directory:")
print(compact_path(result_dir))

print("\nFiles:")
for f in [
    "data.tsv.gz",
    f"bin{bin_size}_data.tsv.gz",
    "sample_sheet.tsv.gz",
    "build_summary.json",
    "cell_region_qc.tsv",
    f"selected_regions_bin{bin_size}.tsv.gz",
]:
    p = result_dir / f
    size_gb = p.stat().st_size / 1024 / 1024 / 1024 if p.exists() and p.is_file() else 0
    print(f"{f}: {p.exists()} {size_gb:.3f} GB {compact_path(p)}")

sample_check = pd.read_csv(
    result_dir / "sample_sheet.tsv.gz",
    sep="\t",
    compression="gzip"
)

data_cell_ids = set()
for chunk in pd.read_csv(
    result_dir / "data.tsv.gz",
    sep="\t",
    compression="gzip",
    chunksize=5_000_000
):
    data_cell_ids.update(chunk["cell_id"].astype(str).unique())

print("\nsample_sheet cells:", sample_check["sample"].nunique())
print("data cells:", len(data_cell_ids))

assert set(sample_check["sample"].astype(str)) == data_cell_ids

print("\nMethylTree input is ready.")
print("\nMain program required files:")
print_path("DATA_FILE", result_dir / "data.tsv.gz")
print_path("SAMPLE_SHEET", result_dir / "sample_sheet.tsv.gz")
print("CLONE_KEY =", "Clone_key")

output

Final output directory:
../output/01_methyltree_input/WTJW969_CellAnnotation_T_cells__resolution.0.5_d30_0_bin500_cov5pct_cell500

Files:
data.tsv.gz True 1.631 GB ../output/01_methyltree_inp...esolution.0.5_d30_0_bin500_cov5pct_cell500/data.tsv.gz
bin500_data.tsv.gz True 1.631 GB ../output/01_methylt...on.0.5_d30_0_bin500_cov5pct_cell500/bin500_data.tsv.gz
sample_sheet.tsv.gz True 0.000 GB ../output/01_methyl...n.0.5_d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz
build_summary.json True 0.000 GB ../output/01_methylt...on.0.5_d30_0_bin500_cov5pct_cell500/build_summary.json
cell_region_qc.tsv True 0.000 GB ../output/01_methylt...on.0.5_d30_0_bin500_cov5pct_cell500/cell_region_qc.tsv
selected_regions_bin500.tsv.gz True 0.037 GB ../outpu..._bin500_cov5pct_cell500/selected_regions_bin500.tsv.gz

sample_sheet cells: 609
data cells: 609

MethylTree input is ready.

Main program required files:
DATA_FILE = ../output/01_methyltree_input/WTJW969_Cel...esolution.0.5_d30_0_bin500_cov5pct_cell500/data.tsv.gz
SAMPLE_SHEET = ../output/01_methyltree_input/WTJW969_...n.0.5_d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz
CLONE_KEY = Clone_key

2.7 输入构建结果说明

输入构建完成后，结果目录中主要包含以下文件：

文件	说明
`data.tsv.gz`	MethylTree 主程序使用的甲基化输入表
`sample_sheet.tsv.gz`	MethylTree 主程序使用的细胞注释表
`cell_region_qc.tsv`	每个细胞覆盖的 region 数及过滤结果
`build_summary.json`	本次构建参数、输入细胞数、保留细胞数和 region 数
`selected_regions_bin{bin_size}.tsv.gz`	每个 genomic bin 的覆盖细胞数及是否保留

其中，MethylTree 主程序真正需要的是 data.tsv.gz 和 sample_sheet.tsv.gz。

data.tsv.gz 为 long format，包含三列：cell_id、genomic_region_id 和 value。
sample_sheet.tsv.gz 至少需要包含 sample 和 Clone_key，也可以包含 Celltype_key、Cluster_key 等热图注释列。

3. 运行 MethylTree 主程序

本节使用 data.tsv.gz 和 sample_sheet.tsv.gz 作为输入，运行 MethylTree 主程序，完成细胞间相似性计算、谱系树构建、clone inference 和结果可视化。

如果已提前准备好符合格式要求的 MethylTree 输入文件，也可以从本节开始运行。

3.1 设置主程序输入目录

主程序输入目录需包含以下两个文件：

data.tsv.gz
sample_sheet.tsv.gz

如果已完成第 2 节的输入构建，可直接使用生成的 result_dir 作为输入目录；如果使用已有输入文件，请手动指定 INPUT_DIR。

python

# ============================================================
# 3.1 设置主程序输入目录
# ============================================================

from pathlib import Path
import pandas as pd
import json

# 如果已经运行过第 2 节，优先使用 result_dir
# 如果只运行主程序，请修改默认输入目录
try:
    INPUT_DIR = Path(result_dir)
except NameError:
    INPUT_DIR = Path(
        "../output/01_methyltree_input/"
        "WTJW969_CellAnnotation_T_cells__resolution.0.5_d30_0_bin500_cov5pct_cell500"
    )

DATA_FILE = INPUT_DIR / "data.tsv.gz"
SAMPLE_SHEET = INPUT_DIR / "sample_sheet.tsv.gz"
SUMMARY_FILE = INPUT_DIR / "build_summary.json"

def show_file(path):
    if path.exists() and path.is_file():
        return f"exists, {path.stat().st_size / 1024**3:.3f} GB"
    return "missing"

print_path("INPUT_DIR", INPUT_DIR)
print("DATA_FILE:", show_file(DATA_FILE))
print("SAMPLE_SHEET:", show_file(SAMPLE_SHEET))
print("SUMMARY_FILE:", show_file(SUMMARY_FILE))

if SAMPLE_SHEET.exists():
    sample_preview = pd.read_csv(SAMPLE_SHEET, sep="\t", compression="gzip")

    print("\nsample shape:")
    print(sample_preview.shape)

    if "Clone_key" in sample_preview.columns:
        print("\nClone_key counts:")
        print(sample_preview["Clone_key"].value_counts(dropna=False))

    if "Subset_key" in sample_preview.columns:
        print("\nSubset_key counts:")
        print(sample_preview["Subset_key"].value_counts(dropna=False))

    print("\nColumns:")
    for col in sample_preview.columns:
        print("  -", col)
else:
    print("\nPlease check INPUT_DIR before continuing.")

output

DATA_FILE: ../output/01_methyltree_input/WTJW969_Cell...d30_0_bin500_cov5pct_cell500/data.tsv.gz True 1.631 GB
SAMPLE_SHEET: ../output/01_methyltree_input/WTJW969_C..._d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz True
SUMMARY_FILE: ../output/01_methyltree_input/WTJW969_C...5_d30_0_bin500_cov5pct_cell500/build_summary.json True

sample shape:
(609, 15)

Clone_key counts:
Clone_key
c0 609
Name: count, dtype: int64

Subset_key counts:
Subset_key
T_cells_c0 609
Name: count, dtype: int64

resolution.0.5_d30 counts:
resolution.0.5_d30
0 609
Name: count, dtype: int64

tumor_program counts:
tumor_program
T_cells 609
Name: count, dtype: int64

Columns:
['sample', 'HQ', 'Clone_key', 'Celltype_key', 'tumor_...0.5_d30', 'gex_cb', 'm_cb', 'allc_path', 'Subset_key']

3.2 设置主程序参数

本步骤用于配置 MethylTree 主程序运行参数。运行前需根据输入目录、输出路径和分析目的确认以下参数：

INPUT_DIR：MethylTree 输入目录，需包含 data.tsv.gz 和 sample_sheet.tsv.gz；
OUT_ROOT：主程序结果输出目录；
CLONE_KEY：用于 clone 分组和展示的列，一般保持为 Clone_key；
HEATMAP_ADDITIONAL_KEY_LIST：热图侧边展示的额外注释列；
REMOVE_CELLTYPE_SIGNAL：是否去除细胞类型相关信号；
CLONE_INFERENCE_THRESHOLD：clone inference 使用的阈值。

单一细胞类型样本一般无需开启 REMOVE_CELLTYPE_SIGNAL；多细胞类型混合样本可根据是否需要削弱细胞类型主导信号决定是否开启。

python

# ============================================================
# 3.2 设置 MethylTree 主程序参数
# ============================================================

import os
import sys
from pathlib import Path

import numpy as np
import pandas as pd
import anndata as ad
import methyltree


# ============================================================
# 输入文件
# ============================================================

# 如果上一格已经设置 INPUT_DIR，这里会直接复用
# 如果直接从本格开始运行，则优先尝试使用第 2 节生成的 result_dir
if "INPUT_DIR" not in globals():
    try:
        INPUT_DIR = Path(result_dir)
    except NameError:
        INPUT_DIR = Path(
            "../output/01_methyltree_input/"
            "WTJW969_CellAnnotation_T_cells__resolution.0.5_d30_0_bin500_cov5pct_cell500"
        ).resolve()

DATA_FILE = INPUT_DIR / "data.tsv.gz"
SAMPLE_SHEET = INPUT_DIR / "sample_sheet.tsv.gz"
SUMMARY_FILE = INPUT_DIR / "build_summary.json"


# ============================================================
# 输出目录
# ============================================================

OUT_ROOT = Path("../output/03_methyltree_result")

DATA_DES = INPUT_DIR.name
CLONE_KEY = "Clone_key"

CLONE_INFERENCE_THRESHOLD = 0.6

threshold_tag = f"{CLONE_INFERENCE_THRESHOLD:.2f}".replace(".", "")
SAVE_DATA_DES = f"{DATA_DES}_{CLONE_KEY}_thr{threshold_tag}"
OUT_DIR = OUT_ROOT / SAVE_DATA_DES
FIGURE_PATH = OUT_DIR / "figure"

OUT_DIR.mkdir(parents=True, exist_ok=True)
FIGURE_PATH.mkdir(parents=True, exist_ok=True)


# ============================================================
# MethylTree 分析参数
# ============================================================

COMPUTE_SIMILARITY = True
SIMILARITY_CORRECTION = True
UPDATE_SAMPLE_INFO = True

REMOVE_CELLTYPE_SIGNAL = False
CELL_TYPE_KEY = None

OPTIMIZE_TREE = False

PERFORM_CLONE_INFERENCE = True
PERFORM_COARSE_GRAINING = False
PERFORM_MEMORY_ANALYSIS = False
PERFORM_DEPTH_ANALYSIS = False

SAVE_ADATA = True


# ============================================================
# 热图参数
# ============================================================

HEATMAP_VMAX_PERCENTILE = 99.9
HEATMAP_VMIN_PERCENTILE = 60
HEATMAP_FIGSIZE = (10, 9.5)
HEATMAP_SHOW_LABEL = False
HEATMAP_SHOW_LEGEND = True
HEATMAP_FONTSIZE = 8

HEATMAP_ADDITIONAL_KEY_LIST = [
    "Cluster_key",
    "CellAnnotation",
]


# ============================================================
# 基础检查
# ============================================================

if not DATA_FILE.exists():
    raise FileNotFoundError(f"DATA_FILE not found: {DATA_FILE}")

if not SAMPLE_SHEET.exists():
    raise FileNotFoundError(f"SAMPLE_SHEET not found: {SAMPLE_SHEET}")

print("python:", sys.executable)
print("methyltree:", methyltree.__file__)
print_path("INPUT_DIR", INPUT_DIR)
print(f"DATA_FILE: {compact_path(DATA_FILE)} {DATA_FILE.stat().st_size / 1024**3:.3f} GB")
print_path("SAMPLE_SHEET", SAMPLE_SHEET)
print_path("OUT_DIR", OUT_DIR)
print("CLONE_INFERENCE_THRESHOLD:", CLONE_INFERENCE_THRESHOLD)

output

python: python3
methyltree: methyltree/__init__.py
INPUT_DIR: ../output/01_methyltree_input/WTJW969_Cell...n_T_cells__resolution.0.5_d30_0_bin500_cov5pct_cell500
DATA_FILE: ../output/01_methyltree_input/WTJW969_Cell....0.5_d30_0_bin500_cov5pct_cell500/data.tsv.gz 1.631 GB
SAMPLE_SHEET: ../output/01_methyltree_input/WTJW969_C...n.0.5_d30_0_bin500_cov5pct_cell500/sample_sheet.tsv.gz
OUT_DIR: ../output/03_methyltree_result/WTJW969_CellA...tion.0.5_d30_0_bin500_cov5pct_cell500_Clone_key_thr060
CLONE_INFERENCE_THRESHOLD: 0.6

python

# ============================================================
# 3.2.1 读取 sample_sheet 和 build_summary
# ============================================================

sample = pd.read_csv(SAMPLE_SHEET, sep="\t", compression="gzip")
sample["sample"] = sample["sample"].astype(str)

required_cols = [CLONE_KEY] + HEATMAP_ADDITIONAL_KEY_LIST
missing_cols = [x for x in required_cols if x not in sample.columns]

if len(missing_cols) > 0:
    raise ValueError(f"sample_sheet missing required columns: {missing_cols}")

print("sample shape:", sample.shape)
print("sample cells:", sample["sample"].nunique())

print("\nClone_key counts:")
print(sample[CLONE_KEY].value_counts(dropna=False).sort_index())

print("\nHeatmap annotation columns:")
print(HEATMAP_ADDITIONAL_KEY_LIST)

if SUMMARY_FILE.exists():
    with open(SUMMARY_FILE) as f:
        summary = json.load(f)

    print("\nbuild_summary:")
    print("kept_cells:", summary.get("kept_cells"))
    print("kept_region_n:", summary.get("kept_region_n"))
else:
    summary = None

output

sample shape: (609, 15)
sample cells: 609

Clone_key counts:
Clone_key
c0 609
Name: count, dtype: int64

Heatmap annotation columns:
['Cluster_key', 'CellAnnotation']

build_summary:
kept_cells: 609
kept_region_n: 3855725

python

# ============================================================
# 3.2.2 可选：检查 data.tsv.gz 与 sample_sheet.tsv.gz 的细胞是否一致
# ============================================================

RUN_INPUT_CELL_CHECK = False

if RUN_INPUT_CELL_CHECK:
    data_cell_ids = set()

    for chunk in pd.read_csv(
        DATA_FILE,
        sep="\t",
        compression="gzip",
        usecols=["cell_id"],
        chunksize=5_000_000
    ):
        data_cell_ids.update(chunk["cell_id"].astype(str).unique())

    sample_cell_ids = set(sample["sample"].astype(str))

    print("data cells:", len(data_cell_ids))
    print("sample cells:", len(sample_cell_ids))
    print("sample cells missing in data:", len(sample_cell_ids - data_cell_ids))
    print("data cells missing in sample:", len(data_cell_ids - sample_cell_ids))

    assert sample_cell_ids == data_cell_ids

    print("input cell check passed")
else:
    print("Skip input cell check. Set RUN_INPUT_CELL_CHECK = True if needed.")

output

Skip input cell check. Set RUN_INPUT_CELL_CHECK = True if needed.

3.3 构建 AnnData 并运行 MethylTree

本步骤将 data.tsv.gz 转换为 cell × genomic region 甲基化矩阵，并运行 MethylTree 主程序。对于细胞数或 region 数较大的数据，矩阵构建可能占用较多内存，建议先使用小规模数据测试流程。

python

# ============================================================
# 3.3.1 读取 data.tsv.gz 并构建 AnnData
# ============================================================

import gc

print("[INFO] Reading data.tsv.gz...")

df = pd.read_csv(
    DATA_FILE,
    sep="\t",
    compression="gzip",
    usecols=["cell_id", "genomic_region_id", "value"]
)

df["cell_id"] = df["cell_id"].astype(str)
df["genomic_region_id"] = df["genomic_region_id"].astype(str)
df["value"] = df["value"].astype("float32")

print("long table shape:", df.shape)
print("data cells:", df["cell_id"].nunique())
print("data regions:", df["genomic_region_id"].nunique())


print("\n[INFO] Pivot to cell x region matrix...")

mat_df = df.pivot(
    index="cell_id",
    columns="genomic_region_id",
    values="value"
).astype("float32")

del df
gc.collect()


# 按 sample_sheet 顺序排列细胞
sample_order = sample["sample"].astype(str).tolist()
mat_df = mat_df.reindex(sample_order)


# 检查 sample_sheet 中的细胞是否都有甲基化数据
missing_cells = mat_df.index[mat_df.isna().all(axis=1)].tolist()

if len(missing_cells) > 0:
    raise ValueError(
        f"{len(missing_cells)} cells in sample_sheet have no methylation data. "
        f"Preview: {missing_cells[:10]}"
    )


mat_values = mat_df.to_numpy(dtype=np.float32)

print("matrix shape:", mat_df.shape)
print("missing value fraction:", np.isnan(mat_values).mean())


print("\n[INFO] Building AnnData...")

adata = ad.AnnData(X=mat_values)
adata.obs_names = mat_df.index.astype(str)
adata.var_names = mat_df.columns.astype(str)
adata.obs = sample.set_index("sample").loc[adata.obs_names].copy()


# 检查必要注释列
required_obs_cols = [CLONE_KEY] + HEATMAP_ADDITIONAL_KEY_LIST
missing_obs_cols = [x for x in required_obs_cols if x not in adata.obs.columns]

if len(missing_obs_cols) > 0:
    raise ValueError(f"Required columns not found in adata.obs: {missing_obs_cols}")


print(f"AnnData object: {adata.n_obs} cells x {adata.n_vars} regions")
print("obs columns:")
for col in adata.obs.columns:
    print("  -", col)

print("\nClone_key counts:")
print(adata.obs[CLONE_KEY].value_counts(dropna=False).sort_index())

output

[INFO] Reading data.tsv.gz...n long table shape: (276326052, 3)
data cells: 609
data regions: 3855725

[INFO] Pivot to cell x region matrix...n matrix shape: (609, 3855725)
missing value fraction: 0.8823211303695384

[INFO] Building AnnData...n AnnData object with n_obs × n_vars = 609 × 3855725
obs: 'HQ', 'Clone_key', 'Celltype_key', 'tumor_program...ution.0.5_d30', 'gex_cb', 'm_cb', 'allc_path', 'Subset_key'

Clone_key counts:
Clone_key
c0 609
Name: count, dtype: int64

python

# ============================================================
# 3.3.2 运行 MethylTree 主程序
# ============================================================

from contextlib import redirect_stdout
from io import StringIO

_methyltree_log = StringIO()
with redirect_stdout(_methyltree_log):
    adata_final, my_tree = methyltree.analysis.comprehensive_lineage_analysis(
        out_dir=str(OUT_DIR),
        data_path=str(INPUT_DIR),
        save_data_des=SAVE_DATA_DES,
        clone_key=CLONE_KEY,
        adata_orig=adata,
    
        compute_similarity=COMPUTE_SIMILARITY,
        update_sample_info=UPDATE_SAMPLE_INFO,
        similarity_correction=SIMILARITY_CORRECTION,
    
        remove_celltype_signal=REMOVE_CELLTYPE_SIGNAL,
        cell_type_key=CELL_TYPE_KEY,
    
        optimize_tree=OPTIMIZE_TREE,
    
        heatmap_vmax_percentile=HEATMAP_VMAX_PERCENTILE,
        heatmap_vmin_percentile=HEATMAP_VMIN_PERCENTILE,
        heatmap_figsize=HEATMAP_FIGSIZE,
        heatmap_show_label=HEATMAP_SHOW_LABEL,
        heatmap_show_legend=HEATMAP_SHOW_LEGEND,
        heatmap_additional_key_list=HEATMAP_ADDITIONAL_KEY_LIST,
        heatmap_fontsize=HEATMAP_FONTSIZE,
    
        fig_dir=str(FIGURE_PATH),
        data_des=DATA_DES,
    
        perform_coarse_graining=PERFORM_COARSE_GRAINING,
    
        perform_clone_inference=PERFORM_CLONE_INFERENCE,
        clone_inference_threshold=CLONE_INFERENCE_THRESHOLD,
    
        # 当前示例中 Clone_key 全部为 c0，因此不计算 accuracy
        perform_accuracy_estimation=False,
    
        perform_memory_analysis=PERFORM_MEMORY_ANALYSIS,
        perform_depth_analysis=PERFORM_DEPTH_ANALYSIS,
    
        save_adata=SAVE_ADATA,
    )

print_compact_log(_methyltree_log.getvalue())

print("MethylTree analysis finished.")
print_path("OUT_DIR", OUT_DIR)
print_path("FIGURE_PATH", FIGURE_PATH)

if "inferred_clones" in adata_final.obs.columns:
    print("\ninferred_clones:")
    print(adata_final.obs["inferred_clones"].value_counts(dropna=False))

if "filtered_clone" in adata_final.obs.columns:
    print("\nfiltered_clone:")
    print(adata_final.obs["filtered_clone"].value_counts(dropna=False))

output

use provided adata
adata shape: (609, 3855725)
update sample
X_similarity_correlation_raw not found in adata.obsm
re-compute similarity matrix
Use correlation for similarity
duration: 269.16495990753174
correct similarity: outer loop 0; current epsilon 0.05
Use fast/analytical correction method

1%| | 8/1000 [00:00<00:06, 160.51it/s]

std: 0.020
similarity normalize----

Reconstruction method: UPGMA



00%|██████████| 100/100 [00:48<00:00,  2.07it/s]

Inferred clone number: 7
save data at 
../output/03_methyltree_result/WTJW969_CellAnnotation...0.5_d30_0_bin500_cov5pct_cell500_Clone_key_thr060.h5ad
MethylTree analysis finished.
OUT_DIR: ../output/03_methyltree_result/WTJW969_CellA...tion.0.5_d30_0_bin500_cov5pct_cell500_Clone_key_thr060
FIGURE_PATH: ../output/03_methyltree_result/WTJW969_C...5_d30_0_bin500_cov5pct_cell500_Clone_key_thr060/figure

inferred_clones:
inferred_clones
NaN        571
clone_0     27
clone_2      3
clone_1      2
clone_3      2
clone_4      2
clone_5      2
Name: count, dtype: int64

filtered_clone:
filtered_clone
NaN    571
c0      38
Name: count, dtype: int64

3.4 保存 clone 结果并重新绘图

本步骤用于查看 MethylTree 的 clone inference 结果，并将每个细胞对应的 inferred clone 信息保存为表格文件。随后，脚本会将 inferred clone 作为新的注释列加入热图侧边注释，并重新绘制结果图。

python

# ============================================================
# 3.4.1 查看 clone inference 结果
# ============================================================

print("Clone_key:")
print(adata_final.obs["Clone_key"].value_counts(dropna=False).sort_index())

if "inferred_clones" in adata_final.obs.columns:
    print("\ninferred_clones:")
    print(adata_final.obs["inferred_clones"].value_counts(dropna=False))

if "filtered_clone" in adata_final.obs.columns:
    print("\nfiltered_clone:")
    print(adata_final.obs["filtered_clone"].value_counts(dropna=False))

output

Clone_key:
Clone_key
c0 609
Name: count, dtype: int64

inferred_clones:
inferred_clones
NaN 571
clone_0 27
clone_2 3
clone_1 2
clone_3 2
clone_4 2
clone_5 2
Name: count, dtype: int64

filtered_clone:
filtered_clone
NaN 571
c0 38
Name: count, dtype: int64

python

# ============================================================
# 3.4.2 保存 clone inference 结果
# ============================================================

if "inferred_clones" not in adata_final.obs.columns:
    raise ValueError("inferred_clones not found.")

# 1. clone 汇总表
vc = adata_final.obs["inferred_clones"].value_counts(dropna=True)

clone_summary = pd.DataFrame([{
    "total_cells": adata_final.n_obs,
    "assigned_cells": int(adata_final.obs["inferred_clones"].notna().sum()),
    "assigned_fraction": float(adata_final.obs["inferred_clones"].notna().mean()),
    "clone_n": int(adata_final.obs["inferred_clones"].dropna().nunique()),
    "size_distribution": str(vc.to_dict()),
}])

clone_summary_out = OUT_DIR / "clone_inference_summary.tsv"
clone_summary.to_csv(clone_summary_out, sep="\t", index=False)


# 2. 细胞级 clone 结果表
assignment_cols = [
    "Clone_key",
    "Celltype_key",
    "Cluster_key",
    "CellAnnotation",
    "inferred_clones",
]
assignment_cols = [x for x in assignment_cols if x in adata_final.obs.columns]

clone_assignment = adata_final.obs[assignment_cols].copy()
clone_assignment.insert(0, "cell_id", adata_final.obs_names.astype(str))

clone_assignment_out = OUT_DIR / "clone_assignment_by_cell.tsv"
clone_assignment.to_csv(clone_assignment_out, sep="\t", index=False)

print("saved:", clone_summary_out)
print("saved:", clone_assignment_out)

print(clone_summary.to_string(index=False))

output

saved: ../output/03_methyltree_result/WTJW969_CellAnn...t_cell500_Clone_key_thr060/clone_inference_summary.tsv
saved: ../output/03_methyltree_result/WTJW969_CellAnn..._cell500_Clone_key_thr060/clone_assignment_by_cell.tsv

total_cells assigned_cells assigned_fraction clone_n \\
0 609 38 0.062397 6

size_distribution
0 {'clone_0': 27, 'clone_2': 3, 'clone_1': 2, 'c...

python

# ============================================================
# 3.4.3 生成热图展示用 clone 注释列
# ============================================================

if "inferred_clones" not in adata_final.obs.columns:
    raise ValueError("inferred_clones not found.")

adata_final.obs["inferred_clone_for_plot"] = adata_final.obs["inferred_clones"].astype("object")

print("inferred_clone_for_plot:")
print(adata_final.obs["inferred_clone_for_plot"].value_counts(dropna=False))

output

inferred_clone_for_plot:
inferred_clone_for_plot
NaN 571
clone_0 27
clone_2 3
clone_4 2
clone_5 2
clone_1 2
clone_3 2
Name: count, dtype: int64

python

# ============================================================
# 3.4.4 重新绘制热图：展示细胞注释和 inferred clone
# ============================================================

PLOT_OUT_DIR = OUT_ROOT / f"{SAVE_DATA_DES}_plot_clone_annotation"
PLOT_FIGURE_PATH = PLOT_OUT_DIR / "figure"

PLOT_OUT_DIR.mkdir(parents=True, exist_ok=True)
PLOT_FIGURE_PATH.mkdir(parents=True, exist_ok=True)

from contextlib import redirect_stdout
from io import StringIO

_replot_log = StringIO()
with redirect_stdout(_replot_log):
    adata_plot, tree_plot = methyltree.analysis.comprehensive_lineage_analysis(
        out_dir=str(PLOT_OUT_DIR),
        data_path=str(INPUT_DIR),
        save_data_des=f"{SAVE_DATA_DES}_plot_clone_annotation",
        clone_key=CLONE_KEY,
        adata_orig=adata_final.copy(),
    
        # 复用已有 similarity，只重新画图
        compute_similarity=False,
        update_sample_info=False,
        similarity_correction=False,
        similarity_normalize=False,
    
        remove_celltype_signal=False,
        cell_type_key=None,
        optimize_tree=False,
    
        heatmap_vmax_percentile=HEATMAP_VMAX_PERCENTILE,
        heatmap_vmin_percentile=HEATMAP_VMIN_PERCENTILE,
        heatmap_figsize=HEATMAP_FIGSIZE,
        heatmap_show_label=HEATMAP_SHOW_LABEL,
        heatmap_show_legend=True,
        heatmap_additional_key_list=[
            "Celltype_key",
            "inferred_clone_for_plot",
        ],
        heatmap_fontsize=HEATMAP_FONTSIZE,
    
        fig_dir=str(PLOT_FIGURE_PATH),
        data_des=f"{DATA_DES}_clone_annotation",
    
        perform_coarse_graining=False,
        perform_clone_inference=False,
        perform_accuracy_estimation=False,
        perform_memory_analysis=False,
        perform_depth_analysis=False,
    
        save_adata=True,
    )

print_compact_log(_replot_log.getvalue())

print("Re-plot finished.")
print_path("PLOT_FIGURE_PATH", PLOT_FIGURE_PATH)

output

use provided adata
adata shape: (609, 3855725)
Warning: 'Celltype_key' is not a valid key. Force to update df_sample
update sample

Reconstruction method: UPGMA
save data at 
../output/03_methyltree_result/WTJW969_CellAnnotation...ct_cell500_Clone_key_thr060_plot_clone_annotation.h5ad
Re-plot finished.
PLOT_FIGURE_PATH: ../output/03_methyltree_result/WTJW..._cell500_Clone_key_thr060_plot_clone_annotation/figure

python

# ============================================================
# 3.4.5 查看输出文件
# ============================================================

print("Clone result files:")
for f in ["clone_inference_summary.tsv", "clone_assignment_by_cell.tsv"]:
    print(f"- {f}: {compact_path(OUT_DIR / f)}")

print("\nFigure files:")
for p in sorted(PLOT_FIGURE_PATH.glob("*")):
    print(f"- {p.name}: {compact_path(p)}")

output

Clone result files:
../output/03_methyltree_result/WTJW969_CellAnnotation...t_cell500_Clone_key_thr060/clone_inference_summary.tsv
../output/03_methyltree_result/WTJW969_CellAnnotation..._cell500_Clone_key_thr060/clone_assignment_by_cell.tsv

Figure files:
../output/03_methyltree_result/WTJW969_CellAnnotation...ne_annotation_similarity_matrix_multiple_colorbars.pdf

最终输出包括 clone inference 汇总表、细胞级 clone assignment 表，以及带细胞注释和 inferred clone 注释条的相似性热图。

4. 主程序结果解读

MethylTree 主程序运行结束后，重点查看 OUT_DIR 和 FIGURE_PATH 中的输出结果。本教程当前主要输出包括相似性热图、clone inference 汇总表、细胞级 clone assignment 表以及保存后的 AnnData 文件。

4.1 相似性热图

相似性热图是最主要的可视化结果，用于观察细胞之间的甲基化相似性结构。

重点看：

是否存在明显的相似性 block；
block 是否与已有注释标签对应；
是否存在明显离群细胞。

如果热图中出现清晰 block，说明部分细胞在甲基化信号上更接近。
如果热图整体较散，说明当前数据中的谱系信号可能较弱，也可能与覆盖度、缺失值或过滤参数有关。

4.2 树重建结果

MethylTree 主程序会基于细胞间相似性进行树重建，例如日志中出现：

text

Reconstruction method: UPGMA

说明程序已经进行了树重建。

需要注意，本教程当前没有额外导出单独的树结构图文件。树重建结果主要用于细胞排序、热图展示和后续 clone inference。

因此在当前教程中，建议主要结合相似性热图和 clone inference 结果进行解释。

如果后续需要单独保存或展示树结构，可以进一步导出 my_tree 或根据 MethylTree 返回的 tree 对象单独作图。

4.3 clone inference 结果

如果开启：

python

perform_clone_inference = True

主程序会根据相似性结构和树重建结果推断潜在 clone 或细胞模块。

建议重点查看：

text

clone_inference_summary.tsv
clone_assignment_by_cell.tsv

其中：

文件	说明
`clone_inference_summary.tsv`	汇总推断 clone 数量、被分配细胞数和 clone 大小分布
`clone_assignment_by_cell.tsv`	记录每个细胞对应的 inferred clone，可用于后续合并到 RNA 注释表

如果 inferred clone 与热图 block 基本一致，说明该结果相对更可信。
如果大量细胞没有被分配到 clone，说明当前阈值下可稳定识别的 clone 较少，不一定代表程序运行失败。

4.4 结果解释边界

MethylTree 结果需要结合输入标签来源进行解释。

如果使用的是真实 clone 标签，可以用于评估 clone 重建效果。
如果使用的是 RNA cluster、细胞类型或其他注释标签，则只能说明甲基化谱系结构与这些标签之间的一致性，不能直接解释为真实 clone 准确率。

5. 常见问题

5.1 输入细胞数和细胞类型怎么控制？

MethylTree 可以同时分析多个相关细胞类型，尤其适用于比较不同细胞类型之间的谱系关系或分化关系。但不建议一开始把所有细胞无选择地全部放入分析。

运行时间主要受输入细胞数、保留 region 数和矩阵缺失情况影响。细胞数越多、region 数越多，输入构建和主程序都会更慢。

实际使用时，建议根据研究问题提取重点细胞群，例如某一类细胞、几个相关细胞类型，或候选 RNA cluster。流程测试阶段可以先用较小子集跑通，再扩大分析范围。

5.2 最终保留细胞很少

优先查看 cell_region_qc.tsv 和 build_summary.json。如果大量细胞的 kept_region_n 低于 min_regions_per_cell，通常说明单细胞甲基化覆盖度偏低，或过滤阈值偏严格。

可以适当降低：

python

min_regions_per_cell = 100

或：

python

min_cell_frac_per_region = 0.005

但需要注意，降低阈值只能增加进入分析的细胞数，不能改善原始数据质量。如果覆盖度本身不足，后续热图和 clone inference 结果仍可能不稳定。

5.3 region 数量过多、运行太慢

如果 data.tsv.gz 文件很大，或主程序构建 cell × region 矩阵很慢，通常是输入细胞数或保留 region 数过多。

优先建议缩小分析范围，只保留与当前问题相关的细胞群。也可以适当提高 region 过滤阈值：

python

min_cell_frac_per_region = 0.02

或：

python

min_cell_frac_per_region = 0.05

阈值越高，保留 region 越少，运行越快；但过高可能丢失谱系信号，需要结合 build_summary.json 中的保留细胞数和 region 数判断。

5.4 热图结构不明显

如果相似性热图没有明显 block，不一定是程序运行失败。常见原因包括覆盖度不足、缺失值较多、输入细胞过杂，或该细胞群本身谱系差异较弱。

这时不建议直接反复调整主程序参数。应先检查输入质量、保留细胞数、保留 region 数，并确认当前输入细胞是否符合分析目的。如果测序深度不足，即使放宽过滤阈值，也不一定能得到稳定结构。

5.5 `missing allc` 很多

missing allc 通常是 barcode 匹配或路径设置问题。优先检查 barcode_col、remove_barcode_suffix 和 allc_dir。

如果 metadata 中的 barcode 带有 _1、_2 等后缀，而 allc 文件名不带后缀，可以尝试：

python

remove_barcode_suffix = True

同时确认 *_allc.gz 文件名去掉 _allc.gz 后，能够与 metadata 中用于匹配的 barcode 对应。

更多功能：MethylTree 的更多高级可视化功能请详见 MethylTree 官方文档。

甲基化 + RNA 双组学：MethylTree 谱系分析 ​

1. 模块简介 ​

2. 构建 MethylTree 输入文件 ​

2.1 准备输入数据 ​

2.2 allc 文件格式 ​

2.3 导入依赖包 ​

2.4 设置输入构建参数 ​

2.5 检查输入文件并匹配 barcode ​

2.6 生成 MethylTree 输入文件 ​

2.7 输入构建结果说明 ​

3. 运行 MethylTree 主程序 ​

3.1 设置主程序输入目录 ​

3.2 设置主程序参数 ​

3.3 构建 AnnData 并运行 MethylTree ​

3.4 保存 clone 结果并重新绘图 ​

4. 主程序结果解读 ​

4.1 相似性热图 ​

4.2 树重建结果 ​

4.3 clone inference 结果 ​

4.4 结果解释边界 ​

5. 常见问题 ​

5.1 输入细胞数和细胞类型怎么控制？ ​

5.2 最终保留细胞很少 ​

5.3 region 数量过多、运行太慢 ​

5.4 热图结构不明显 ​

5.5 missing allc 很多 ​

甲基化 + RNA 双组学：MethylTree 谱系分析

1. 模块简介

2. 构建 MethylTree 输入文件

2.1 准备输入数据

2.2 allc 文件格式

2.3 导入依赖包

2.4 设置输入构建参数

2.5 检查输入文件并匹配 barcode

2.6 生成 MethylTree 输入文件

2.7 输入构建结果说明

3. 运行 MethylTree 主程序

3.1 设置主程序输入目录

3.2 设置主程序参数

3.3 构建 AnnData 并运行 MethylTree

3.4 保存 clone 结果并重新绘图

4. 主程序结果解读

4.1 相似性热图

4.2 树重建结果

4.3 clone inference 结果

4.4 结果解释边界

5. 常见问题

5.1 输入细胞数和细胞类型怎么控制？

5.2 最终保留细胞很少

5.3 region 数量过多、运行太慢

5.4 热图结构不明显

5.5 `missing allc` 很多