From da5432f65ec4ca580f171bb0393b0798849dee71 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 00:38:15 +0200 Subject: [PATCH] fix(solo): port STAR's actual CellRanger4 clip rules Both ends of --clipAdapterType CellRanger4 were approximations of rules STAR implements exactly, so the clip lengths differed from STAR's in both directions. 5' TSO. STAR aligns the TSO against the read with Opal in overlap mode (ClipCR4.cpp: match +1, mismatch -2, gap open and extend 2) and clips endLocationTarget + 1 bases unless the alignment is too weak (ClipMate_clipChunk.cpp: S < 20, or S == 20 with L > 26, or S == 21 with L > 30). Here it was a fixed-length prefix comparison with a mismatch budget, which cannot see a TSO that starts a few bases into the read and fires on full-length matches STAR scores below the floor. Ported as an affine-gap DP: the query is 30 bases and the window is 91, so it is a small fixed cost. 3' poly-A. Ported ClipCR4::polyTail3p verbatim: +1 per A, -2 otherwise, remember the longest prefix scoring at least 70% of its length, stop once the score falls more than 27 behind, and require a final score of 20. The previous rule trimmed only a literal run of A of length >= 8, so a tail with one sequencing error kept about half its length in the alignment. Also adds test/cr4_clip_diff.py, a synthetic STAR-vs-rustar differential for this flag combination that runs in seconds instead of needing the 10x mouse chr19 dataset, and ignores the local .claude/ and *.code-workspace scratch files that keep landing in `git add -A`. --- .gitignore | 4 + CHANGELOG.md | 9 ++ src/solo/mod.rs | 274 +++++++++++++++++++++++++++++++++++++----- test/cr4_clip_diff.py | 156 ++++++++++++++++++++++++ 4 files changed, 412 insertions(+), 31 deletions(-) create mode 100644 test/cr4_clip_diff.py diff --git a/.gitignore b/.gitignore index 042ce739..3715a1fa 100644 --- a/.gitignore +++ b/.gitignore @@ -47,3 +47,7 @@ TODO.md # macOS Finder metadata .DS_Store + +# Local editor/agent scratch +.claude/ +*.code-workspace diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346b..c9890bd7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,15 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `--clipAdapterType CellRanger4` now uses STAR's own two rules instead of + approximations. The 5' TSO clip is an overlap alignment (`ClipCR4.cpp` + + `ClipMate_clipChunk.cpp`: match +1, mismatch -2, gap open/extend 2, keep + `endLocationTarget + 1` unless the score is below the floor), where it was a + fixed-length prefix comparison that could neither see a TSO starting partway + into the read nor reject a full-length match STAR scores too low. The 3' + poly-A trim is a port of `ClipCR4::polyTail3p`, a scored scan that walks + through sequencing errors, where it was a run of literal `A`s and so kept + roughly half of any tail carrying one error. - Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`. diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e02..7b4f8de8 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -391,44 +391,155 @@ pub fn open_paired_reader(params: &Parameters) -> Result i32 { + if a == 4 && b == 4 { + 0 + } else if a == b { + 1 + } else { + -2 + } +} + +/// How many 5' bases the CR4 TSO clip removes from `seq`. +/// +/// STAR aligns the TSO against the read with Opal in `OPAL_MODE_OV` (overlap: +/// gaps before the start and after the end of either sequence are free), then +/// keeps `endLocationTarget + 1` bases as the clip unless the alignment is too +/// weak (`ClipMate_clipChunk.cpp`): +/// +/// ```text +/// L0 = S < 20 || (S == 20 && L > 26) || (S == 21 && L > 30) +/// ``` +/// +/// A fixed-length prefix comparison is not the same rule: it misses a TSO that +/// starts partway into the read or carries an indel, and it fires on a +/// full-length match whose score STAR would reject. Both directions were +/// visible as reads misplaced by a few bases (#199). +pub fn cr4_tso_clip_len(seq: &[u8]) -> usize { + let query: Vec = TSO_SEQ + .iter() + .map(|&b| crate::io::fastq::encode_base(b)) + .collect(); + + // The target is the read padded with N to STAR's fixed window. + let mut target: Vec = Vec::with_capacity(CR4_OPAL_READ_LEN); + target.extend(seq.iter().copied().take(CR4_OPAL_READ_LEN)); + target.resize(CR4_OPAL_READ_LEN, 4); + + let m = query.len(); + let n = target.len(); + + // Affine-gap DP. Row 0 and column 0 are zero: leading gaps are free on + // both sequences, which is what OV mode means. + const NEG: i32 = i32::MIN / 4; + let mut h_prev = vec![0i32; n + 1]; + let mut e_prev = vec![NEG; n + 1]; + let mut h_cur = vec![0i32; n + 1]; + let mut e_cur = vec![NEG; n + 1]; + + // Best over the last row and the last column: trailing gaps are free too. + let mut best_score = i32::MIN; + let mut best_end = 0usize; + + for i in 1..=m { + h_cur[0] = 0; + e_cur[0] = NEG; + let mut f = NEG; // gap in the query (consuming target bases) + for j in 1..=n { + // Gap in the target (consuming query bases). + e_cur[j] = (h_prev[j] - CR4_GAP_OPEN).max(e_prev[j] - CR4_GAP_EXT); + f = (h_cur[j - 1] - CR4_GAP_OPEN).max(f - CR4_GAP_EXT); + let diag = h_prev[j - 1] + cr4_score(query[i - 1], target[j - 1]); + h_cur[j] = diag.max(e_cur[j]).max(f); + + // Last column: the query may end early with a free trailing gap. + if j == n && h_cur[j] > best_score { + best_score = h_cur[j]; + best_end = j - 1; + } + } + if i == m { + // Last row: every end position in the target is a candidate. + for (j, &h) in h_cur.iter().enumerate().skip(1) { + if h > best_score { + best_score = h; + best_end = j - 1; + } + } + } + std::mem::swap(&mut h_prev, &mut h_cur); + std::mem::swap(&mut e_prev, &mut e_cur); + } + + let l = best_end + 1; + let s = best_score; + let too_weak = s < 20 || (s == 20 && l > 26) || (s == 21 && l > 30); + if too_weak { 0 } else { l.min(seq.len()) } +} + +/// How many 3' bases the CR4 poly-A trim removes from `seq`. +/// +/// Direct port of `ClipCR4::polyTail3p`: walking in from the 3' end, an `A` +/// scores +1 and anything else -2; the longest prefix of that walk whose score +/// is at least 70% of its length is remembered, the scan stops once the score +/// has dropped more than 27 behind the length, and a final score below 20 +/// means no trim at all. Reads shorter than 20 bases are never trimmed. +/// +/// A plain "trailing run of A" rule is stricter: it stops at the first +/// non-`A`, so a tail with one sequencing error keeps ~half its length. +pub fn cr4_polya_clip_len(seq: &[u8]) -> usize { + let seq_len = seq.len(); + if seq_len < 20 { + return 0; + } + let mut ib1 = seq_len - 1; + let mut score: i32 = 0; + let mut score1: i32 = 0; + for ib in 1..=seq_len { + if seq[seq_len - ib] == 0 { + score += 1; + if score * 10 >= (ib as i32) * 7 { + ib1 = ib; + score1 = score; + } + } else { + score -= 2; + if (ib as i32) - score > 27 { + break; + } + } + } + if score1 < 20 { 0 } else { ib1 } +} + /// Clip the 10x TSO from the 5' end and trim a 3' polyA tail of the cDNA read, /// matching `--clipAdapterType CellRanger4`. Operates on encoded bases -/// (0=A..3=T,4=N) with parallel quality bytes. Returns the clipped read. +/// (0=A..3=T,4=N) with parallel quality bytes. +/// +/// Both ends follow STAR's own rules: see [`cr4_tso_clip_len`] and +/// [`cr4_polya_clip_len`]. /// -/// Conservative thresholds (full-length TSO match ≤ 3 mismatches at the 5' -/// anchor; trailing polyA run ≥ 8) keep this a no-op on adapter-free reads. /// Returns `(clipped_seq, clipped_qual, clip5p, clip3p)` — the CR4-clipped read plus /// the bases trimmed from the 5' (TSO) and 3' (polyA) ends, so the caller can soft-clip /// them (STARsolo keeps them in SEQ as soft-clips, e.g. `60M30S`, not dropped). pub fn clip_adapter_cr4(seq: &[u8], qual: &[u8]) -> (Vec, Vec, usize, usize) { - let mut start = 0usize; - let mut end = seq.len(); - - // 5' TSO: compare the read prefix against the full TSO; clip on a match. - if seq.len() >= TSO_SEQ.len() { - let tso: Vec = TSO_SEQ - .iter() - .map(|&b| crate::io::fastq::encode_base(b)) - .collect(); - let mismatches = seq[..tso.len()] - .iter() - .zip(&tso) - .filter(|(a, b)| a != b) - .count(); - if mismatches <= 3 { - start = tso.len(); - } - } - - // 3' polyA: trim a trailing run of A (encoded 0) of length >= 8. - let mut run = 0usize; - while end > start && seq[end - 1] == 0 { - run += 1; - end -= 1; - } - if run < 8 { - end += run; // not a real polyA tail; keep those bases - } + let start = cr4_tso_clip_len(seq); + // STAR trims the poly-A from the read that is left after the 5' clip. + let after5p = &seq[start..]; + let trim3 = cr4_polya_clip_len(after5p); + let end = seq.len() - trim3; if start == 0 && end == seq.len() { return (seq.to_vec(), qual.to_vec(), 0, 0); @@ -443,6 +554,107 @@ pub fn clip_adapter_cr4(seq: &[u8], qual: &[u8]) -> (Vec, Vec, usize, us ) } +#[cfg(test)] +mod cr4_clip_tests { + use super::{TSO_SEQ, cr4_polya_clip_len, cr4_tso_clip_len}; + use crate::io::fastq::encode_base; + + fn enc(s: &[u8]) -> Vec { + s.iter().map(|&b| encode_base(b)).collect() + } + + /// A read that is nothing but cDNA must not be clipped: the whole point of + /// the score floor is that a weak, incidental match is not an adapter. + #[test] + fn a_read_without_the_tso_is_not_clipped() { + let read = b"TTTTGCACTGCACGTGTCGATCGGCATCGGATCGATCGGCATTTACGCTACGTACGATCGATCGGCATCGATCGATCGTACGGCATCGAT"; + assert_eq!(cr4_tso_clip_len(&enc(read)), 0); + } + + /// The full TSO at the very start is clipped exactly, and nothing more. + #[test] + fn a_full_tso_prefix_is_clipped_to_its_length() { + let mut read = TSO_SEQ.to_vec(); + read.extend_from_slice(b"GCACTGCACGTGTCGATCGGCATCGGATCGATCGGCATTTACGCTACGTACGATCGATCGG"); + assert_eq!(cr4_tso_clip_len(&enc(&read)), TSO_SEQ.len()); + } + + /// A short TSO overlap scores below STAR's floor of 20, so it stays. + /// This is the case a fixed-prefix comparison and STAR's rule agree on, + /// and it is why the fix does not simply clip more. + #[test] + fn a_five_base_tso_overlap_is_below_the_score_floor() { + let mut read = TSO_SEQ[TSO_SEQ.len() - 5..].to_vec(); + read.extend_from_slice( + b"GCACTGCACGTGTCGATCGGCATCGGATCGATCGGCATTTACGCTACGTACGATCGATCGGCATCGAT", + ); + assert_eq!(cr4_tso_clip_len(&enc(&read)), 0); + } + + /// A TSO carrying mismatches is still found, and the clip covers it: the + /// score is 30 - 3*3 = 21 with L == 30, which STAR keeps. + #[test] + fn a_tso_with_three_mismatches_is_still_clipped() { + let mut tso = TSO_SEQ.to_vec(); + for i in [3usize, 11, 19] { + tso[i] = if tso[i] == b'A' { b'C' } else { b'A' }; + } + let mut read = tso; + read.extend_from_slice(b"GCACTGCACGTGTCGATCGGCATCGGATCGATCGGCATTTACGCTACGTACGATCGATCGG"); + assert_eq!(cr4_tso_clip_len(&enc(&read)), TSO_SEQ.len()); + } + + /// The clip is an overlap alignment, so a TSO that begins a few bases into + /// the read is clipped up to its end, not missed. A prefix comparison + /// anchored at position 0 cannot see this at all. + #[test] + fn a_tso_starting_inside_the_read_is_clipped_through_its_end() { + let mut read = b"GATC".to_vec(); + read.extend_from_slice(TSO_SEQ); + read.extend_from_slice(b"GCACTGCACGTGTCGATCGGCATCGGATCGATCGGCATTTACGCTACGTACGATCG"); + assert_eq!(cr4_tso_clip_len(&enc(&read)), 4 + TSO_SEQ.len()); + } + + /// Poly-A: a clean tail is trimmed, and only the tail. The cDNA before it + /// deliberately ends in non-`A` bases, because STAR's scan keeps walking + /// past the tail and will absorb an `A` a base or two upstream. + #[test] + fn a_clean_polya_tail_is_trimmed() { + let mut read = b"GCTCTGCTCGTGTCGCTCGGCTTCGGCTCGCTCGGCTTTTCGCTCCGTTCG".to_vec(); + let tail = 30; + read.extend(std::iter::repeat_n(b'A', tail)); + assert_eq!(cr4_polya_clip_len(&enc(&read)), tail); + } + + /// A tail with one sequencing error keeps its full length: STAR's scan + /// pays -2 for the mismatch and carries on, where a "trailing run of A" + /// rule would stop at the error and keep half the tail in the alignment. + #[test] + fn a_polya_tail_with_one_error_is_still_trimmed_whole() { + let mut read = b"GCTCTGCTCGTGTCGCTCGGCTTCGGCTCGCTCGGCTTTTCGCTCCGTTCG".to_vec(); + let mut tail = vec![b'A'; 30]; + tail[10] = b'G'; + read.extend_from_slice(&tail); + assert_eq!(cr4_polya_clip_len(&enc(&read)), 30); + } + + /// A short run of A is ordinary sequence, not a tail: the final score has + /// to reach 20. + #[test] + fn a_short_a_run_is_not_a_tail() { + let mut read = b"GCTCTGCTCGTGTCGCTCGGCTTCGGCTCGCTCGGCTTTTCGCTCCGTTCG".to_vec(); + read.extend_from_slice(b"AAAAAAAA"); + assert_eq!(cr4_polya_clip_len(&enc(&read)), 0); + } + + /// Reads under 20 bases are never trimmed (`ClipCR4::polyTail3p`). + #[test] + fn a_very_short_read_is_never_trimmed() { + let read = vec![b'A'; 19]; + assert_eq!(cr4_polya_clip_len(&enc(&read)), 0); + } +} + // --------------------------------------------------------------------------- // Solo counting context + per-read processing (Phase 14.3) // --------------------------------------------------------------------------- diff --git a/test/cr4_clip_diff.py b/test/cr4_clip_diff.py new file mode 100644 index 00000000..26a6eac6 --- /dev/null +++ b/test/cr4_clip_diff.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 +"""Differential check for --clipAdapterType CellRanger4 + --clip5pNbases. + +Builds a small synthetic solo fixture, runs STAR 2.7.11b and rustar-aligner +with the same flags, and reports the per-read POS and leading-soft-clip +differences. This is the measurement behind issue #199, shrunk to something +that runs in seconds instead of needing the 10x mouse chr19 dataset. +""" +import os +import random +import subprocess +import sys +from pathlib import Path + +TSO = "AAGCAGTGGTATCAACGCAGAGTACATGGG" +CB_LEN, UMI_LEN = 16, 12 +READ_LEN = 90 +N_READS = 400 + +RUSTAR = sys.argv[1] if len(sys.argv) > 1 else "./target/release/rustar-aligner" +OUT = Path(sys.argv[2] if len(sys.argv) > 2 else "/tmp/cr4diff") + + +def lcg(seed, n): + bases = "ACGT" + state = seed + out = [] + for _ in range(n): + state = (state * 1103515245 + 12345) & 0xFFFFFFFF + out.append(bases[(state >> 16) & 3]) + return "".join(out) + + +def main(): + rng = random.Random(20260827) + OUT.mkdir(parents=True, exist_ok=True) + genome = lcg(88888, 20000) + (OUT / "genome.fa").write_text(">chr1\n" + genome + "\n") + + # A minimal annotation: one long exon, so gene assignment never filters. + (OUT / "genes.gtf").write_text( + 'chr1\tsyn\texon\t1\t20000\t.\t+\t.\tgene_id "G1"; transcript_id "G1_T1";\n' + ) + + cdna, barcode, whitelist = [], [], [] + for i in range(N_READS): + start = rng.randrange(200, 19000 - READ_LEN) + body = genome[start : start + READ_LEN] + # Four kinds of read: clean, full TSO prefix, TSO a few bases in, + # TSO with mismatches. Each exercises a different branch of the rule. + kind = i % 4 + if kind == 0: + seq = body + elif kind == 1: + seq = (TSO + body)[:READ_LEN] + elif kind == 2: + seq = ("GATC" + TSO + body)[:READ_LEN] + else: + tso = list(TSO) + for p in (3, 11, 19): + tso[p] = "C" if tso[p] == "A" else "A" + seq = ("".join(tso) + body)[:READ_LEN] + cdna.append((f"r{i}", seq)) + cb = lcg(1000 + (i % 8), CB_LEN) + umi = lcg(7000 + i, UMI_LEN) + barcode.append((f"r{i}", cb + umi)) + whitelist.append(cb) + + def write_fq(path, records): + with open(path, "w") as f: + for name, seq in records: + f.write(f"@{name}\n{seq}\n+\n{'I' * len(seq)}\n") + + write_fq(OUT / "cdna.fq", cdna) + write_fq(OUT / "bc.fq", barcode) + (OUT / "whitelist.txt").write_text("\n".join(sorted(set(whitelist))) + "\n") + + star_idx, rustar_idx = OUT / "star_idx", OUT / "rustar_idx" + for idx, exe in ((star_idx, "STAR"), (rustar_idx, RUSTAR)): + idx.mkdir(exist_ok=True) + subprocess.run( + [exe, "--runMode", "genomeGenerate", "--genomeDir", str(idx), + "--genomeFastaFiles", str(OUT / "genome.fa"), + "--genomeSAindexNbases", "7", + "--sjdbGTFfile", str(OUT / "genes.gtf"), "--sjdbOverhang", "89", + "--outFileNamePrefix", str(OUT / f"{idx.name}_")], + check=True, capture_output=True, + ) + + common = [ + "--readFilesIn", str(OUT / "cdna.fq"), str(OUT / "bc.fq"), + "--soloType", "CB_UMI_Simple", + "--soloCBwhitelist", str(OUT / "whitelist.txt"), + "--soloCBstart", "1", "--soloCBlen", str(CB_LEN), + "--soloUMIstart", str(CB_LEN + 1), "--soloUMIlen", str(UMI_LEN), + "--soloFeatures", "Gene", + "--sjdbGTFfile", str(OUT / "genes.gtf"), + *(["--clipAdapterType", "CellRanger4"] if os.environ.get("CR4", "1") == "1" else []), + "--clip5pNbases", "5", + "--clip3pNbases", "3", + "--outSAMtype", "SAM", + ] + subprocess.run(["STAR", "--genomeDir", str(star_idx), *common, + "--outFileNamePrefix", str(OUT / "star_")], + check=True, capture_output=True) + subprocess.run([RUSTAR, "--runMode", "alignReads", "--genomeDir", str(rustar_idx), + *common, "--outFileNamePrefix", str(OUT / "rustar_")], + check=True, capture_output=True) + + def primary(path): + rows = {} + for line in open(path): + if line.startswith("@"): + continue + f = line.split("\t") + flag = int(f[1]) + if flag & 0x900 or flag & 0x4: + continue + rows[f[0]] = (int(f[3]), f[5]) + return rows + + a = primary(OUT / "star_Aligned.out.sam") + b = primary(OUT / "rustar_Aligned.out.sam") + shared = sorted(set(a) & set(b)) + + def lead_clip(cigar): + n = "" + for c in cigar: + if c.isdigit(): + n += c + else: + return int(n) if c == "S" else 0 + return 0 + + deltas = {} + clip_deltas = {} + for r in shared: + d = b[r][0] - a[r][0] + deltas[d] = deltas.get(d, 0) + 1 + cd = lead_clip(b[r][1]) - lead_clip(a[r][1]) + clip_deltas[cd] = clip_deltas.get(cd, 0) + 1 + + print(f"STAR primary: {len(a)} rustar primary: {len(b)} shared: {len(shared)}") + print("POS delta (rustar - STAR):", dict(sorted(deltas.items()))) + print("leading soft-clip delta:", dict(sorted(clip_deltas.items()))) + agree = deltas.get(0, 0) + print(f"identical POS: {agree}/{len(shared)}") + if len(shared) and agree == len(shared): + print("CR4 CLIP DIFF: NONE") + else: + for r in shared[:5]: + if b[r][0] != a[r][0]: + print(f" {r}: STAR {a[r]} rustar {b[r]}") + + +main()