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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

//! Generative tests for dupblaster's pair-level duplicate semantics:
//! ordering invariance, three-deep duplicate sets, --remove-dups, and
//! cross-chromosome signature uniqueness.

mod helpers;

use helpers::*;

#[test]
fn pair_order_swap_does_not_change_duplicate_decision() {
    let env = TestEnv::new();
    // Pair A: R1 first, then R2. Pair B: R2 first, then R1. Same
    // coordinates and strands — the signature must canonicalize so
    // both pairs are detected as duplicates regardless of arrival order.
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("rA", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("rA", 147, "chr1", 200, "50M", "=", 100, -150)
        // rB written R2 first, then R1.
        .rec_simple("rB", 147, "chr1", 200, "50M", "=", 100, -150)
        .rec_simple("rB", 99, "chr1", 100, "50M", "=", 200, 150)
        .write_to(&env.input);
    let bam_out = env._tmp.path().join("out.bam");
    let recs = run_and_capture(&env.input, &bam_out, &[]).flags;

    assert_eq!(recs[0], ("rA".into(), 99));
    assert_eq!(recs[1], ("rA".into(), 147));
    // Output order matches input order: rB's R2 came first, then R1.
    assert_eq!(recs[2], ("rB".into(), 147 | FLAG_DUPLICATE));
    assert_eq!(recs[3], ("rB".into(), 99 | FLAG_DUPLICATE));
}

#[test]
fn three_identical_pairs_only_first_kept_as_non_duplicate() {
    let env = TestEnv::new();
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("r1", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r1", 147, "chr1", 200, "50M", "=", 100, -150)
        .rec_simple("r2", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r2", 147, "chr1", 200, "50M", "=", 100, -150)
        .rec_simple("r3", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r3", 147, "chr1", 200, "50M", "=", 100, -150)
        .write_to(&env.input);
    let bam_out = env._tmp.path().join("out.bam");
    let recs = run_and_capture(&env.input, &bam_out, &[]).flags;

    // r1 is the representative (no dup flag); r2 and r3 are duplicates.
    assert_eq!(recs[0], ("r1".into(), 99));
    assert_eq!(recs[1], ("r1".into(), 147));
    assert_eq!(recs[2], ("r2".into(), 99 | FLAG_DUPLICATE));
    assert_eq!(recs[3], ("r2".into(), 147 | FLAG_DUPLICATE));
    assert_eq!(recs[4], ("r3".into(), 99 | FLAG_DUPLICATE));
    assert_eq!(recs[5], ("r3".into(), 147 | FLAG_DUPLICATE));
}

#[test]
fn remove_dups_drops_duplicate_records_from_output() {
    let env = TestEnv::new();
    // Three identical pairs. With -r, only r1's two records should remain.
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .rec_simple("r1", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r1", 147, "chr1", 200, "50M", "=", 100, -150)
        .rec_simple("r2", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r2", 147, "chr1", 200, "50M", "=", 100, -150)
        .rec_simple("r3", 99, "chr1", 100, "50M", "=", 200, 150)
        .rec_simple("r3", 147, "chr1", 200, "50M", "=", 100, -150)
        .write_to(&env.input);
    let bam_out = env._tmp.path().join("out.bam");
    let recs = run_and_capture(&env.input, &bam_out, &["-r"]).flags;

    // Only the kept pair survives. No record carries FLAG_DUPLICATE
    // because the dup records have been removed.
    assert_eq!(recs.len(), 2);
    assert_eq!(recs[0], ("r1".into(), 99));
    assert_eq!(recs[1], ("r1".into(), 147));
}

#[test]
fn cross_chromosome_pairs_at_distinct_positions_are_not_dups() {
    let env = TestEnv::new();
    // Two distinct cross-chrom pairs. The two-key signature includes the
    // chromosome of each end, so different (chr,pos) combinations must
    // not collide.
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .sq("chr2", 1_000_000)
        // Pair A: R1 on chr1:100 fwd, R2 on chr2:200 rev.
        .rec_simple("rA", 97, "chr1", 100, "50M", "chr2", 200, 0)
        .rec_simple("rA", 145, "chr2", 200, "50M", "chr1", 100, 0)
        // Pair B: R1 on chr1:300 fwd, R2 on chr2:400 rev — different positions.
        .rec_simple("rB", 97, "chr1", 300, "50M", "chr2", 400, 0)
        .rec_simple("rB", 145, "chr2", 400, "50M", "chr1", 300, 0)
        .write_to(&env.input);
    let bam_out = env._tmp.path().join("out.bam");
    let recs = run_and_capture(&env.input, &bam_out, &[]).flags;

    for (qname, flag) in &recs {
        assert_eq!(
            flag & FLAG_DUPLICATE,
            0,
            "{qname} (flag {flag}) should not be marked duplicate"
        );
    }
}

#[test]
fn cross_chromosome_pairs_at_identical_positions_are_dups() {
    let env = TestEnv::new();
    // Same as the "distinct positions" test but pair B mirrors A's
    // coordinates exactly — must be detected as a dup.
    SamBuilder::new()
        .sq("chr1", 1_000_000)
        .sq("chr2", 1_000_000)
        .rec_simple("rA", 97, "chr1", 100, "50M", "chr2", 200, 0)
        .rec_simple("rA", 145, "chr2", 200, "50M", "chr1", 100, 0)
        .rec_simple("rB", 97, "chr1", 100, "50M", "chr2", 200, 0)
        .rec_simple("rB", 145, "chr2", 200, "50M", "chr1", 100, 0)
        .write_to(&env.input);
    let bam_out = env._tmp.path().join("out.bam");
    let recs = run_and_capture(&env.input, &bam_out, &[]).flags;

    assert_eq!(recs[0], ("rA".into(), 97));
    assert_eq!(recs[1], ("rA".into(), 145));
    assert_eq!(recs[2], ("rB".into(), 97 | FLAG_DUPLICATE));
    assert_eq!(recs[3], ("rB".into(), 145 | FLAG_DUPLICATE));
}