metagenomics/contamination-controls/SKILL.md
Cleans a shotgun metagenome of everything that is not the target community before profiling - host-read depletion (Hostile, bowtie2/T2T-CHM13), reagent/kitome contamination control with blanks and decontam, mock-community validation, and depth-adequacy checks (Nonpareil). Covers why a metagenomic result is a position in a choice-chain rather than a direct observation, why extraction is the experiment, why a low-biomass community can be entirely kitome, why absence means not-detectable-by-this-chain, and why a confident classifier call can still be wrong when the reference is contaminated. Use when designing controls, removing host reads, identifying reagent contaminants, validating with mocks, or judging whether a low-biomass result is real. For adapter/quality trimming see read-qc; for MAG-level decontamination see genome-assembly/metagenome-assembly.
npx skillsauth add GPTomics/bioSkills bio-metagenomics-contamination-controlsInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
Reference examples tested with: decontam 1.22+, Hostile 1.1+, Bowtie2 2.5+, Nonpareil 3.4+, pandas 2.2+, R 4.3+.
Before using code patterns, verify installed versions match. If versions differ:
packageVersion('decontam') then ?isContaminant to verify parametershostile --version, nonpareil -h to confirm flags and indexespip show <package> then help(module.function) to check signaturesIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
The controls define the result: an extraction blank defines the kitome for its batch/lot, a mock community defines the limit of detection and the extraction-lysis bias, and the host reference (prefer T2T-CHM13 over GRCh38) defines what host is removed. Record the extraction kit and lot, the host index, the blanks, the mock version (whole-cell vs DNA), and the reads removed at each step.
"Is this signal real, or did my pipeline create it?" -> Remove host reads, define the kitome with blanks, validate with a mock, and confirm depth - because a low-biomass community can be entirely reagent contamination.
decontam::isContaminant(seqtab, conc=, neg=, method='combined') on the classifier output tablehostile clean --fastq1 R1.fq.gz --fastq2 R2.fq.gz --index human-t2t-hlaScope: sample-level pre-analysis cleanup and controls - host depletion, kitome/blank/mock controls, decontam, depth adequacy. Adapter/quality trimming mechanics -> read-qc/adapter-trimming, read-qc/quality-filtering. MAG-level decontamination (CheckM2/GUNC chimerism, FCS-GX foreign sequence) -> genome-assembly/metagenome-assembly - a different, genome-level problem. Classification -> kraken-classification, metaphlan-profiling.
A metagenomic profile is the product of a chain of choices - extraction, host/contaminant depletion, depth, read-vs-assembly, classifier, database, normalization - and each link silently sets what is observable. The community is never observed directly; the report is the community as refracted by this pipeline. Three consequences a newcomer misses:
Lysis efficiency is taxon-dependent: tough-walled Gram-positives (Firmicutes, Staphylococcus, Enterococcus), endospores (Bacillus, Clostridium), acid-fast Mycobacterium, and fungi/archaea resist lysis and are under-represented unless bead-beating is used. Gentle/enzymatic kits inflate easy-to-lyse Gram-negatives - so a Firmicutes:Bacteroidetes shift can be an extraction artifact. Extraction had the largest effect on observed composition across 21 protocols (Costea 2017 Nat Biotechnol 35:1069). Use bead-beating, hold one method constant across a study, validate lysis with a whole-cell mock, and report kit and lot. Note the tradeoff: aggressive bead-beating shears DNA and hurts long-read assembly, so the best extraction depends on the read-vs-assembly choice.
| Scenario | Recommended | Why | |----------|-------------|-----| | Host-associated sample (gut, oral, skin, tissue) | host-deplete first (Hostile, T2T-CHM13) | host reads waste depth and leak false calls; also a data-sharing/ethics requirement | | Low-biomass sample (skin, BAL, CSF, tissue, blood) | blanks + DNA quantification + decontam mandatory | the lower the biomass, the larger the kitome fraction | | Novel taxa claimed in low biomass | treat as kitome until proven (canonical genera) | placenta/tumor "microbiomes" were largely kitome | | Need limit of detection / lysis check | run a mock (ZymoBIOMICS whole-cell) | the only sample with a known answer | | Is my depth enough for this question? | Nonpareil coverage curve | depth sets the detection limit; host depletion halves usable depth | | Confident classifier call, odd taxon | suspect a contaminated reference | confidence is not correctness if the reference is mislabeled | | Adapter/quality trimming | -> read-qc | this skill owns metagenomics-specific cleanup, not generic trimming | | MAG chimerism / foreign sequence in a bin | -> genome-assembly/metagenome-assembly | genome-level decontamination is a different problem |
# Hostile removes >99.5% of human reads while discarding far fewer microbial reads than naive mapping.
# Prefer the T2T-CHM13-based index over GRCh38; high-sensitivity Bowtie2 drives removal more than the reference.
hostile clean --fastq1 sample_R1.fq.gz --fastq2 sample_R2.fq.gz \
--index human-t2t-hla --aligner bowtie2
# Report the reads removed - it is a QC metric, not a footnote. For long reads use --aligner minimap2.
Remove host for two reasons: analytical (depth, false positives, runtime) and ethical (raw human-associated reads carry identifiable host genotype; depleting before deposit is increasingly required). Wet-lab depletion (saponin/DNase, methyl-CpG capture) saves sequencing but adds its own bias - a genuine tradeoff.
Goal: Separate real low-abundance taxa from the kitome using blanks and DNA concentration.
Approach: Run decontam on the classifier output table (taxa x samples) using the frequency signal (contaminants scale inversely with input DNA) and the prevalence signal (contaminants are enriched in blanks); raise the prevalence threshold for low biomass.
library(decontam)
# seqtab: samples x taxa from the Bracken/MetaPhlAn table; conc: per-sample DNA concentration; neg: TRUE for blanks.
contam <- isContaminant(seqtab, conc = dna_conc, neg = is_blank, method = 'combined', threshold = 0.1)
# Low-biomass studies: use the prevalence method at the more aggressive threshold 0.5, and inspect the calls.
contam_lowbio <- isContaminant(seqtab, neg = is_blank, method = 'prevalence', threshold = 0.5, batch = batch_id)
seqtab_clean <- seqtab[, !contam$contaminant]
decontam runs per batch (batch=) because the kitome differs by lot/run. Always inspect the called contaminants against the canonical kitome genera (Bradyrhizobium, Ralstonia, Burkholderia, Pseudomonas, Acinetobacter, Sphingomonas, Methylobacterium, Stenotrophomonas) rather than applying blindly - over-aggressive removal deletes real taxa.
# Nonpareil estimates how much of the community's sequence space you have sampled, without assembly or a DB.
nonpareil -s reads.fasta -T kmer -f fasta -b sample_np
# Plot the coverage-vs-effort curve in R (Nonpareil.curve); a non-detection below the implied limit is meaningless.
Depth is set by the question: dominant taxa need a few million reads, rare-pathogen detection sets a limit of detection, and strain SNVs need high per-genome coverage. Host depletion can silently halve usable depth - budget for it.
Trigger: a Firmicutes:Bacteroidetes shift or "low Gram-positive" community from a gentle-lysis kit. Mechanism: taxon-dependent lysis under-represents tough-walled organisms. Symptom: composition differences tracking the kit, not the sample. Fix: bead-beating, one method held constant, a whole-cell mock to prove hard taxa are lysed.
Trigger: reporting novel low-abundance taxa from skin/BAL/tissue/blood without blanks. Mechanism: reagent DNA is a fixed dose; at low biomass it dominates the signal. Symptom: canonical kitome genera presented as discovery. Fix: extraction blanks through the full workflow, DNA quantification, decontam (prevalence + frequency), skepticism toward the kitome genera.
Trigger: "taxon/function not present" or "low diversity." Mechanism: detection is bounded by depth, database, extraction, and depletion. Symptom: a negative interpreted as biology. Fix: report the classified fraction and the limit of detection; state which link is responsible before interpreting absence.
Trigger: trusting a high-confidence classifier assignment. Mechanism: >2 million GenBank/RefSeq entries carry mislabeled or chimeric sequence (Steinegger & Salzberg 2020 Genome Biol 21:115). Symptom: a confident, systematic wrong assignment (the classic stray human/vector in a microbial genome). Fix: treat confidence as not equal to correctness; cross-check surprising calls against a cleaner database.
| Threshold | Source | Rationale | |-----------|--------|-----------| | decontam threshold 0.1 default; 0.5 prevalence for low biomass | Davis 2018 Microbiome 6:226 | aggressive prevalence call needed when the kitome dominates | | >= 1 extraction blank per batch | Salter 2014 BMC Biol 12:87 | blanks define the kitome for that lot; more for very low biomass | | Hostile removes > 99.5% host | Constantinides 2023 Bioinformatics 39:btad728 | high host removal with low microbial loss | | Prefer T2T-CHM13 over GRCh38 + mask rDNA | host-removal practice | GRCh38 gaps let host reads escape; rDNA masking spares microbial reads | | Whole-cell vs DNA mock | mock-standard practice | whole-cell tests extraction/lysis; DNA tests classifier/library only | | Nonpareil coverage before interpreting absence | Rodriguez-R 2018 mSystems 3:e00039-18 | a non-detection below the limit of detection is uninformative |
| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| decontam finds nothing useful | no blanks or DNA concentration supplied | add blanks (neg=) and/or DNA quant (conc=); run per batch |
| Host removal leaves human reads | GRCh38 with gaps, low-sensitivity aligner | use a T2T-CHM13 index and high-sensitivity Bowtie2 |
| Real microbes deleted in host removal | rDNA / conserved regions not masked | mask host rDNA; check microbial reads removed |
| Low-biomass "novel taxon" not reproducible | kitome | blanks + decontam; check canonical kitome genera |
| Cross-study profiles disagree | different extraction/depth/DB chains | hold the chain constant; do not meta-analyze across links |
tools
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
development
Detect, date, and contextualize whole-genome duplication (WGD / paleopolyploidy) events using wgd v2 (Chen et al 2024), KsRates (Sensalari 2022 substitution-rate-corrected Ks dating), DupGen_finder (Qiao 2019), MAPS (Li 2018 phylogenomic), POInT (Conant 2008 ordered-block), SLEDGe (2024 ML-based), Whale.jl (Bayesian DL+WGD), and synteny-anchored paranome construction. Use when identifying ancient polyploidy from Ks distributions and synteny block analysis, positioning WGD events relative to speciation, distinguishing tandem from segmental from WGD duplications, dating the 2R/3R vertebrate / fish / salmonid WGDs, building paranome and Ks-age mixture models, applying KsRates substitution-rate correction across lineages, or testing alternative biased-fractionation / dosage-balance models post-WGD.
tools
Build whole-genome alignments using Progressive Cactus (Armstrong 2020 reference-free clade-level WGA), Minigraph-Cactus (Hickey 2024 pangenome-aware), LASTZ chain/net (UCSC pipeline), MUMmer4 (Marçais 2018 pairwise), minimap2 -x asm5/10/20 (Li 2018 fast pairwise), AnchorWave (Song 2022 WGD-aware), and Mauve / progressiveMauve (bacterial). Operates the HAL toolkit (Hickey 2013) for downstream extraction including halSynteny, halLiftover, halBranchMutations, and hal2maf. Use when constructing multi-species alignments for comparative-annotation projection (TOGA), synteny detection, conservation analyses (phyloP / PhastCons), or pangenome graph construction; selecting between reference-free (Cactus) and reference-anchored (LASTZ chains/nets) approaches; tuning sensitivity for closely vs distantly related genomes; or producing HAL files for genome-wide downstream tools.
development
Detect syntenic blocks and structural rearrangements between genomes using MCScanX (Wang 2012), JCVI/MCScan (Tang 2008 Python), GENESPACE (Lovell 2022) for orthology-anchored riparian visualization, SyRI for structural variation, AnchorWave for sequence-level synteny, i-ADHoRe 3.0 for highly diverged species, SynNet for synteny networks, and ntSynt for multi-genome macrosynteny. Use when identifying collinear gene blocks across species, distinguishing macrosynteny from microsynteny, detecting inversions/translocations/duplications, anchoring orthology in WGD lineages, producing publication riparian plots, computing synteny block age via Ks (cross-references whole-genome-duplication), or running synteny-aware ortholog inference in polyploids.