scMethyl + RNA Multi-omics: MethylTree Lineage Analysis

Author: SeekGene

Time: 60 min

Words: 11.8k words

Updated: 2026-07-15

Reads: 0 times

1. Module Introduction

This module is developed based on MethylTree. It is used to extract lineage-related signals from single-cell DNA methylation data and combine them with RNA annotation information to analyze clone structure, differentiation relationships, and state heterogeneity within cell populations.

Conventional single-cell RNA analysis mainly reflects the current transcriptional state of cells, but it is difficult to directly determine whether these cells originate from the same ancestral cell or share a common clonal origin. MethylTree uses gradually accumulated epimutation signals in DNA methylation as a natural lineage record to infer lineage relationships among cells. Compared with methods that rely on artificial genetic barcodes or somatic mutations, this type of method is better suited for non-invasive lineage analysis in human samples or datasets without a pre-labeling system.

This tutorial is organized according to the actual execution order:

1. Build MethylTree input files: convert per-cell *_allc.gz files and a cell annotation table into data.tsv.gz and sample_sheet.tsv.gz;
2. Run the MethylTree main program: perform similarity calculation, lineage tree reconstruction, heatmap visualization, and clone inference;
3. Review and interpret the output results.

Biologically, MethylTree can help address the following questions:

• Whether cells share a common lineage origin: identify cell groups with similar lineage histories by calculating methylation similarity between cells.
• Whether clone or subclone structures exist within a cell population: detect potential methyl clones or lineage-related cell modules based on methylation similarity and clone inference.
• Whether differentiation relationships or cell fate choices have a lineage basis: in development, differentiation, or tumor evolution studies, assess whether different cell states share common ancestors, branching relationships, or clonal expansion patterns.

Before You Start

The performance of MethylTree strongly depends on the quality of single-cell methylation data. According to the software authors, single-cell genome coverage is preferably no less than approximately 2% in practice. When coverage is too low, the resulting cell × genomic region matrix will contain many missing values, making cell-cell similarity calculation unstable and reducing the reliability of the similarity heatmap, lineage tree, and clone inference.

This tutorial aggregates CpG methylation signals into 500 bp genomic bins by default. The bin size can also be adjusted to 1 kb depending on data quality. Both 500 bp and 1 kb are commonly used and suitable window sizes for MethylTree analysis, balancing reduced sparsity with retention of local epimutation signals. Larger windows can further reduce missing values, but may average out local methylation differences and weaken lineage signal detection.

2. Build MethylTree Input Files

This section converts per-cell methylation files and a cell annotation table into input files that can be directly used by the MethylTree main program. The following files will be generated:

• data.tsv.gz: methylation matrix input file;
• sample_sheet.tsv.gz: cell annotation and grouping information file.

Subsequent MethylTree lineage reconstruction and visualization are based on these two files.

---

2.1 Prepare Input Data

Before running the analysis, prepare the following files:

• Per-cell methylation files: *_allc.gz, with one allc file for each cell. These files are usually obtained by extracting allcools.tar.gz from the standard analysis results.
• Cell annotation table: meta.tsv, which must contain at least one barcode column that can be matched to allc filenames. If cell type, clone, or cluster information needs to be shown, the table may also include columns such as CellAnnotation, clone_id, and resolution.0.5_d30.

After removing the _allc.gz suffix, each allc filename should match the barcode in meta.tsv. For example, AAAGAAGAAGGATTGTT_allc.gz corresponds to the cell ID AAAGAAGAAGGATTGTT.

If RNA barcodes contain sample suffixes, such as AAAGAAGAAGGATTGTT_7, while the allc filenames do not contain _7, set remove_barcode_suffix = True in the parameter section. The script will remove the suffix after the last underscore to match allc filenames.

---

2.2 allc File Format

This tutorial assumes that the input allc files use the standard format, where each row corresponds to one methylation site:

chrom    pos    strand    context    mc    cov    methylated

The fields are defined as follows:

Field	Description
chrom	Chromosome name
pos	Genomic coordinate of the site
strand	Strand direction
context	Methylation context, such as CG, CHG, or CHH
mc	UMI/read count supporting methylation at this site
cov	UMI/read coverage at this site
methylated	Indicator of significant methylation; usually 1 if no significance test is performed

The downstream calculation mainly uses the chrom, pos, context, mc, and cov fields. strand and methylated are retained in the original allc files but are not used for bin-level methylation rate calculation.

This tutorial only uses CpG methylation signals, retaining sites whose context starts with CG.

The script aggregates CpG sites into fixed-size genomic bins according to the specified bin_size, and calculates the methylation rate for each cell in each bin:

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

Therefore, the value column in the final data.tsv.gz represents the CpG methylation rate of a given cell in a given genomic bin, typically ranging from 0–100.

2.3 Import Dependencies

# ============================================================
# 2.3 Import dependencies
# ============================================================

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 Set Input-Building Parameters

This step configures the key parameters for building MethylTree input files, including input paths, barcode matching rules, annotation columns, and cell/region filtering thresholds. Modify the parameters according to the actual dataset before running the workflow.

# ============================================================
# 2.4 Set input-building parameters
# ============================================================

# Sample name, used for output directory naming
sample_name = "WTJW969"

# Input files
# meta_file: cell annotation table
# allc_dir: directory containing per-cell *_allc.gz files; searched recursively by the script
meta_file = "../input/meta.tsv"
allc_dir = "../input/allcools"

# Output directory
# This section generates data.tsv.gz, sample_sheet.tsv.gz, and QC files in this directory
outdir = "../output/01_methyltree_input"


# ============================================================
# Helper functions for document-friendly outputs
# ============================================================

def compact_path(path, max_len=100):
    """Return a shortened one-line path for notebook/document display."""
    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):
    """Print captured logs while shortening very long path-like lines."""
    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 column configuration
# ============================================================

# Barcode column in meta.tsv used to match allc filenames
barcode_col = "barcode"

# Set to True if barcodes in meta.tsv contain suffixes such as _1 or _2 but allc filenames do not
remove_barcode_suffix = False


# ============================================================
# Cell filtering configuration
# ============================================================

# None means analyzing all cells
# A dictionary filters cells by specified columns, for example only cluster 0 within T cells
subset_filters = {
    "CellAnnotation": ["T_cells"],
    "resolution.0.5_d30": [0]
}


# ============================================================
# MethylTree grouping label configuration
# ============================================================

# clone_key_col: column used to generate Clone_key in sample_sheet.tsv.gz
# If using RNA clusters, set clone_key_prefix = "c" to generate labels such as c0, c1, and c2
# If using true clone labels or cell types, clone_key_prefix can be set to None
clone_key_col = "resolution.0.5_d30"
clone_key_prefix = "c"

# celltype_key_col: column used to generate Celltype_key
celltype_key_col = "CellAnnotation"

# cluster_key_col: column used to generate Cluster_key; set to None if not needed
cluster_key_col = "resolution.0.5_d30"


# ============================================================
# Bin aggregation and filtering parameters
# ============================================================

# bin_size: window size for aggregating CpG methylation signals; commonly 500 or 1000
bin_size = 500

# min_cell_frac_per_region: filtering threshold for the fraction of cells covered per region
# 0.05 means that a region must be covered in at least 5% of input cells to be retained
min_cell_frac_per_region = 0.05

# min_regions_per_cell: cell filtering threshold
# A cell must cover at least this number of retained regions to enter the final input
min_regions_per_cell = 500

# Number of parallel threads
threads = 32

2.5 Check Input Files and Match Barcodes

This step checks whether the input files are complete and matches cell barcodes in meta.tsv to the corresponding *_allc.gz files. The script scans the allc directory, reads the cell annotation table, generates the initial sample_sheet, and filters cells without matched allc files.

If many cells are reported as missing allc files, check whether barcode_col, remove_barcode_suffix, and allc_dir are set correctly.

# ============================================================
# 2.5.1 Path normalization and output directory creation
# ============================================================

def clean_label(x):
    """Convert any label into a string suitable for directory names."""
    x = str(x)
    x = x.replace("/", "-").replace("\\", "-")
    x = re.sub(r"\s+", "_", x)
    return x


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


# Path normalization
meta_file = Path(meta_file).expanduser()
allc_dir = Path(allc_dir).expanduser()
outdir = Path(outdir).expanduser()


# Generate output labels according to cell filtering conditions
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)


# Automatically generate a cov label from min_cell_frac_per_region
# For example, 0.01 -> cov1pct, 0.05 -> cov5pct, and 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: final output directory for the MethylTree input files in this run
result_dir = outdir / f"{sample_name}_{out_label}"

# tmp_dir: directory for per-cell bin-level intermediate files
# If rerun with unchanged parameters, existing per-cell intermediate files will be reused
# To rerun from scratch, manually delete this directory
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)


# Parameter checks
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

# ============================================================
# 2.5.2 Check allc files
# ============================================================

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)}")

# Extract cell IDs from allc filenames for matching to barcodes in meta.tsv
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

# ============================================================
# 2.5.3 Read metadata and generate the initial sample_sheet
# ============================================================

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

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

# Check required columns
required_cols = [barcode_col, clone_key_col, celltype_key_col]
if cluster_key_col is not None:
    required_cols.append(cluster_key_col)

# Also check columns used in 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())


# ============================================================
# Filter cells according to 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])


# ============================================================
# Generate barcodes used to match allc files
# ============================================================

# gex_cb: keep the original barcode
sub_meta["gex_cb"] = sub_meta[barcode_col].astype(str)

# m_cb: barcode used to match allc filenames
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)


# ============================================================
# Generate MethylTree annotation columns
# ============================================================

# 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)

    # If the original value already has a prefix such as c0/c1, do not add it again
    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)


# ============================================================
# Generate the initial 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")


# ============================================================
# Output check
# ============================================================

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

# ============================================================
# 2.5.4 Scan allc files and match sample_sheet
# ============================================================

def get_cell_id_from_allc_path(p):
    """Extract cell ID from a *_allc.gz filename."""
    name = Path(p).name
    return name.replace("_allc.gz", "")


# Scan allc files
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],
})

# Check duplicated allc cell IDs
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())

# If multiple allc files exist for the same cell ID, keep the first one by default
allc_df = allc_df.drop_duplicates(subset=["sample"], keep="first").copy()

# Match with 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."
    )

# Save the matched sample_sheet
sample_out = result_dir / "sample_sheet.tsv.gz"
matched_sample_sheet.to_csv(sample_out, sep="\t", index=False, compression="gzip")


# ============================================================
# Output matching results
# ============================================================

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 Generate MethylTree Input Files

This step converts per-cell *_allc.gz files into MethylTree input format, including CpG site selection, bin-level methylation rate calculation, and filtering of low-coverage regions and cells.

Existing intermediate files will be reused automatically to avoid repeated computation.

# ============================================================
# 2.6.1 Define the function for converting allc to 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,
    }

# ============================================================
# 2.6.2 Generate bin-level intermediate files per cell
# ============================================================

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

input cells: 613
threads: 32
workers used: 32
tmp_dir: ../output/01_methyltree_input/WTJW969_CellAn...30_0_bin500_cov5pct_cell500/intermediate_500bp_by_cell
[500bp] binned 50/613 cells | elapsed 0.2 min | speed 3.98 cells/s | ETA 2.4 min
[500bp] binned 100/613 cells | elapsed 0.4 min | speed 4.03 cells/s | ETA 2.1 min
[500bp] binned 150/613 cells | elapsed 0.6 min | speed 4.07 cells/s | ETA 1.9 min
[500bp] binned 200/613 cells | elapsed 0.8 min | speed 4.12 cells/s | ETA 1.7 min
[500bp] binned 250/613 cells | elapsed 1.0 min | speed 4.10 cells/s | ETA 1.5 min
[500bp] binned 300/613 cells | elapsed 1.2 min | speed 4.06 cells/s | ETA 1.3 min
[500bp] binned 350/613 cells | elapsed 1.5 min | speed 4.02 cells/s | ETA 1.1 min
[500bp] binned 400/613 cells | elapsed 1.7 min | speed 4.00 cells/s | ETA 0.9 min
[500bp] binned 450/613 cells | elapsed 1.9 min | speed 4.01 cells/s | ETA 0.7 min
[500bp] binned 500/613 cells | elapsed 2.1 min | speed 4.00 cells/s | ETA 0.5 min
[500bp] binned 550/613 cells | elapsed 2.3 min | speed 3.99 cells/s | ETA 0.3 min
[500bp] binned 600/613 cells | elapsed 2.5 min | speed 3.98 cells/s | ETA 0.1 min
[500bp] binned 613/613 cells | elapsed 2.6 min | speed 3.95 cells/s | ETA 0.0 min

Finished per-cell binning.
per-cell intermediate files: 613
raw regions before filtering: 5401639

Preview raw region counts per cell:



TGTGGAAGTTATAGGTA 1881198
TGAAGGGTGAGGAATAA 1797073
AGTAATGTAGGTTGGGA 1083892
AATTATGGTTATGGATA 1065599
TGAGTGGATTATGATAG 1063672
Name: raw_region_n, dtype: int64

# ============================================================
# 2.6.3 Filter 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

# ============================================================
# 2.6.4 Merge intermediate files and filter low-coverage cells
# ============================================================

# The actual data file contains bin information in its name; a standard entry file data.tsv.gz will be created later
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")

    # Keep only genomic regions that passed region filtering
    df_cell = df_cell[
        df_cell["genomic_region_id"].astype(str).isin(kept_regions)
    ].copy()

    # If no regions remain for this cell after filtering, recover cell_id from the filename
    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]

    # After filtering by cell-level region count, write to the final data file
    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

# ============================================================
# 2.6.5 Generate final sample_sheet, QC files, and summary
# ============================================================

# Keep only cells that passed the min_regions_per_cell filter
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()],
})

# Add major annotation columns to facilitate checking which groups were filtered
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)


# ============================================================
# Standard entry file data.tsv.gz
# ============================================================

# The actual data file is bin{bin_size}_data.tsv.gz
# Create an additional standard entry file data.tsv.gz for the main program
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:
    # If symlinks are not supported by the filesystem, fall back to copying
    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)


# ============================================================
# Output check
# ============================================================

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

# ============================================================
# 2.6.6 Final check: confirm that MethylTree input is ready
# ============================================================

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 Input-Building Output Files

After input construction, the result directory mainly contains the following files:

File	Description
data.tsv.gz	Methylation input table used by the MethylTree main program
sample_sheet.tsv.gz	Cell annotation table used by the MethylTree main program
cell_region_qc.tsv	Number of covered regions per cell and cell filtering results
build_summary.json	Parameters, input cell count, retained cell count, and retained region count for this run
selected_regions_bin{bin_size}.tsv.gz	Coverage cell count and retention status for each genomic bin

The MethylTree main program only requires data.tsv.gz and sample_sheet.tsv.gz.

data.tsv.gz is a long-format table with three columns: cell_id, genomic_region_id, and value.
sample_sheet.tsv.gz must contain at least sample and Clone_key, and may also include heatmap annotation columns such as Celltype_key and Cluster_key.

3. Run the MethylTree Main Program

This section uses data.tsv.gz and sample_sheet.tsv.gz as inputs to run the MethylTree main program, including cell-cell similarity calculation, lineage tree reconstruction, clone inference, and result visualization.

If MethylTree-compatible input files have already been prepared, the workflow can also be started from this section.

3.1 Set the Main-Program Input Directory

The main-program input directory must contain the following two files:

• data.tsv.gz
• sample_sheet.tsv.gz

If Section 2 has been completed, the generated result_dir can be used directly as the input directory. If using existing input files, manually specify INPUT_DIR.

# ============================================================
# 3.1 Set the main-program input directory
# ============================================================

from pathlib import Path
import pandas as pd
import json

# If Section 2 has been run, use result_dir first
# If running only the main program, modify the default input directory
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 Set Main-Program Parameters

This step configures the MethylTree main-program parameters. Before running, confirm the input directory, output path, and analysis-related parameters:

• INPUT_DIR: MethylTree input directory, which must contain data.tsv.gz and sample_sheet.tsv.gz;
• OUT_ROOT: output directory for the main-program results;
• CLONE_KEY: column used for clone grouping and display, usually kept as Clone_key;
• HEATMAP_ADDITIONAL_KEY_LIST: additional annotation columns shown next to the heatmap;
• REMOVE_CELLTYPE_SIGNAL: whether to remove cell type-associated signals;
• CLONE_INFERENCE_THRESHOLD: threshold used for clone inference.

For samples containing a single cell type, REMOVE_CELLTYPE_SIGNAL is usually not needed. For mixed cell-type samples, decide whether to enable it depending on whether cell type-driven signals should be reduced.

# ============================================================
# 3.2 Set MethylTree main-program parameters
# ============================================================

import os
import sys
from pathlib import Path

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


# ============================================================
# Input files
# ============================================================

# If INPUT_DIR has been set in the previous cell, it will be reused here
# If starting from this cell, try to use result_dir generated in Section 2 first
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"


# ============================================================
# Output directory
# ============================================================

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 analysis parameters
# ============================================================

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 parameters
# ============================================================

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",
]


# ============================================================
# Basic checks
# ============================================================

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

# ============================================================
# 3.2.1 Read sample_sheet and 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

# ============================================================
# 3.2.2 Optional: check consistency between cells in data.tsv.gz and 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 Build AnnData and Run MethylTree

This step converts data.tsv.gz into a cell × genomic region methylation matrix and runs the MethylTree main program. For datasets with many cells or regions, matrix construction may require substantial memory. It is recommended to test the workflow on a smaller subset first.

# ============================================================
# 3.3.1 Read data.tsv.gz and build 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()


# Order cells according to sample_sheet
sample_order = sample["sample"].astype(str).tolist()
mat_df = mat_df.reindex(sample_order)


# Check whether all cells in sample_sheet have methylation data
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()


# Check required annotation columns
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

# ============================================================
# 3.3.2 Run the MethylTree main program
# ============================================================

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,
    
        # In this example, all Clone_key values are c0, so accuracy is not calculated
        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 Save Clone Results and Replot the Heatmap

This step reviews the MethylTree clone inference results and saves the inferred clone assignment for each cell as a table. The inferred clone information is then added as a new side annotation column for the heatmap, and the result figure is replotted.

# ============================================================
# 3.4.1 Review clone inference results
# ============================================================

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

# ============================================================
# 3.4.2 Save clone inference results
# ============================================================

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

# 1. Clone summary table
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. Cell-level clone assignment table
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...

# ============================================================
# 3.4.3 Generate clone annotation column for heatmap display
# ============================================================

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

# ============================================================
# 3.4.4 Replot heatmap: show cell annotations and inferred clones
# ============================================================

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(),
    
        # Reuse existing similarity and only replot the figure
        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

# ============================================================
# 3.4.5 Review output files
# ============================================================

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

The final outputs include a clone inference summary table, a cell-level clone assignment table, and similarity heatmaps with cell annotations and inferred clone annotation bars.

4. Interpret Main-Program Results

After the MethylTree main program finishes, focus on the output files in OUT_DIR and FIGURE_PATH. The main outputs in this tutorial include similarity heatmaps, the clone inference summary table, the cell-level clone assignment table, and the saved AnnData object.

4.1 Similarity Heatmap

The similarity heatmap is the most important visualization result and is used to inspect the methylation similarity structure among cells.

Focus on:

1. Whether clear similarity blocks are present;
2. Whether the blocks correspond to existing annotation labels;
3. Whether obvious outlier cells are present.

If clear blocks are present in the heatmap, some cells are more similar in methylation signals.
If the heatmap appears diffuse, the lineage signal in the current dataset may be weak, or the result may be affected by low coverage, missing values, or filtering parameters.

4.2 Tree Reconstruction Result

The MethylTree main program performs tree reconstruction based on cell-cell similarity. For example, the log may contain:

Reconstruction method: UPGMA

This indicates that tree reconstruction has been performed.

Note that this tutorial does not export a separate tree figure. The reconstructed tree is mainly used for cell ordering, heatmap visualization, and downstream clone inference.

Therefore, in this tutorial, interpretation should mainly rely on the similarity heatmap and clone inference results.

If a standalone tree structure needs to be saved or displayed in the future, my_tree or the tree object returned by MethylTree can be exported and plotted separately.

4.3 Clone Inference Result

If the following option is enabled:

perform_clone_inference = True

the main program will infer potential clones or cell modules based on the similarity structure and reconstructed tree.

Check the following files:

clone_inference_summary.tsv
clone_assignment_by_cell.tsv

File	Description
clone_inference_summary.tsv	Summary of inferred clone number, assigned cell count, and clone size distribution
clone_assignment_by_cell.tsv	Inferred clone assignment for each cell; can be merged with RNA annotations for downstream analysis

If inferred clones are largely consistent with heatmap blocks, the result is more likely to be reliable.
If many cells are not assigned to any clone, it means that few stable clones can be identified under the current threshold; this does not necessarily indicate that the program failed.

4.4 Interpretation Boundaries

MethylTree results should be interpreted together with the source of the input labels.

If true clone labels are used, the results can be used to evaluate clone reconstruction.
If RNA clusters, cell types, or other annotation labels are used, the results only indicate consistency between methylation lineage structure and those labels, and should not be directly interpreted as true clone reconstruction accuracy.

5. Frequently Asked Questions

5.1 How should input cell number and cell types be controlled?

MethylTree can analyze multiple related cell types simultaneously and is particularly useful for comparing lineage or differentiation relationships across cell types. However, it is not recommended to include all cells without selection at the beginning.

Runtime is mainly affected by the number of input cells, the number of retained regions, and matrix missingness. More cells and more regions will make both input construction and the main program slower.

In practice, select cell populations according to the research question, such as one specific cell type, several related cell types, or candidate RNA clusters. During workflow testing, first use a small subset to confirm that the workflow runs successfully, then expand the analysis scope.

5.2 Very Few Cells Are Retained

First check cell_region_qc.tsv and build_summary.json. If many cells have kept_region_n below min_regions_per_cell, the single-cell methylation coverage may be low, or the filtering threshold may be too strict.

You can lower:

min_regions_per_cell = 100

or:

min_cell_frac_per_region = 0.005

However, lowering thresholds only increases the number of cells entering the analysis and does not improve the raw data quality. If the coverage itself is insufficient, downstream heatmap and clone inference results may still be unstable.

5.3 Too Many Regions or Slow Runtime

If data.tsv.gz is large, or building the cell × region matrix is slow, the input may contain too many cells or too many retained regions.

The first recommendation is to reduce the analysis scope and retain only cell populations relevant to the current question. You can also increase the region filtering threshold:

min_cell_frac_per_region = 0.02

or:

min_cell_frac_per_region = 0.05

Higher thresholds retain fewer regions and speed up the workflow, but overly strict thresholds may lose lineage signals. Check the retained cell and region counts in build_summary.json when adjusting these parameters.

5.4 Heatmap Structure Is Not Clear

If the similarity heatmap does not show clear blocks, it does not necessarily mean that the program failed. Common reasons include insufficient coverage, many missing values, overly heterogeneous input cells, or weak lineage differences within the analyzed population.

In this case, do not repeatedly tune main-program parameters first. Instead, check input quality, retained cell count, retained region count, and whether the selected cells match the analysis goal. If sequencing depth is insufficient, relaxing filtering thresholds may still not produce a stable structure.

5.5 Many missing allc Records

Many missing allc records usually indicate barcode matching or path configuration issues. First check barcode_col, remove_barcode_suffix, and allc_dir.

If barcodes in metadata contain suffixes such as _1 or _2, while allc filenames do not, try:

remove_barcode_suffix = True

Also make sure that each *_allc.gz filename, after removing _allc.gz, can be matched to the barcode column used in metadata.

More Features: For additional advanced visualization functions in MethylTree, refer to the MethylTree official documentation.