Skip to content

bamRemoveDuplicates: a mateless unique read keeps the cluster open, as STAR - #370

Open
BenjaminDEMAILLE wants to merge 1 commit into
perf/bam-dedup-streamingfrom
fix/bam-dedup-ties
Open

BenjaminDEMAILLE wants to merge 1 commit into
perf/bam-dedup-streamingfrom
fix/bam-dedup-ties

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

Part of #354 (--bamRemoveDuplicatesType: which pair STAR keeps).

PR base: perf/bam-dedup-streaming (#363). That PR rewrote src/bam_dedup.rs to stream, so this fix is on the new code. Please merge #363 first; only the last commit is this PR.

Problem

The task was: rustar keeps another pair than STAR among duplicates with the same AS. I tested that against STAR 2.7.11b built on Linux, and ties already match. The tie difference only shows up against the macOS build of STAR (see below). What did differ is something else: a unique read without a mate (single-end, PNEXT 0 in SAM, -1 in BAM) next to paired reads. STAR keeps it; rustar marked it as a duplicate.

What STAR does (2.7.11b source)

All in bamRemoveDuplicates.cpp.

  • Which pair is kept. Unique records of a cluster are sorted by name with qsort (line 206; funCompareNames, line 13: name length, then bytes, then the 0x80 bit), then pairs are sorted with qsort (line 207) by funCompareCoordFlagCigarSeq (line 72). In each run of equal pairs, the first pair with a strictly higher AS wins (lines 209, 223-226). So on a tie, the winner is the first pair in the sorted order.
  • glibc qsort is a stable merge sort (the Linux image runs glibc 2.36). So equal pairs stay in name order, and the tie goes to the smallest name. rustar's sort_by is stable too, so it gives the same answer.
  • macOS libc qsort is not stable. /opt/homebrew/bin/STAR keeps another pair on ties (for example hi instead of b9, both AS 30, same key). I took the Linux build as the reference.
  • Clusters. leftE, rightE and rightMax are uint32 (lines 141, 173-174, 199-200). A cluster ends when leftE > rightMax (line 205). A unique read with rightE > leftE raises rightMax (line 258). With PNEXT -1, rightE is 0xFFFFFFFF, so rightMax becomes 0xFFFFFFFF and the cluster runs to the end of the reference.
  • Pairs are taken two by two after the name sort (line 207, grN/2 pairs), so a mateless read pairs with whatever comes next by name.

Change

  • src/bam_dedup.rs: compare POS/PNEXT as u32 when building clusters, as STAR does. Three lines.
  • Module doc now says how ties are broken and that the macOS STAR build can differ.
  • Two unit tests (below) and a CHANGELOG entry.

Measured against STAR 2.7.11b

STAR built for Linux (Docker, glibc 2.36). Each BAM was run with UniqueIdentical and UniqueIdenticalNotMulti. I compared samtools view md5. For sorted1 I also compared the md5 of the decompressed BAM (gzcat): identical. Records that differ = records whose name/flag/chr/pos differ from STAR.

Input Records Before (both types) After (both types)
case1: AS ties, names out of order and of different lengths, a multimapper, 2S18M same key, other flags, 2 refs 26 identical identical
case2: pairs plus 2 mateless unique reads 14 2 records differ identical
case3: single-end only, ties, a multimapper 6 identical identical
sorted1 / sorted2 / sorted3 (paired, many ties) 112183 / 112507 / 112752 identical identical
se1 / se2 / se3 (single-end, many ties) 49855 / 49748 / 50182 identical identical
mixed (sorted1 + se1 merged) 162038 identical identical
big (paired) 5633297 identical identical

case2, before:

STAR:    se1  0     chr1  150      se2  0     chr1  1050
rustar:  se1  1024  chr1  150      se2  1024  chr1  1050

STAR has one cluster for all of chr1, so se1 and se2 pair up by name and are kept. rustar split the cluster at the gap, so each of them was left alone and stayed marked.

Tests

  • a_mateless_unique_read_keeps_the_cluster_open: the case2 shape. It fails without the fix ([99, 1024, 147, 99, 1024, 147]) and passes with it.
  • an_as_tie_keeps_the_first_pair_by_name: pins the tie rule (b9 wins over alpha and zeta on the same AS). It passes before and after.
  • cargo test --release, cargo clippy --all-targets -- -D warnings and cargo fmt --check are clean on macOS.

Not checked

  • Names with bytes above 127: STAR compares char, which is signed on x86. rustar compares unsigned bytes. SAM does not allow such names, so I did not change this.
  • AS of -999 or lower: STAR starts bScore at -999, so such a run keeps a stale pair. rustar starts at the minimum. Not matched.
  • Only the macOS binary of rustar was run; STAR ran on Linux in Docker.

🤖 Generated with Claude Code

STAR compares POS and PNEXT as uint32 when it builds the clusters, so a
unique read without a mate (PNEXT -1) sets the cluster's right end to
0xFFFFFFFF and holds it open to the end of the reference. rustar compared
them as signed, ended the cluster earlier, and could leave such a read
unpaired and marked where STAR keeps it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant