Bio Methylation Based Detection

Version Compatibility

Reference examples tested with: MethylDackel 0.6+, Bismark 0.24+, numpy 1.26+, pandas 2.2+, scipy 1.12+, statsmodels 0.14+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Notes specific to this skill: MethylDackel extract bedGraph column order is fixed (chrom / start / end / methylation-% rounded to integer / count-methylated / count-unmethylated); always run MethylDackel mbias first and feed its suggested --OT/--OB trimming into extract. cfMeDIP data are coverage, not conversion — do not feed them into per-CpG bisulfite pipelines.

Methylation-Based Detection

"Detect cancer and find where it came from using cfDNA methylation" -> Call per-CpG (and read-level) methylation from converted plasma DNA, then deconvolve tissue-of-origin against a reference atlas.

• CLI: MethylDackel extract for per-CpG methylation from bisulfite/EM-seq BAMs
• CLI: MethylDackel mbias to choose strand-specific trimming before extraction
• R: MEDIPS / QSEA for cfMeDIP enrichment (coverage, not conversion)
• Python: scipy.optimize.nnls for atlas-based tissue deconvolution

The Single Most Important Modern Insight -- methylation beats mutations, and read-level haplotypes beat averaged beta

Methylation is the right altitude for multi-cancer early detection (MCED) and tissue-of-origin (TOO) in a single assay because the genome carries thousands of stable, cell-type-specific differentially methylated regions, whereas somatic mutations are sparse, recurrent only at a few driver loci, and carry no tissue label. On the same cfDNA inside CCGA, the methylation assay outperformed WGS-SNV/CNV approaches, which is why GRAIL down-selected to a targeted methylation panel (Liu 2020). One panel answers both "is there cancer?" and "where is it?" — mutations answer neither well.

The sensitivity engine is read-level, not site-level. Averaging beta across reads at a CpG discards phasing. A single tumor-derived fragment that is concordantly methylated across the k CpGs of a methylation haplotype block (Guo 2017) has a background probability of roughly p^k of arising from the hematopoietic ocean; for a block of 5-8 CpGs that is small enough that ONE such fragment is strong evidence, independent of tumor fraction. Per-site beta dilutes that signal into sampling noise and clonal-hematopoiesis variance and has essentially no power at parts-per-million tumor fraction. The correct primitive for detection is molecule counting over haplotype blocks, not site averaging.

Conversion Chemistry Tradeoffs

The conversion step is chosen on the worst possible substrate — already-fragmented, low-input plasma DNA — so destructiveness is load-bearing, not a footnote.

Bisulfite is gold standard for cell-line gDNA, not for low-input fragmented plasma; EM-seq and TAPS exist precisely to recover ctDNA molecules bisulfite destroys. Neither bisulfite nor EM-seq separates 5mC from 5hmC without added oxBS/TAB steps; TAPS variants (TAPSbeta, CAPS) can split the marks. cfMeDIP coverage is enrichment, not quantitation — density bias must be modeled before any absolute-methylation claim.

Decision Tree by Scenario

Methodology evolves; verify current atlas versions and panel-marker coverage against live tool docs before committing — atlas markers are platform-specific and do not transfer across assays.

Extract Per-CpG Methylation with MethylDackel

Goal: Produce per-CpG methylation calls from a bisulfite/EM-seq cfDNA BAM, with end-repair artifacts trimmed.

Approach: Run mbias first to read the suggested strand-specific trimming, then extract with --mergeContext and that trimming so each CpG is one row; parse the fixed 6-column bedGraph.

# Step 1: choose trimming. mbias prints a suggestion like --OT 2,0,0,98 and writes M-bias SVGs.
MethylDackel mbias ref.fa sample.bam sample_mbias

# Step 2: extract per-CpG (one row per CpG) with the suggested trimming.
MethylDackel extract ref.fa sample.bam -o sample --mergeContext --minDepth 1 --OT 2,0,0,98
# Output sample_CpG.bedGraph columns (fixed order):
#   chrom  start  end  methylation%(integer, rounded)  count_methylated  count_unmethylated

Add --CHG --CHH only to audit non-CpG methylation (a conversion-failure check); CpG is the default context. Low per-CpG depth gates can erase cfDNA signal — prefer region/molecule aggregation over a high --minDepth.

Tissue-of-Origin Deconvolution Against an Atlas

Goal: Attribute cfDNA to its cell types of origin and surface a solid-tissue coefficient elevated above the hematopoietic baseline.

Approach: Model the observed methylation vector m ~ A*w with atlas A (rows = markers, cols = cell types), solve for non-negative mixing fractions w with NNLS over atlas-covered markers, then renormalize so the fractions sum to one.

from scipy.optimize import nnls

def deconvolve_tissue(sample_beta, atlas):
    'sample_beta: Series indexed by marker; atlas: DataFrame markers x cell_types.'
    markers = sample_beta.index.intersection(atlas.index)
    w, _ = nnls(atlas.loc[markers].values, sample_beta.loc[markers].values)
    w = w / w.sum()
    return dict(zip(atlas.columns, w))

The simplex constraint (w >= 0, sum w = 1) is mandatory — unconstrained regression gives nonsense fractions (Moss 2018). Use Loyfer 2023's fragment-level WGBS atlas (39 cell types from 205 healthy samples) where the assay covers its markers; a generic atlas does not transfer, because WGBS-fragment markers differ from 450K/EPIC probes and from capture-panel coverage.

Region-Level DMR Discovery

Goal: Define a discriminating panel by finding regions (not single CpGs) that separate cancer from normal cfDNA.

Approach: Aggregate per-CpG beta into pre-defined regions/blocks, test cancer vs normal per region, and control FDR with Benjamini-Hochberg specified explicitly — naming method='fdr_bh' rather than relying on the statsmodels default ('hs', Holm-Sidak).

from scipy import stats
from statsmodels.stats.multitest import multipletests

def region_dmrs(cancer, normal, region_col='region'):
    'cancer/normal: long DataFrames with [region_col, beta]; one row per sample-region.'
    out = []
    for region, c in cancer.groupby(region_col)['beta']:
        n = normal.loc[normal[region_col] == region, 'beta'].dropna()
        c = c.dropna()
        if len(c) < 3 or len(n) < 3:
            continue
        _, p = stats.mannwhitneyu(c, n, alternative='two-sided')
        out.append((region, c.mean() - n.mean(), p))
    import pandas as pd
    res = pd.DataFrame(out, columns=['region', 'delta_beta', 'pvalue'])
    res['fdr'] = multipletests(res['pvalue'], method='fdr_bh')[1]
    return res.sort_values('fdr')

cfMeDIP Enrichment Analysis

Goal: Get density-corrected differential methylation from antibody-enrichment coverage rather than conversion data.

Approach: Use MEDIPS (Lienhard 2014) for CpG-density-corrected differential coverage, or QSEA (Lienhard 2017) when absolute, BS-comparable methylation levels are needed — QSEA adds a Bayesian CpG-density model, CNV correction, and TMM effective-library-size normalization. Both are R/Bioconductor; do not pass cfMeDIP coverage through MethylDackel.