Merge pull request #3 from ewels/detect-duplicates

Add duplicate-marking validation for BAM files

Author: ewels

AGENTS.md

@@ -164,6 +164,22 @@ All three must pass. Uses `dtolnay/rust-toolchain@stable` and `Swatinem/rust-cac
 | `env_logger`   | Log output backend               |
 | `indexmap`     | Insertion-order-preserving maps  |
 
+## Duplicate Marking Validation
+
+RustQC verifies that BAM files have been processed by a duplicate-marking tool before
+running the RNA duplication analysis. This is implemented in two layers:
+
+1. **Pre-flight `@PG` header check** (`counting::verify_duplicates_marked`): Parses the
+   BAM header text for `@PG` lines and checks the `ID:` and `PN:` fields against
+   `KNOWN_DUP_MARKERS` (Picard MarkDuplicates, samblaster, sambamba, biobambam, etc.).
+   Exits with `anyhow::bail!()` if no known tool is found.
+2. **Post-hoc zero-duplicates check**: After counting, if `total_dup == 0` with
+   `total_mapped > 0`, bails with an error — catches cases where the header is present
+   but duplicates weren't actually flagged.
+
+Both checks are skipped when `--skip-dup-check` is passed (stored as `RnaArgs.skip_dup_check`,
+forwarded to `count_reads()` as the `skip_dup_check: bool` parameter).
+
 ## Notes for Agents
 
 - The codebase is a pure binary crate with no library target.

README.md

@@ -136,6 +136,15 @@ rustqc rna <BAM> <GTF> [OPTIONS]
 | `--threads <N>` | `1` | Number of threads for parallel BAM processing |
 | `--outdir <DIR>` | `.` | Output directory |
 | `--config <FILE>` | none | Path to a YAML configuration file (see [Configuration](#configuration)) |
+| `--skip-dup-check` | `false` | Skip verification that duplicates have been marked in the BAM file (see [Duplicate marking](#duplicate-marking)) |
+
+### Duplicate marking
+
+RustQC requires that the input BAM file has been processed by a duplicate-marking tool such as [Picard MarkDuplicates](https://broadinstitute.github.io/picard/command-line-overview.html#MarkDuplicates), [samblaster](https://github.com/GregoryFaust/samblaster), or [sambamba markdup](https://lomereiter.github.io/sambamba/). These tools set the SAM flag `0x400` on PCR/optical duplicate reads, which RustQC uses to compute duplication rates.
+
+Before processing, RustQC checks the BAM `@PG` header lines for known duplicate-marking programs. If none are found, it exits with an error explaining how to mark duplicates. As a secondary safeguard, if processing completes but zero duplicate-flagged reads are found among mapped reads, RustQC also exits with an error.
+
+If you are confident that your BAM file has duplicates correctly flagged despite the header check failing, you can bypass the verification with `--skip-dup-check`.
 
 ### Example
 

src/cli.rs

@@ -60,6 +60,14 @@ pub struct RnaArgs {
     /// Path to a YAML configuration file (e.g. chromosome name mapping)
     #[arg(short, long, value_name = "CONFIG")]
     pub config: Option<String>,
+
+    /// Skip the check for duplicate-marking tools in the BAM header.
+    ///
+    /// By default, RustQC verifies that the BAM file has been processed by a
+    /// duplicate-marking tool (e.g. Picard MarkDuplicates, samblaster) before
+    /// running. Use this flag to bypass that check.
+    #[arg(long, default_value_t = false)]
+    pub skip_dup_check: bool,
 }
 
 /// Parse command-line arguments and return the Cli struct.

src/counting.rs

@@ -67,6 +67,111 @@ impl GeneIdInterner {
     }
 }
 
+// ===================================================================
+// BAM header validation
+// ===================================================================
+
+/// Known duplicate-marking tool identifiers.
+///
+/// These strings are matched (case-insensitively) against the `ID:` and `PN:`
+/// (program name) fields of `@PG` header entries. Matching is restricted to these
+/// fields to avoid false positives from command-line strings or unrelated tools
+/// (e.g., `picard SortSam` or `sambamba sort`).
+const KNOWN_DUP_MARKERS: &[&str] = &[
+    "markduplicates",
+    "samblaster",
+    "sambamba markdup",
+    "sambamba_markdup",
+    "biobambam",
+    "estreamer",
+    "fgbio",
+    "umis",
+    "umi_tools",
+    "umi-tools",
+    "gencore",
+    "gatk markduplicates",
+    "sentieon dedup",
+];
+
+/// Extracts the value of a specific tag from a SAM header line.
+///
+/// For example, given the line `@PG\tID:MarkDuplicates\tPN:MarkDuplicates` and
+/// the tag `"ID"`, returns `Some("MarkDuplicates")`.
+fn extract_header_tag<'a>(line: &'a str, tag: &str) -> Option<&'a str> {
+    let prefix = format!("{}:", tag);
+    line.split('\t')
+        .find(|field| field.starts_with(&prefix))
+        .map(|field| &field[prefix.len()..])
+}
+
+/// Checks whether a SAM header text contains evidence of a duplicate-marking tool.
+///
+/// Only the `ID:` and `PN:` fields of `@PG` lines are inspected to avoid false
+/// positives from command-line arguments or unrelated tools that happen to share
+/// a name (e.g., `picard SortSam`).
+///
+/// Returns `true` if a known duplicate-marking tool is found.
+fn header_has_dup_marker(header_text: &str) -> bool {
+    for line in header_text.lines() {
+        // Only inspect @PG (program) header lines
+        if !line.starts_with("@PG") {
+            continue;
+        }
+
+        // Extract only the ID and PN fields for matching
+        let id = extract_header_tag(line, "ID").unwrap_or("");
+        let pn = extract_header_tag(line, "PN").unwrap_or("");
+        let id_lower = id.to_lowercase();
+        let pn_lower = pn.to_lowercase();
+
+        if KNOWN_DUP_MARKERS
+            .iter()
+            .any(|marker| id_lower.contains(marker) || pn_lower.contains(marker))
+        {
+            debug!("Found duplicate-marking tool in @PG header: {}", line);
+            return true;
+        }
+    }
+    false
+}
+
+/// Checks the BAM `@PG` header lines for evidence that a duplicate-marking tool has been run.
+///
+/// Returns `Ok(())` if a known duplicate-marking tool is found, or an error with
+/// a descriptive message if none is detected.
+///
+/// # Arguments
+///
+/// * `header` - The BAM header view to inspect
+/// * `bam_path` - The BAM file path (used in the error message)
+fn verify_duplicates_marked(header: &bam::HeaderView, bam_path: &str) -> Result<()> {
+    let header_text = String::from_utf8_lossy(header.as_bytes());
+    if header_has_dup_marker(&header_text) {
+        return Ok(());
+    }
+
+    anyhow::bail!(
+        "No duplicate-marking tool found in BAM header of '{}'.\n\
+         \n\
+         RustQC requires that BAM files have duplicates marked (SAM flag 0x400)\n\
+         but NOT removed. The BAM @PG header lines do not contain evidence of a\n\
+         known duplicate-marking tool.\n\
+         \n\
+         Please run one of the following tools before using RustQC:\n\
+         \n\
+           - Picard MarkDuplicates: picard MarkDuplicates I=input.bam O=marked.bam M=metrics.txt\n\
+           - samblaster:            samtools view -h input.bam | samblaster | samtools view -bS - > marked.bam\n\
+           - sambamba markdup:      sambamba markdup input.bam marked.bam\n\
+         \n\
+         If you are certain that duplicates are already marked, use --skip-dup-check to bypass this check.",
+        bam_path
+    )
+}
+
+// ===================================================================
+// BAM flag constants
+// ===================================================================
+
 /// Flag indicating the read is a PCR or optical duplicate (0x400).
 const BAM_FDUP: u16 = 0x400;
 /// Flag indicating the read is unmapped (0x4).
@@ -681,6 +786,8 @@ fn process_chromosome_batch(
 /// * `stranded` - Library strandedness (0, 1, or 2)
 /// * `paired` - Whether the library is paired-end
 /// * `threads` - Number of threads for BAM processing
+/// * `skip_dup_check` - If true, skip the BAM header check for duplicate-marking tools
+#[allow(clippy::too_many_arguments)]
 pub fn count_reads(
     bam_path: &str,
     genes: &IndexMap<String, Gene>,
@@ -689,18 +796,28 @@ pub fn count_reads(
     threads: usize,
     chrom_mapping: &HashMap<String, String>,
     chrom_prefix: Option<&str>,
+    skip_dup_check: bool,
 ) -> Result<CountResult> {
     // Build gene ID interner for allocation-free lookups in the hot loop
     let interner = GeneIdInterner::from_genes(genes);
 
     // Build spatial index (uses interned gene indices)
     let index = build_index(genes, &interner);
 
-    // Get chromosome names from header using a temporary reader
+    // Get chromosome names from header using a temporary reader,
+    // and verify that duplicates have been marked in the BAM file.
     let tid_to_name: Vec<String> = {
         let bam = bam::Reader::from_path(bam_path)
             .with_context(|| format!("Failed to open BAM file: {}", bam_path))?;
         let header = bam.header().clone();
+
+        // Check for evidence of duplicate-marking in @PG header lines
+        if skip_dup_check {
+            info!("Skipping duplicate-marking verification (--skip-dup-check)");
+        } else {
+            verify_duplicates_marked(&header, bam_path)?;
+        }
+
         (0..header.target_count())
             .map(|tid| String::from_utf8_lossy(header.tid2name(tid)).to_string())
             .collect()
@@ -1108,6 +1225,27 @@ pub fn count_reads(
         merged.total_fragments
     );
 
+    // Verify that at least some reads were flagged as duplicates.
+    // Even if the @PG header check passed (or was skipped), it's possible that
+    // duplicates were not actually marked. Warn the user in that case.
+    if merged.total_dup == 0 && merged.total_mapped > 0 && !skip_dup_check {
+        anyhow::bail!(
+            "No duplicate-flagged reads found among {} mapped reads in '{}'.\n\
+             \n\
+             Although the BAM header suggests a duplicate-marking tool was run,\n\
+             no reads have the duplicate flag (0x400) set. This likely means\n\
+             duplicates were removed rather than marked, or the tool did not\n\
+             flag any duplicates.\n\
+             \n\
+             RustQC requires duplicates to be marked (flagged) but NOT removed.\n\
+             Please re-run your duplicate-marking tool without removing duplicates.\n\
+             \n\
+             If this is expected, use --skip-dup-check to bypass this check.",
+            merged.total_mapped,
+            bam_path
+        );
+    }
+
     // Detect chromosome name mismatch
     let genes_with_reads = merged
         .gene_counts
@@ -1206,4 +1344,117 @@ mod tests {
         // Read2 on - strand: effective + -> matches + gene
         assert!(strand_matches(true, false, true, '+', 1));
     }
+
+    // ---------------------------------------------------------------
+    // Duplicate marking validation tests
+    // ---------------------------------------------------------------
+
+    #[test]
+    fn test_extract_header_tag() {
+        let line = "@PG\tID:MarkDuplicates\tPN:MarkDuplicates\tVN:2.27.4\tCL:picard MarkDuplicates I=in.bam O=out.bam";
+        assert_eq!(extract_header_tag(line, "ID"), Some("MarkDuplicates"));
+        assert_eq!(extract_header_tag(line, "PN"), Some("MarkDuplicates"));
+        assert_eq!(extract_header_tag(line, "VN"), Some("2.27.4"));
+        assert_eq!(
+            extract_header_tag(line, "CL"),
+            Some("picard MarkDuplicates I=in.bam O=out.bam")
+        );
+        assert_eq!(extract_header_tag(line, "XX"), None);
+    }
+
+    #[test]
+    fn test_dup_check_picard_markduplicates() {
+        let header = "@HD\tVN:1.6\tSO:coordinate\n\
+                       @PG\tID:MarkDuplicates\tPN:MarkDuplicates\tVN:2.27.4";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_samblaster() {
+        let header = "@HD\tVN:1.6\n\
+                       @PG\tID:samblaster\tPN:samblaster\tVN:0.1.26";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_sambamba_markdup() {
+        let header = "@HD\tVN:1.6\n\
+                       @PG\tID:sambamba_markdup\tPN:sambamba markdup";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_biobambam() {
+        let header = "@HD\tVN:1.6\n\
+                       @PG\tID:biobambam2\tPN:biobambam";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_case_insensitive() {
+        // MarkDuplicates with different casing
+        let header = "@PG\tID:MARKDUPLICATES\tPN:markduplicates";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_no_dup_marker() {
+        // Header with only alignment tools, no dup marker
+        let header = "@HD\tVN:1.6\tSO:coordinate\n\
+                       @PG\tID:bwa\tPN:bwa\tVN:0.7.17\n\
+                       @PG\tID:samtools\tPN:samtools\tVN:1.17";
+        assert!(!header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_empty_header() {
+        assert!(!header_has_dup_marker(""));
+        assert!(!header_has_dup_marker("@HD\tVN:1.6\tSO:coordinate"));
+    }
+
+    #[test]
+    fn test_dup_check_picard_sortsam_no_false_positive() {
+        // Picard SortSam should NOT match — only the CL field mentions "picard"
+        let header = "@PG\tID:SortSam\tPN:SortSam\tCL:picard SortSam I=in.bam O=out.bam";
+        assert!(!header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_picard_collect_metrics_no_false_positive() {
+        // Another Picard tool that should not match
+        let header =
+            "@PG\tID:CollectInsertSizeMetrics\tPN:CollectInsertSizeMetrics\tCL:picard CollectInsertSizeMetrics";
+        assert!(!header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_sambamba_sort_no_false_positive() {
+        // sambamba sort should NOT match
+        let header = "@PG\tID:sambamba\tPN:sambamba sort";
+        assert!(!header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_multiple_pg_lines() {
+        // MarkDuplicates appears as a later @PG entry in the chain
+        let header = "@HD\tVN:1.6\tSO:coordinate\n\
+                       @PG\tID:bwa\tPN:bwa\tVN:0.7.17\n\
+                       @PG\tID:samtools\tPN:samtools\tVN:1.17\n\
+                       @PG\tID:MarkDuplicates\tPN:MarkDuplicates\tPP:samtools";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_gatk_markduplicates() {
+        let header = "@PG\tID:GATK MarkDuplicates\tPN:GATK MarkDuplicates\tVN:4.4.0";
+        assert!(header_has_dup_marker(header));
+    }
+
+    #[test]
+    fn test_dup_check_dup_marker_only_in_cl_no_match() {
+        // A tool that mentions MarkDuplicates only in the CL field should not match
+        // (e.g., a wrapper script that calls MarkDuplicates but has its own ID/PN)
+        let header = "@PG\tID:my_pipeline\tPN:my_pipeline\tCL:java -jar picard.jar MarkDuplicates";
+        assert!(!header_has_dup_marker(header));
+    }
 }

src/main.rs

@@ -86,6 +86,7 @@ fn run_rna(args: cli::RnaArgs) -> Result<()> {
         args.threads,
         &chrom_mapping,
         config.chromosome_prefix(),
+        args.skip_dup_check,
     )?;
     info!(
         "Counting complete in {:.2}s",

tests/integration_test.rs

@@ -39,6 +39,7 @@ fn run_rustqc(outdir: &str) -> std::process::Output {
             "tests/data/test.gtf",
             "--outdir",
             outdir,
+            "--skip-dup-check",
         ])
         .output()
         .expect("Failed to execute rustqc")