bio-fragment-analysis

What this skill does

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

| Feature family | Primary method | What it physically measures | Citation |
|---|---|---|---|
| Genome-wide short/long ratio | DELFI: ~5 Mb bins, GC-corrected short(100-150)/long(151-220), boosted classifier | Coarse fragmentation state across the genome (entangled with CNA + GC) | Cristiano 2019 Nature 570:385 |
| Nucleosome positioning | WPS: spanning fragments minus end-containing fragments in a sliding window | Where nucleosomes sit; promoter/gene-body phasing encodes tissue + expression | Snyder 2016 Cell 164:57 |
| GC-corrected accessibility | Griffin: length-specific GC correction then composite coverage around site sets | TF/DHS accessibility in the tissue of origin; robust at low tumor fraction | Doebley 2022 Nat Commun 13:7475 |
| End motifs / MDS | 4-mer at the 5' cut site; MDS = normalized Shannon entropy of the 256 motifs | Nuclease-cleavage signature (DNASE1L3 sculpts the normal CC-ending spectrum) | Jiang 2020 Cancer Discov 10:664 |
| Orientation-aware ends | OCF: phase offset between upstream- and downstream-end peaks in open chromatin | Tissue-of-origin via end orientation, not coverage | Sun 2019 Genome Res 29:418 |

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

| Scenario | Recommended | Why |
|---|---|---|
| Cancer detection, pan-genome screen | DELFI genome-wide short/long profile (FinaleToolkit `delfi`) | Coarse genome-wide signal is what the boosted classifier was built on; cheap at low-pass WGS |
| Tissue / subtype of origin (e.g. ER status, NEPC) | Griffin nucleosome profiling around TF/accessible site sets | Accessibility composite scales with the contributing tissue; the GC correction makes it portable |
| Low tumor fraction (TF < ~0.03) | Griffin (GC-correction removes dominant technical signal) and/or in-silico size selection (90-150 bp) | Raw ratios are GC-dominated at low TF; Griffin holds AUC ~0.92 vs ~0.99 at TF >= 0.05 |
| Nuclease / cleavage biology, single-scalar comparison | End-motifs + MDS (FinaleToolkit `end-motifs` then `mds`) | MDS is one number per sample; rises in cancer as orderly DNASE1L3 cleavage is lost |
| Nucleosome positions / TF footprints directly | WPS (FinaleToolkit `wps` + `adjust-wps`) | WPS peaks recover nucleosome positions; S-WPS exposes TF footprints |
| Any cross-batch or cross-protocol comparison | GC-correct AND co-process through one pipeline, or do not compare | Uncorrected, cross-protocol fragmentomics is uninterpretable (see Failure Modes) |

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

```python
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