dupblaster 0.3.0

Fast duplicate marking for query-grouped SAM/BAM files, inspired by samblaster and Picard MarkDuplicates
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

//! Tests for dupblaster's handling of reads whose 5' clip extends past a
//! contig edge by more than `--max-read-length`. Such reads can't be
//! represented exactly on dupblaster's synthetic single-axis genome, so their
//! coordinate is clamped into the contig. Clamping is a correctness
//! compromise (potential false duplicates among clamped reads), so dupblaster
//! still writes its output, warns, and exits non-zero — and bumping
//! `--max-read-length` past the clip makes the run clean again.

mod helpers;

use helpers::*;

/// A read mapped near a contig start with thousands of bases of leading
/// hard-clip clips past the default (1000 bp) padding. dupblaster clamps it,
/// warns, and exits non-zero — but still writes the output.
#[test]
fn clip_beyond_padding_clamps_warns_and_exits_nonzero() {
    let env = TestEnv::new();
    let out_bam = env._tmp.path().join("out.bam");
    // A proper pair whose first end has a 5000-base leading hard clip at
    // chr1:1. forward 5' pos = 1 - 5000 ≪ 0, clipping ~4000 bp past the
    // 1000 bp padding.
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("r1", 99, "chr1", 1, "5000H50M", "=", 200, 249)
        .rec_simple("r1", 147, "chr1", 200, "50M", "=", 1, -249)
        .write_to(&env.input);

    let r = dupblaster()
        .args(["-i"])
        .arg(&env.input)
        .args(["-o"])
        .arg(&out_bam)
        .args(["--metrics-prefix"])
        .arg(metrics_prefix_for(&out_bam))
        .output()
        .unwrap();
    assert!(!r.status.success(), "clamping should make dupblaster exit non-zero");
    let stderr = String::from_utf8_lossy(&r.stderr);
    assert!(stderr.contains("clamped"), "stderr should mention clamping, got: {stderr}");
    assert!(
        stderr.contains("max-read-length"),
        "stderr should point at --max-read-length, got: {stderr}"
    );
    // Output is still written despite the non-zero exit: the pair round-trips.
    let records = read_records(&out_bam);
    assert_eq!(records.len(), 2, "output should still be written when clamping");
}

/// The same input runs clean (exit 0, no clamp warning) once
/// `--max-read-length` is bumped past the clip.
#[test]
fn bumping_max_read_length_avoids_clamping() {
    let env = TestEnv::new();
    let out_bam = env._tmp.path().join("out.bam");
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("r1", 99, "chr1", 1, "5000H50M", "=", 200, 249)
        .rec_simple("r1", 147, "chr1", 200, "50M", "=", 1, -249)
        .write_to(&env.input);

    let r = dupblaster()
        .args(["-i"])
        .arg(&env.input)
        .args(["-o"])
        .arg(&out_bam)
        .args(["--metrics-prefix"])
        .arg(metrics_prefix_for(&out_bam))
        .args(["--max-read-length", "6000"])
        .output()
        .unwrap();
    assert!(
        r.status.success(),
        "no clamping expected with a larger --max-read-length: {}",
        String::from_utf8_lossy(&r.stderr)
    );
    let stderr = String::from_utf8_lossy(&r.stderr);
    assert!(!stderr.contains("clamped"), "no clamp warning expected, got: {stderr}");
    let records = read_records(&out_bam);
    assert_eq!(records.len(), 2);
}

/// A normal pair well inside the contig neither clamps nor exits non-zero.
#[test]
fn normal_reads_do_not_clamp() {
    let env = TestEnv::new();
    let out_bam = env._tmp.path().join("out.bam");
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("r1", 99, "chr1", 5000, "50M", "=", 5200, 250)
        .rec_simple("r1", 147, "chr1", 5200, "50M", "=", 5000, -250)
        .write_to(&env.input);

    let r = dupblaster()
        .args(["-i"])
        .arg(&env.input)
        .args(["-o"])
        .arg(&out_bam)
        .args(["--metrics-prefix"])
        .arg(metrics_prefix_for(&out_bam))
        .output()
        .unwrap();
    assert!(r.status.success(), "stderr: {}", String::from_utf8_lossy(&r.stderr));
    assert!(!String::from_utf8_lossy(&r.stderr).contains("clamped"));
}