[BenjaminDEMAILLE]

Part of #354 (item 25).

Not stacked; based on main.

Problem

bam_remove_duplicates read every record of the input into memory, then grouped clusters, then wrote everything. On a 7.6 million record BAM (175 MB) that is 6.4 GB of RSS, 36 times the file and about four times what STAR uses on the same input.

Change

The first loop now does all of it. When a record falls outside the current cluster (another reference, or past the cluster's last mate), write_cluster collapses that cluster, writes all its records in input order and clears the buffers. The tests for a boundary are the same as before, against the same values, applied to the record as it arrives instead of in a second pass, and every record since the cluster began is written at that point, multimappers included. So the bytes do not change.

collapse_group is untouched. DedupRec loses pos, nh and mate_pos, which only the second pass read.

What changes for a user

• Memory follows the largest cluster, not the file. A cluster only ends past a mate on a later coordinate. If no record of the cluster has one, it runs to the end of the reference. So a single-end BAM still holds its largest chromosome in memory, measured below. That is how clusters are defined in the existing code, which I kept in order to keep the output.
• A run that fails leaves more behind. A missing NH or AS tag is still a fatal error with the same message, and exit status is still non-zero. But Processed.out.bam now holds the clusters that ended before the bad record, and it looks finished: on a test BAM with a bad 11th record, main left 0 bytes and this branch leaves 8 records that pass samtools quickcheck, both with exit status 1. A pipeline that ignores the exit status would take it for a complete file. Removing the file on error would be a small change, if you prefer that; I kept this PR to what was asked. A file with several bad records can report a different one first: NH used to be checked on the whole file before any AS.

Measured

Same machine, same input, no --runThreadN. macOS arm64, /usr/bin/time -l, maximum resident set size. The machine was loaded (load average 65 to 120), so wall-clock numbers are only a rough guide; I make no claim about speed.

main and branch are binaries built from 32fbbc1 and from this branch.

Peak RSS

Input	Records	BAM size	main	This branch	STAR 2.7.11b
generated PE, seed 1	76,118	1.8 MB	89.5 MB	9.8 MB	20.5 MB
generated PE, seed 2	77,136	1.8 MB	89.6 MB	10.1 MB	20.8 MB
generated PE, seed 3	7,599,752	175 MB	6,370.7 MB	10.3 MB	1,644.6 MB
same generator, first mates only, single-end	1,898,073	38 MB	1,510.2 MB	1,123.6 MB	not run

The numbers are for UniqueIdentical; UniqueIdenticalNotMulti differs by less than 0.4 MB on every paired-end row. Wall clock on the 7.6M row was 13.2 s against 12.8 s (UniqueIdentical) and 20.1 s against 18.7 s (UniqueIdenticalNotMulti).

Byte-identical to main

samtools view -h --no-PG | md5. --no-PG because without it samtools adds its own command line, with the file path, to the header. The file md5 column is the md5 of Processed.out.bam itself, which is identical too, compression included.

Input	Type	main	This branch	File md5 same
seed 1	UniqueIdentical	125a2a220682	125a2a220682	yes
seed 1	UniqueIdenticalNotMulti	9199a8cdba8b	9199a8cdba8b	yes
seed 2	UniqueIdentical	fb291205bbe8	fb291205bbe8	yes
seed 2	UniqueIdenticalNotMulti	c7f5962c7777	c7f5962c7777	yes
seed 3	UniqueIdentical	fd4b061a3090	fd4b061a3090	yes
seed 3	UniqueIdenticalNotMulti	732ef6625816	732ef6625816	yes
STAR-aligned PE (below)	UniqueIdentical	e74c065bb5ec	e74c065bb5ec	yes
STAR-aligned PE	UniqueIdenticalNotMulti	a7ee2b9a88d0	a7ee2b9a88d0	yes
STAR-aligned PE	UniqueIdentical, --bamRemoveDuplicatesMate2basesN 20	71581f5b4d00	71581f5b4d00	yes
single-end	UniqueIdentical	file md5 1dfa11a08741	1dfa11a08741	yes

Against STAR 2.7.11b (--runMode inputAlignmentsFromBAM)

Observed output, source not read. Homebrew's STAR. This is the same on main: this PR changes none of it.

Input	Type	Duplicates, branch / STAR	Header	Records identical to STAR
seed 1	UniqueIdentical	39,322 / 39,322	same	yes, whole view -h md5 125a2a220682
seed 1	UniqueIdenticalNotMulti	36,120 / 36,120	same	yes, 9199a8cdba8b
seed 2	UniqueIdentical	40,248 / 40,248	same	no, 1 of 18,444 alignment keys keeps another pair
seed 2	UniqueIdenticalNotMulti	37,138 / 37,138	same	no
seed 3	UniqueIdentical	3,920,410 / 3,920,410	same	no, 2,675 of 1,839,671 keys
seed 3	UniqueIdenticalNotMulti	3,599,762 / 3,599,762	same	no
STAR-aligned PE	UniqueIdentical	9,506 / 9,506	same	no, 1,202 of 7,996 keys
STAR-aligned PE	UniqueIdenticalNotMulti	9,498 / 9,498	same	no, 1,202 of 7,996 keys
STAR-aligned PE	UniqueIdentical, --bamRemoveDuplicatesMate2basesN 20	9,020 / 9,020	same	no

The duplicate counts are always equal, and the differences are of one kind. I grouped the pairs by the dedup key (both mates' flag, soft-clip-folded position and CIGAR) and compared the pair each tool kept: in every group that differs, in the four sets checked (seed 2 and seed 3 with UniqueIdentical, STAR-aligned with both types), the two kept pairs have the same AS. So rustar-aligner and STAR agree on what is a duplicate and on the best score, and differ on which of several equally scored pairs survives. rustar-aligner takes the first in a stable sort. STAR's choice looks like an artefact of its sort, but that is a guess from the output; I have not read the source and did not find out what order it uses. This is an existing parity gap and is not part of this PR. The mate-2 row was not put through the grouping check, because that key would need the mate's SEQ as well.

The STAR-aligned BAM: 8,000 simulated templates (1 to 3 copies each, 20 % of reads with one changed base) on a random 20 kb genome, aligned by STAR with --outSAMtype BAM SortedByCoordinate --outFilterMismatchNmax 3: 25,498 records, 8 of them multimappers.

Reproducing the generated BAMs

python3 bam_gen.py 20000 1 | samtools view -b -o small.bam -
python3 bam_gen.py 20000 2 | samtools view -b -o small2.bam -
python3 bam_gen.py 2000000 3 | samtools view -b -o big.bam -
/usr/bin/time -l rustar-aligner --runMode inputAlignmentsFromBAM --inputBAMfile big.bam \
  --bamRemoveDuplicatesType UniqueIdentical --outFileNamePrefix out/

bam_gen.py: coordinate-sorted paired-end SAM with PCR duplicates, soft clips, splices and multimappers

#!/usr/bin/env python3
"""Write a coordinate-sorted, paired-end SAM to stdout: unique pairs with PCR
duplicates (same or different AS), soft clips, splices, chains of overlapping
clusters, multimappers, half-mapped and unpaired records.

usage: bam_gen.py N_TEMPLATES SEED
"""
import heapq, random, sys

n, seed = int(sys.argv[1]), int(sys.argv[2])
R = random.Random(seed)
CHRS = [('chr1', 200_000_000), ('chr2', 100_000_000), ('chr3', 50_000_000)]
L = 100
QUAL = 'I' * L
BASES = [''.join(R.choices('ACGT', k=L)) for _ in range(64)]


def line(name, flag, chrom, pos, cigar, mchr, mpos, tlen, nh, score):
    return (pos, f'{name}\t{flag}\t{chrom}\t{pos}\t{255 if nh == 1 else 3}\t{cigar}\t{mchr}\t{mpos}\t{tlen}\t'
                 f'{R.choice(BASES)}\t{QUAL}\tNH:i:{nh}\tHI:i:1\tAS:i:{score}\tnM:i:0\n')


def cigar():
    x = R.random()
    if x < 0.8:
        return f'{L}M', L
    if x < 0.9:
        s = R.randint(1, 20)
        return f'{s}S{L - s}M', L - s
    if x < 0.95:
        s = R.randint(1, 20)
        return f'{L - s}M{s}S', L - s
    a, gap = R.randint(10, 90), R.randint(50, 5000)
    return f'{a}M{gap}N{L - a}M', L + gap


print('@HD\tVN:1.6\tSO:coordinate')
for c, ln in CHRS:
    print(f'@SQ\tSN:{c}\tLN:{ln}')
print('@PG\tID:bam_gen\tPN:bam_gen')

total = sum(ln for _, ln in CHRS)
serial = 0
for c, ln in CHRS:
    per_chr = int(n * ln / total)
    # Starts in order, with gaps so that clusters end now and then.
    starts = sorted(R.randrange(1, ln - 20_000) for _ in range(per_chr))
    out = []
    for start in starts:
        serial += 1
        kind = R.random()
        if kind < 0.08:  # a multimapper pair
            cig, span = cigar()
            m = start + R.randint(0, 400)
            h = R.randint(2, 5)
            out.append(line(f't{serial}', 99 | 256 * R.randint(0, 1), c, start, cig, '=', m, m + L - start, h, 150))
            out.append(line(f't{serial}', 147, c, m, f'{L}M', '=', start, start - m - L, h, 150))
            continue
        copies = 1 if kind < 0.7 else R.randint(2, 6)  # PCR duplicates
        cig, span = cigar()
        mate_gap = R.randint(0, 400)
        m = start + mate_gap
        rev = R.random() < 0.5
        f1, f2 = (99, 147) if not rev else (83, 163)
        mcig, _ = cigar()
        for copy in range(copies):
            serial += 1
            score = R.randint(140, 198)
            # Some copies differ by a clip only: same key once S is folded.
            out.append(line(f't{serial}', f1, c, start, cig, '=', m, m + L - start, 1, score))
            out.append(line(f't{serial}', f2, c, m, mcig, '=', start, start - m - L, 1, score))
    out.sort(key=lambda t: t[0])
    sys.stdout.writelines(t[1] for t in out)

The single-end input is this with seed 4 and 1,000,000 templates, keeping lines without the 128 flag bit and setting FLAG to 0 or 16, RNEXT to *, and PNEXT and TLEN to 0.

Tests

In src/bam_dedup.rs. There was no test of bam_remove_duplicates as a whole before this.

• clusters_are_collapsed_and_written_in_input_order: an in-memory BAM with two duplicate pairs and a multimapper in one cluster, a later cluster, and the same coordinates on another reference. It checks each record's flag and the order, with mark_multi on and off.
• a_missing_tag_leaves_the_earlier_clusters_written: the error is returned, and the output holds the cluster that had ended and not the one in progress. This is also the test that fails if the function goes back to reading everything first.

cargo test --release (815 passed), cargo clippy --all-targets and cargo fmt --check are clean on macOS.

Not checked

• test/parity_suite.sh was not run. It does not cover this run mode as far as I can tell from grep: nothing under tests/ or test/ passes --bamRemoveDuplicatesType.
• No real sequencing data: the BAMs are generated, or aligned from simulated reads.
• The stdin path (--inputBAMfile -) was not run end to end. It calls the same function.
• A real paired-end BAM can have larger clusters than these. A long chain of overlapping pairs, as in a very highly expressed gene, is one cluster, and its memory is what remains. I did not measure such a case.
• Single-end RSS against STAR was not measured.

🤖 Generated with Claude Code