$npx -y skills add GPTomics/bioSkills --skill atac-qcATAC-seq library quality control -- TSS enrichment, FRiP, fragment-size periodicity, library complexity (NRF/PBC1/PBC2), mitochondrial fraction, and ENCODE 4 thresholds. Use when assessing whether an ATAC-seq library passes ENCODE acceptance criteria, diagnosing transposition art
| 1 | ## Version Compatibility |
| 2 | |
| 3 | Reference examples tested with: deepTools 3.5+, Picard 3.1+, samtools 1.19+, bedtools 2.31+, ATACseqQC 1.26+, pysam 0.22+, pyBigWig 0.3+, numpy 1.26+, pandas 2.2+, MultiQC 1.21+. |
| 4 | |
| 5 | Before using code patterns, verify installed versions match. If versions differ: |
| 6 | - Python: `pip show <package>` then `help(module.function)` to check signatures |
| 7 | - R: `packageVersion('<pkg>')` then `?function_name` to verify parameters |
| 8 | - CLI: `<tool> --version` then `<tool> --help` to confirm flags |
| 9 | |
| 10 | If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt. |
| 11 | |
| 12 | # ATAC-seq Quality Control |
| 13 | |
| 14 | **"Does my ATAC library pass ENCODE quality criteria?"** -> Compute the seven canonical metrics (depth, alignment rate, mitochondrial fraction, library complexity, fragment-size periodicity, TSS enrichment, FRiP) and compare against ENCODE 4 thresholds, then diagnose failures. |
| 15 | |
| 16 | - CLI: `picard CollectInsertSizeMetrics`, `samtools flagstat`, `samtools idxstats` |
| 17 | - CLI: `deeptools plotFingerprint`, `computeMatrix reference-point` + `plotProfile` |
| 18 | - R: `ATACseqQC::TSSEscore`, `ATACseqQC::fragSizeDist`, `ATACseqQC::PTscore` |
| 19 | - Python: custom NRF/PBC from coordinate hash; pyBigWig for TSS enrichment |
| 20 | |
| 21 | ## ENCODE 4 ATAC-seq Acceptance Thresholds |
| 22 | |
| 23 | | Metric | Definition | Ideal | Acceptable | Reject | Source | |
| 24 | |--------|-----------|-------|------------|--------|--------| |
| 25 | | Nuclear reads (after dedup, no chrM) | Mapped, MAPQ >= 30, non-chrM, deduped | >= 50M | 25-50M | < 25M | ENCODE 4 ATAC-seq Standards | |
| 26 | | Alignment rate | Mapped / total reads | >= 95% | 80-95% | < 80% | ENCODE 4 | |
| 27 | | Mitochondrial fraction | chrM / total mapped | < 5% (Omni-ATAC), < 20% (standard) | 20-50% | > 50% | Corces 2017 (Omni-ATAC) | |
| 28 | | NRF (Non-Redundant Fraction) | Distinct positions / total reads | >= 0.9 | 0.7-0.9 | < 0.7 | Landt 2012 | |
| 29 | | PBC1 (PCR Bottlenecking Coefficient 1) | Positions w/ 1 read / Positions w/ >= 1 read | >= 0.9 | 0.7-0.9 | < 0.7 | Landt 2012 | |
| 30 | | PBC2 | Positions w/ 1 read / Positions w/ 2 reads | >= 3.0 | 1.0-3.0 | < 1.0 | Landt 2012 | |
| 31 | | TSS enrichment (hg38, GENCODE v29) | Avg signal at TSS / avg flanking | >= 7 | 5-7 | < 5 | ENCODE 4 | |
| 32 | | FRiP (Fraction Reads in Peaks) | Reads in MACS peaks / total | >= 0.3 | 0.2-0.3 | < 0.2 | ENCODE 4, Landt 2012 | |
| 33 | | Insert-size periodicity | NFR + mono-nuc + di-nuc peaks visible | Clear 3+ peaks | NFR + mono only | Flat / single peak | Buenrostro 2013 | |
| 34 | |
| 35 | ENCODE thresholds are organism-specific. Mouse (mm10, GENCODE M21) TSS enrichment >= 5 is acceptable; non-model organisms have no published threshold (use cohort percentile rank instead). Methodology evolves; verify against the current ENCODE ATAC-seq Standards before reporting. |
| 36 | |
| 37 | ## TSS Enrichment: ENCODE Method vs ATACseqQC Method |
| 38 | |
| 39 | The two most common implementations DO NOT produce identical scores. |
| 40 | |
| 41 | | Method | Numerator | Denominator | Scaling | |
| 42 | |--------|-----------|-------------|---------| |
| 43 | | ENCODE pyTSSe / Kundaje gtsse | Mean signal in 100 bp window centered at TSS | Mean signal in 100 bp window at +/- 1900 to +/- 2000 bp (flanks) | Per-base normalization to flanks; reported as fold-enrichment | |
| 44 | | ATACseqQC TSSEscore | Sum signal in TSS +/- 100 bp | Sum signal at +/- 1000 bp flanking windows | Different window sizes; ratios are larger | |
| 45 | | deeptools plotProfile | Visual; numeric ratio not standardized | Reference-point matrix | No standard score; for visualization only | |
| 46 | |
| 47 | **Trigger:** Comparing a TSS score across studies. |
| 48 | |
| 49 | **Mechanism:** Different normalization windows shift the absolute number; ATACseqQC's TSSEscore is typically 2-3x ENCODE's because of the wider flank. |
| 50 | |
| 51 | **Symptom:** Reported score 21 vs ENCODE-ideal 7 mismatch. Likely the calculator was ATACseqQC; the equivalent ENCODE score might be 8. |
| 52 | |
| 53 | **Fix:** State which implementation was used. For ENCODE comparisons, use `pyTSSe` (Kundaje lab) or implement the ENCODE recipe directly. |
| 54 | |
| 55 | ```python |
| 56 | import numpy as np |
| 57 | import pyBigWig |
| 58 | |
| 59 | def encode_tss_enrichment(bw_path, tss_bed, flank=2000): |
| 60 | """ENCODE-style TSS enrichment: signal at TSS center / signal at flanks.""" |
| 61 | bw = pyBigWig.open(bw_path) |
| 62 | profiles = [] |
| 63 | for line in open(tss_bed): |
| 64 | chrom, start, end, *rest = line.strip().split('\t') |
| 65 | tss = int(start) |
| 66 | strand = rest[2] if len(rest) > 2 else '+' |
| 67 | try: |
| 68 | vals = bw.values(chrom, tss - flank, tss + flank) |
| 69 | if vals is None or len(vals) != 2 * flank: continue |
| 70 | if strand == '-': vals = vals[::-1] |
| 71 | pr |