diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..f5129e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **CellRanger4 3' poly-A trim** (`clip::cellranger4::poly_tail_3p`) and the + `--clipAdapterType CellRanger4` clip plumbing, plus `--clip5pAdapterSeq` and + `--clip5pAdapterMMp`. The 5' TSO trim is an overlap alignment and waits on a + dedicated deterministic-SIMD crate, so `CellRanger4` is still rejected at + parse time rather than half-applied. + - **STARsolo single-cell quantification (`--soloType`)** — the 10x Chromium / plate-based count-matrix pipeline, ported from STAR and verified against real STARsolo (#90). diff --git a/src/clip/cellranger4.rs b/src/clip/cellranger4.rs new file mode 100644 index 0000000..e0d7a45 --- /dev/null +++ b/src/clip/cellranger4.rs @@ -0,0 +1,114 @@ +//! `--clipAdapterType CellRanger4`: the 10x Chromium v4 clipping rules. +//! +//! Two independent trims, ported from STAR's `ClipCR4.cpp` and +//! `ClipMate_clipChunk.cpp`: +//! +//! - a 3' poly-A tail trim, scored base by base from the 3' end; +//! - a 5' template-switch-oligo trim, which is an overlap alignment of the TSO +//! against the first 91 bases of the read. +//! +//! Only the poly-A trim is implemented here. The 5' TSO trim is an overlap +//! alignment, for which STAR links the Opal SIMD library; rustar will take +//! that from a dedicated deterministic-SIMD crate rather than carrying a +//! second aligner in-tree. Until then `--clipAdapterType CellRanger4` is +//! rejected rather than silently doing half the job. + +/// Number of 3' bases to trim as a CellRanger4 poly-A tail. +/// +/// STAR `ClipCR4::polyTail3p`. Walks in from the 3' end scoring `+1` per `A` +/// and `-2` per non-`A`, and remembers the longest prefix of that walk whose +/// running score still clears a 70% density threshold (`score * 10 >= ib * 7`). +/// It gives up once the score has fallen more than 27 behind the position, and +/// returns nothing unless the remembered score reached 20. +/// +/// `seq` is numeric base codes, so `A == 0`. +pub fn poly_tail_3p(seq: &[u8]) -> usize { + let seq_len = seq.len(); + if seq_len < 20 { + return 0; + } + let mut best_len: i64 = seq_len as i64 - 1; + let mut score: i64 = 0; + let mut best_score: i64 = 0; + for ib in 1..=seq_len as i64 { + if seq[seq_len - ib as usize] == 0 { + score += 1; + if score * 10 >= ib * 7 { + best_len = ib; + best_score = score; + } + } else { + score -= 2; + if ib - score > 27 { + break; + } + } + } + if best_score < 20 { + 0 + } else { + best_len as usize + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// ACGT text to base codes. + fn code(s: &str) -> Vec { + s.bytes() + .map(|b| match b { + b'A' => 0, + b'C' => 1, + b'G' => 2, + b'T' => 3, + _ => 4, + }) + .collect() + } + + #[test] + fn cr4_polya_trim_matches_star() { + // A clean 30-base poly-A tail is trimmed whole. The prefix is + // deliberately A-free: the scan does not stop at the tail boundary, so + // an `A` just upstream of it would legitimately extend the trim. + let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30))); + assert_eq!(poly_tail_3p(&read), 30); + + // No tail: nothing to trim. + let read = code("ACGTACGTACGTACGTACGTACGTACGTACGT"); + assert_eq!(poly_tail_3p(&read), 0); + + // A tail shorter than the score-20 floor is not trimmed, however clean. + let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(10))); + assert_eq!(poly_tail_3p(&read), 0); + } + + #[test] + fn poly_tail_needs_twenty_bases_of_read() { + assert_eq!(poly_tail_3p(&code("AAAAAAAAAAAAAAAAAAA")), 0); // 19 + } + + #[test] + fn poly_tail_scan_does_not_stop_at_the_tail_boundary() { + // STAR keeps scoring past the run of A's, so A-rich sequence just + // upstream extends the trim. Worth pinning: it looks like an off-by-one + // otherwise. + let read = code(&format!("ACGTACGTACGTACGTACGT{}", "A".repeat(30))); + assert_eq!(poly_tail_3p(&read), 34); + } + + #[test] + fn poly_tail_tolerates_a_single_interruption() { + // 70% density is the threshold, so one non-A inside a long tail is + // survivable. + let tail = format!("{}C{}", "A".repeat(15), "A".repeat(15)); + let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{tail}")); + assert!( + poly_tail_3p(&read) >= 15, + "one mismatch should not abandon the tail, got {}", + poly_tail_3p(&read) + ); + } +} diff --git a/src/clip/mod.rs b/src/clip/mod.rs index 20901ec..df64b6e 100644 --- a/src/clip/mod.rs +++ b/src/clip/mod.rs @@ -17,9 +17,12 @@ //! `clip5pAfterAdapterNbases`), then 3' on the 5'-clipped read (`clip3pNbases`, //! then the 3' adapter Hamming scan, then `clip3pAfterAdapterNbases`). //! -//! Only `--clipAdapterType Hamming` (STAR's default) is supported; a 5' Hamming -//! adapter is not a thing STAR itself supports either (only `CellRanger4` mode -//! clips a 5' adapter, the 10x TSO) — that mode is out of scope here. +//! `--clipAdapterType Hamming` (STAR's default) uses the 3' Hamming scan above. +//! `CellRanger4` instead trims a 3' poly-A tail and a 5' TSO; see +//! [`cellranger4`]. A 5' *Hamming* adapter is not a thing STAR supports either: +//! the 5' end only ever carries an adapter under CellRanger4. + +pub mod cellranger4; use crate::io::fastq::encode_base; use crate::params::Parameters; @@ -48,6 +51,12 @@ pub struct ClipParams { pub five: ClipEnd, /// 3' end (STAR `ClipMate` type 1). pub three: ClipEnd, + /// `--clipAdapterType CellRanger4`: replace the Hamming rules with the 10x + /// poly-A / TSO trims. See [`cellranger4`]. + pub cellranger4: bool, + /// `--clip5pAdapterSeq` as base codes: the 10x TSO, clipped from the 5' end + /// under CellRanger4. Empty when none is configured. + pub five_adapter: Vec, } /// Build [`ClipParams`] for `mate` (0 or 1) from the run's `--clip{5,3}pNbases` @@ -56,6 +65,15 @@ pub struct ClipParams { /// per-mate (`clip5p(mate)`/`clip3p(mate)`); the adapter / mmp / after-adapter /// clips are single-valued and apply to both mates. Cheap to build per batch. pub fn clip_params_from(params: &Parameters, mate: usize) -> ClipParams { + // The 5' TSO trim of CellRanger4 needs an overlap alignment, which waits on + // a dedicated deterministic-SIMD crate. Say so once, loudly, rather than + // leaving the user to infer from the output that half the mode ran. + if params.clip_adapter_type == "CellRanger4" && params.clip5p_adapter_seq != "-" && mate == 0 { + log::warn!( + "--clipAdapterType CellRanger4: the 3' poly-A trim is applied, but the 5' \ + adapter (TSO) trim is not yet implemented and --clip5pAdapterSeq is ignored" + ); + } // Encode the adapter to base codes (A=0..T=3) ONCE here. The read reaching // clip_mate is already numeric (io::fastq encodes at read time), so the 3' // Hamming scan compares numeric-vs-numeric — re-encoding the read there turned @@ -65,7 +83,14 @@ pub fn clip_params_from(params: &Parameters, mate: usize) -> ClipParams { } else { params.clip3p_adapter_seq.bytes().map(encode_base).collect() }; + let five_adapter = if params.clip5p_adapter_seq == "-" { + Vec::new() + } else { + params.clip5p_adapter_seq.bytes().map(encode_base).collect() + }; ClipParams { + cellranger4: params.clip_adapter_type == "CellRanger4", + five_adapter, five: ClipEnd { n: params.clip5p(mate), adapter: Vec::new(), @@ -120,6 +145,10 @@ fn local_search(x: &[u8], y: &[u8], p_mm: f64) -> usize { pub fn clip_mate(read: &[u8], p: &ClipParams) -> (usize, usize) { let len = read.len(); + if p.cellranger4 { + return clip_mate_cellranger4(read, p); + } + // ---- 5' end (STAR ClipMate type 0) ---- let five_active = p.five.n > 0; let mut c5 = 0; @@ -155,6 +184,84 @@ pub fn clip_mate(read: &[u8], p: &ClipParams) -> (usize, usize) { (c5, c3) } +/// `--clipAdapterType CellRanger4`, which replaces the Hamming rules entirely +/// (STAR `ClipCR4`): a 5' TSO trim and a 3' poly-A trim. +/// +/// The fixed `--clip{5,3}pNbases` still apply first, as in the Hamming path. +/// The 5' TSO is only configured for the first mate, matching STAR, so mate 2 +/// simply has no 5' adapter and the 5' trim reduces to the fixed clip. +fn clip_mate_cellranger4(read: &[u8], p: &ClipParams) -> (usize, usize) { + let len = read.len(); + + // 5': fixed clip only. The TSO trim is an overlap alignment and waits on a + // dedicated deterministic-SIMD crate. A configured TSO is therefore not + // applied, and the caller is warned once rather than left to infer it from + // the output. + let mut c5 = p.five.n.min(len); + if p.five.n_after > 0 && c5 < len { + c5 += p.five.n_after.min(len - c5); + } + + // 3': fixed clip, then the poly-A scan on the 5'-clipped read. + let s = &read[c5..]; + let sl = s.len(); + let mut c3 = p.three.n.min(sl); + let remaining = sl - c3; + if remaining > 0 { + c3 += cellranger4::poly_tail_3p(&s[..remaining]).min(remaining); + } + if p.three.n_after > 0 && c3 < sl { + c3 += p.three.n_after.min(sl - c3); + } + + (c5, c3) +} + +#[cfg(test)] +mod cr4_wiring_tests { + use super::*; + + fn code(s: &str) -> Vec { + s.bytes().map(encode_base).collect() + } + + const TSO: &str = "AAGCAGTGGTATCAACGCAGAGTACATGGG"; + + fn cr4_params(tso: &str) -> ClipParams { + ClipParams { + cellranger4: true, + five_adapter: code(tso), + five: ClipEnd::default(), + three: ClipEnd::default(), + } + } + + #[test] + fn cellranger4_leaves_a_read_without_either_feature_alone() { + let read = code("CGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGTCGT"); + assert_eq!(clip_mate(&read, &cr4_params(TSO)), (0, 0)); + } + + #[test] + fn cellranger4_without_a_tso_still_trims_polya() { + // Mate 2 has no 5' adapter under STAR, so only the 3' trim applies. + let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30))); + let (c5, c3) = clip_mate(&read, &cr4_params("-")); + assert_eq!(c5, 0); + assert_eq!(c3, 30); + } + + #[test] + fn cellranger4_applies_the_fixed_clip_before_the_polya_trim() { + let read = code(&format!("CGTCGTCGTCGTCGTCGTCG{}", "A".repeat(30))); + let mut p = cr4_params("-"); + p.five.n = 5; + let (c5, c3) = clip_mate(&read, &p); + assert_eq!(c5, 5, "the fixed 5' clip still applies"); + assert_eq!(c3, 30, "and the poly-A tail is trimmed from what remains"); + } +} + #[cfg(test)] mod tests { use super::*; @@ -183,6 +290,8 @@ mod tests { #[test] fn fixed_5p_3p() { let p = ClipParams { + cellranger4: false, + five_adapter: Vec::new(), five: ClipEnd { n: 3, ..Default::default() @@ -217,6 +326,8 @@ mod tests { fn after_adapter_alone_is_noop() { // STAR's inactive-end short-circuit: n_after with no fixed clip and no adapter clips nothing. let p = ClipParams { + cellranger4: false, + five_adapter: Vec::new(), five: ClipEnd { n_after: 4, ..Default::default() @@ -240,6 +351,8 @@ mod tests { #[test] fn no_adapter_configured_only_fixed_clips() { let p = ClipParams { + cellranger4: false, + five_adapter: Vec::new(), five: ClipEnd { n: 2, ..Default::default() diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..7f5cf4a 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -555,6 +555,16 @@ pub struct Parameters { #[arg(long = "clipAdapterType", default_value = "Hamming")] pub clip_adapter_type: String, + /// 5' adapter sequence to clip, one per mate. Only used by + /// `--clipAdapterType CellRanger4`, where it is the 10x template switch + /// oligo. `-` (the default) means none. + #[arg(long = "clip5pAdapterSeq", default_value = "-")] + pub clip5p_adapter_seq: String, + + /// Max mismatch fraction for the 5' adapter. + #[arg(long = "clip5pAdapterMMp", default_value_t = 0.1)] + pub clip5p_adapter_mmp: f64, + /// 3' adapter sequence to clip (Hamming scan), `-` = none #[arg(long = "clip3pAdapterSeq", default_value = "-")] pub clip3p_adapter_seq: String, diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..ca9dc53 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2024,3 +2024,86 @@ fn test_wasp_samtag() { "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)" ); } + +// --------------------------------------------------------------------------- +// --clipAdapterType CellRanger4 +// --------------------------------------------------------------------------- + +/// End to end: a read carrying a poly-A tail aligns over its non-A prefix only. +/// +/// Before this, `--clipAdapterType CellRanger4` parsed and then did nothing, so +/// the tail was carried into the alignment as soft-clipped or mismatching +/// bases. The check is that the flag changes the CIGAR at all, and changes it +/// the way a trim should. +#[test] +fn test_clip_adapter_type_cellranger4() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + // 60 genomic bases followed by a 40-base poly-A tail the genome does not + // have. The tail is long enough to clear the trim's score-20 floor. + let start = 3000usize; + let mut read = genome[start..start + 60].to_vec(); + read.extend(std::iter::repeat_n(b'A', 40)); + + let fastq_path = tmpdir.path().join("reads.fq"); + { + let mut f = fs::File::create(&fastq_path).unwrap(); + writeln!(f, "@polya").unwrap(); + f.write_all(&read).unwrap(); + writeln!(f).unwrap(); + writeln!(f, "+").unwrap(); + writeln!(f, "{}", "I".repeat(read.len())).unwrap(); + } + + let run = |extra: &[&str], out_name: &str| -> Vec { + let output_dir = tmpdir.path().join(out_name); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + let mut args = vec![ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq_path.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]; + args.extend_from_slice(extra); + cargo_bin_cmd!("rustar-aligner") + .args(&args) + .assert() + .success(); + fs::read_to_string(output_dir.join("Aligned.out.sam")) + .unwrap() + .lines() + .filter(|l| !l.starts_with('@')) + .map(ToString::to_string) + .collect() + }; + + let plain = run(&[], "out_plain"); + let clipped = run(&["--clipAdapterType", "CellRanger4"], "out_cr4"); + + // Untrimmed, 40 of 100 bases mismatch and the read fails the mismatch + // filters outright. That is the cost of the flag having been inert. + assert!( + plain.is_empty(), + "without trimming the poly-A read should not pass the mismatch filters, got {plain:?}" + ); + + assert_eq!(clipped.len(), 1, "trimmed, the read aligns"); + let cigar = clipped[0].split('\t').nth(5).unwrap(); + assert!( + cigar.ends_with("40S"), + "the 40-base poly-A tail should be clipped, got {cigar}" + ); + assert!( + cigar.starts_with("60M"), + "the genomic prefix should still align in full, got {cigar}" + ); +}