Bio Fragment Analysis

Version Compatibility

Reference examples tested with: numpy 1.26+, pandas 2.2+, pysam 0.22+, finaletoolkit 0.7+

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: FinaleToolkit CLI subcommands are hyphenated (frag-length-bins, end-motifs, delfi-gc-correct); the Python functions are underscored in finaletoolkit.frag. The filter subcommand is filter-file, NOT filter-bam. Griffin is a Snakemake pipeline, not an importable function. DELFI is a methodology and a company (DELFI Diagnostics), not a pip install-able tool.

Fragment Analysis

"Analyze cfDNA fragment patterns for cancer signal" -> Quantify nucleosome-footprint features (size ratios, protection scores, accessibility profiles, end motifs) from plasma WGS for detection or tissue-of-origin.

• Python/CLI: finaletoolkit for DELFI ratios, WPS, end-motifs, MDS, cleavage profiles
• Snakemake: Griffin pipeline for GC-corrected nucleosome profiling at TF/accessible sites
• Python: pysam for a custom binned short/long ratio when a dependency-light readout is wanted

The Single Most Important Modern Insight -- Fragmentomics measures nucleosome positioning, not sequence; the biology is real but it can be destroyed in the wet lab or left un-GC-corrected

Plasma cfDNA is the digestion product of chromatin by apoptotic and intracellular nucleases, so every fragmentomic feature is a re-readout of the same physical object: the nucleosome footprint. The ~167 bp mode is the 147 bp histone-protected core plus ~20 bp of linker; the 10.4 bp sawtooth below it is the helical pitch of DNA on the histone surface (the nuclease cuts only where the minor groove faces out, once per turn). DELFI ratios, WPS, Griffin profiles, end-motifs, and OCF are four views of this one object, not four independent measurements -- their correlation inflates apparent multi-feature performance and leaks across train/test splits.

Two consequences dominate practice. First, the single biggest threat to any fragmentomic feature is GC, library, and batch confounding, not biology: an uncorrected genome-wide short/long ratio tracks GC content and library prep far more strongly than tumor fraction, which is why naive fragmentomics works in discovery and dies in validation. Griffin's actual contribution is its fragment-length-specific GC correction, not the nucleosome plot. Second, DELFI-style ratios are partly entangled with copy-number alteration and coverage (a bin's ratio reflects both fragmentation state and how many genomes contributed): deconvolving the fragmentation-specific signal from CNA is a known open problem, so a raw 5 Mb ratio is a hybrid CNA + fragmentation + GC readout, not pure fragmentation.

Methods / Feature-Family Landscape

DELFI = DNA EvaLuation of Fragments; the end-motif/MDS biology is anchored in DNASE1L3, whose deletion reorders length and end-motif frequencies (Serpas 2019 PNAS 116:641). Methodology here is still evolving (CNA deconvolution, standalone GC correctors like GCparagon) -- verify current best practice against live FinaleToolkit and Griffin docs before committing to one feature family.

Decision Tree by Scenario

Genome-Wide Short/Long Ratio (custom DELFI-style)

Goal: Produce a genome-wide vector of short-to-long fragment ratios in fixed bins as a dependency-light DELFI-style feature, with the explicit caveat that without GC correction it is a GC + CNA readout.

Approach: Walk proper-pair fragments per bin from the BAM template length, classify each as short (100-150 bp) or long (151-220 bp), and emit the per-bin ratio. For a publication-grade profile, prefer FinaleToolkit delfi (which GC-corrects) over this illustrative version.

import pysam
import numpy as np
import pandas as pd

def binned_short_long_ratio(bam_path, bin_size=5_000_000, chroms=None):
    '''Per-bin short(100-150)/long(151-220) ratio. NOT GC-corrected -- illustrative only.'''
    chroms = chroms or [f'chr{i}' for i in range(1, 23)]
    bam = pysam.AlignmentFile(bam_path, 'rb')
    rows = []
    for chrom in chroms:
        if chrom not in bam.references:
            continue
        n_bins = bam.get_reference_length(chrom) // bin_size + 1
        short = np.zeros(n_bins)
        long = np.zeros(n_bins)
        for read in bam.fetch(chrom):
            if not read.is_proper_pair or read.is_secondary or read.template_length <= 0:
                continue
            size = read.template_length
            b = read.reference_start // bin_size
            if 100 <= size <= 150:
                short[b] += 1
            elif 151 <= size <= 220:
                long[b] += 1
        ratio = np.divide(short, long, out=np.full(n_bins, np.nan), where=long > 0)
        rows.extend({'chrom': chrom, 'bin': i, 'short': short[i], 'long': long[i], 'ratio': ratio[i]} for i in range(n_bins))
    bam.close()
    return pd.DataFrame(rows)

GC-Corrected DELFI Score (FinaleToolkit)

Goal: Compute a GC-corrected DELFI score so the genome-wide profile reflects fragmentation rather than base composition.

Approach: FinaleToolkit's delfi corrects short and long bin counts for GC before forming the ratio; the CLI and Python API are equivalent. Run on a BAM/CRAM or a tabix-indexed .frag.gz fragment file.

# CLI (subcommands are hyphenated). delfi positionals: input chrom_sizes reference bins_file.
# GC correction is ON by default (-G disables it); 100kb bins are merged to 5Mb by default.
# -R keeps no-coverage regions when the genome is not hg19.
finaletoolkit delfi sample.bam hg38.chrom.sizes hg38.fa bins_100kb.bed -g gaps.bed -R -o sample.delfi.bed
finaletoolkit end-motifs sample.bam hg38.fa -o sample.end_motifs.tsv
finaletoolkit mds sample.end_motifs.tsv          # Motif Diversity Score (normalized Shannon entropy)
finaletoolkit wps sample.bam sites.bed -c hg38.chrom.sizes -o sites.wps.bw   # per-site, not a single region

from finaletoolkit.frag import delfi, end_motifs, wps  # public finaletoolkit.frag symbols
# delfi() returns GC-corrected short/long per bin; end_motifs() returns an EndMotifFreqs
# object whose .motif_diversity_score() gives the MDS (there is no top-level frag.mds).

Griffin Nucleosome Profiling (Snakemake pipeline)

Goal: Obtain GC-corrected composite coverage around a TF/accessible-site set for tissue-of-origin, robust at low tumor fraction and ~0.1x WGS.

Approach: Griffin is not an importable function; it is three sequential Snakemake modules. Run them in order against samples.yaml, the hg38 reference, and a sites.yaml site list.

# Run each module from Griffin's snakemakes/ dir (config edited per cohort)
snakemake -s griffin_genome_GC_frequency/griffin_genome_GC_frequency.snakefile --cores 8
snakemake -s griffin_GC_and_mappability_correction/griffin_GC_and_mappability_correction.snakefile --cores 8
snakemake -s griffin_nucleosome_profiling/griffin_nucleosome_profiling.snakefile --cores 8
# Output: GC-corrected + uncorrected composite coverage profiles around each site set.

Per-Method Failure Modes

Uncorrected GC dominates the ratio

Trigger: comparing raw short/long ratios across samples without GC correction. Mechanism: PCR and binding-based purification overrepresent GC-balanced fragments, and short vs long fragments have different GC dependence. Symptom: a beautiful discovery-cohort separation that collapses in validation; the profile clusters by sequencing batch. Fix: GC-correct (FinaleToolkit delfi/delfi-gc-correct, Griffin, or GCparagon) and co-process all samples through one pipeline.

Cross-protocol non-comparability

Trigger: combining ssDNA and dsDNA libraries, or two end-repair/PCR chemistries, in one analysis. Mechanism: ssDNA prep recovers the sub-100 bp ultrashort population that dsDNA prep loses at the double-strand ligation step, shifting the entire size distribution and every derived feature. Symptom: a model trained on one chemistry mislabels the other systematically. Fix: a fragmentomic model is conditioned on its library chemistry -- match protocols, never cross them, and state the prep as a precondition.