$npx -y skills add GPTomics/bioSkills --skill duplicate-handlingMark and remove PCR/optical duplicates using samtools fixmate and markdup. Use when preparing alignments for variant calling or when duplicate reads would bias analysis.
| 1 | ## Version Compatibility |
| 2 | |
| 3 | Reference examples tested with: picard 3.1+, pysam 0.22+, samtools 1.19+ |
| 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 | - CLI: `<tool> --version` then `<tool> --help` to confirm flags |
| 8 | |
| 9 | If code throws ImportError, AttributeError, or TypeError, introspect the installed |
| 10 | package and adapt the example to match the actual API rather than retrying. |
| 11 | |
| 12 | # Duplicate Handling |
| 13 | |
| 14 | **"Remove PCR duplicates from my BAM file"** -> Mark or remove duplicate reads using the fixmate-sort-markdup pipeline to prevent duplicate bias in variant calling. |
| 15 | - CLI: `samtools fixmate`, `samtools markdup` (samtools) |
| 16 | - Python: `pysam.fixmate()`, `pysam.markdup()` (pysam) |
| 17 | |
| 18 | Mark and remove PCR/optical duplicates using samtools. |
| 19 | |
| 20 | ## Why Remove Duplicates? |
| 21 | |
| 22 | PCR duplicates are identical copies of the same original molecule, created during library preparation. They inflate coverage, bias allele frequencies, and create false positive variant calls. Optical duplicates are flowcell-proximity artifacts: on unpatterned flowcells they arise when the imaging software splits one cluster into two adjacent calls; on patterned flowcells (NovaSeq, NovaSeq X, NextSeq 1000/2000, HiSeq X/4000) the dominant source is ExAmp (exclusion-amplification) "pad-hopping", where a library molecule re-seeds a nearby nanowell. |
| 23 | |
| 24 | ## When to Mark Duplicates -- and When NOT To |
| 25 | |
| 26 | Standard `samtools markdup` is the right tool for some assays and actively harmful for others. The decision is assay-driven: |
| 27 | |
| 28 | | Assay | Standard markdup? | Recommended approach | |
| 29 | |-------|------------------|----------------------| |
| 30 | | Germline WGS / WES (PCR or PCR-free) | YES | `samtools markdup` (PCR-free still has ~0.5% optical duplicates on patterned flowcells) | |
| 31 | | Somatic tumor/normal (no UMI) | YES | Same | |
| 32 | | Exome / target capture | YES (20-50% expected) | `samtools markdup` | |
| 33 | | ChIP-seq | **MARK, do not remove** | Then use peak caller's auto-dup logic (`macs3 --keep-dup auto`) | |
| 34 | | CUT&RUN / CUT&Tag | MARK, do not remove | Same | |
| 35 | | ATAC-seq | YES, BEFORE Tn5 +4/-5 shift | Then shift coords for footprinting | |
| 36 | | Bulk RNA-seq (no UMIs) | **NO** | Duplicates are biological at highly-expressed loci; removing them biases DE proportional to expression | |
| 37 | | Bulk RNA-seq (with UMIs) | NO | umi_tools dedup | |
| 38 | | scRNA (10x, STARsolo, drop-seq) | **NO** | umi_tools dedup with CB+UB tags, or rely on Cell Ranger UMI counts | |
| 39 | | ctDNA / liquid biopsy / deep panel (UMI) | NO | fgbio GroupReadsByUmi -> CallDuplexConsensusReads | |
| 40 | | Twist / IDT / Roche UMI capture | NO | fgbio or Picard UmiAwareMarkDuplicatesWithMateCigar | |
| 41 | | Amplicon / hotspot panel (no UMI) | **NO** | Every read is a "duplicate" by coordinate; markdup erases the dataset. Use `samtools ampliconclip` instead -- see alignment-amplicon-clipping. | |
| 42 | | Amplicon / hotspot panel (UMI) | NO | fgbio consensus | |
| 43 | | Long-read native (ONT, PacBio HiFi unamplified) | NO | No PCR step; markdup is meaningless | |
| 44 | | PacBio HiFi amplicon | YES (rare) | `pbmarkdup` | |
| 45 | | Ancient DNA (aDNA) | YES + mapDamage | Run markdup, then mapDamage `--rescale` before variant calling | |
| 46 | | Microbiome 16S/ITS | NO | Read counts encode community structure | |
| 47 | |
| 48 | If the BAM came from 10x Cell Ranger / STARsolo and `samtools markdup` produces a 50-95% duplicate rate, it is the wrong tool, not a bug. |
| 49 | |
| 50 | ## Tool Selection: markdup vs Picard vs UMI-aware |
| 51 | |
| 52 | | Tool | Speed | Threading | Optical | UMI | Notes | |
| 53 | |------|-------|-----------|---------|-----|-------| |
| 54 | | `samtools markdup` | Fast | Yes | Yes (`-d`) | Limited (`--barcode-tag` exact-match) | Fast production choice (nf-core/sarek defaults to GATK MarkDuplicates) | |
| 55 | | `picard MarkDuplicates` | Slow | No | Yes | UmiAware variant (BETA, transcriptome bug) | GATK Best Practices reference | |
| 56 | | `biobambam2 bammarkduplicates2` | Fastest | Yes | Yes | No | Sanger / 1KGP pipelines | |
| 57 | | `samblaster` | Streaming, fast | No | Optional | No | Pipe directly from aligner; no name sort | |
| 58 | | `sambamba markdup` | Fast | Yes | Yes | No | Less actively maintained | |
| 59 | | `fgbio GroupReadsByUmi` + `CallMolecularConsensusReads` | Fast | Yes | n/a | **Best UMI tool** | Graph-based; supports duplex | |
| 60 | | `umi_tools dedup` | Slow | No | n/a | Yes (mature) | Reference for scRNA / bulk UMI | |
| 61 | | `pbmarkdup` | Fast | Yes | n/a | n/a | PacBio HiFi amplicons only | |
| 62 | |
| 63 | Picard `UmiAwareMarkDuplicatesWithMateCigar` is BETA and has known bugs on transcriptome-aligned BAMs (silently keeps duplicates). Avoid for RNA-seq UMIs. |
| 64 | |
| 65 | ## Optical Distance Is Platform-Specific |
| 66 | |
| 67 | `samtools markdup` default is `-d 0`, meaning **optical-duplicate detection is disabled by default**. Set explicitly per platform: |
| 68 | |
| 69 | | Platform |