diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..cc8cd2a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -67,6 +67,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- `--outFilterMismatchNoverReadLmax` is accepted (default 1.0), as in STAR + 2.7.11b. It caps the mismatches of a read at this fraction of its length + (both mates, after clipping), together with `--outFilterMismatchNmax`, in + the same places STAR does: while stitching and extending, in the mapped + filter and for `--quantMode TranscriptomeSAM`. On 4 000 reads and 4 000 + pairs with 0 to 9 mismatches each, at 0.04, the SAM records are identical + to STAR's. + - **`--clipAdapterType CellRanger4` now matches STAR exactly**, both halves of the clip (#198). The 5' TSO trim is an overlap alignment against the first 91 bases of the read, replicating `ClipCR4`'s Opal call (ACGTN alphabet; match +1, mismatch -2, @@ -232,6 +240,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `NH` counts every locus of a multimapper, also when `--outSAMmultNmax` + writes fewer. A read with 3 loci under `--outSAMmultNmax 2` had `NH:i:2`, + STAR writes `NH:i:3`. Single-end and paired-end. +- A read whose FASTQ name has the Illumina filter flag (`@r1 1:Y:0:ACGT`) gets + SAM flag 0x200 (not passing quality control) on all its records, mapped or + unmapped, and on both mates (from mate 1's name), in SAM, BAM and + `Aligned.toTranscriptome.out.bam`, as STAR does. It was never set. + - **STARsolo's stats files as STAR writes them.** `Solo.out/Barcodes.stats`, `/Features.stats` and `/UMIperCellSorted.txt` were not written, SJ and Velocyto had no `Summary.csv`, and the Gene/GeneFull diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index ed87d8c..3fc1f9b 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -65,6 +65,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--outFilterMultimapScoreRange` | `1` | Score range for keeping multi-mappers within best score. | | `--outFilterMismatchNmax` | `10` | Max mismatches per pair. | | `--outFilterMismatchNoverLmax` | `0.3` | Max ratio of mismatches to mapped length. | +| `--outFilterMismatchNoverReadLmax` | `1.0` | Max ratio of mismatches to read length (both mates). | | `--outFilterScoreMin` | `0` | Min absolute alignment score. | | `--outFilterScoreMinOverLread` | `0.66` | Min alignment score normalized to read length. | | `--outFilterMatchNmin` | `0` | Min absolute matched bases. | diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 8fc7feb..8a6b819 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -409,7 +409,8 @@ fn align_read_inner( // does NOT filter clusters by seed coverage. Removed to match STAR. // Step 3: Stitch seeds within each cluster - let scorer = AlignmentScorer::from_params(params); + let mut scorer = AlignmentScorer::from_params(params); + scorer.n_mm_max = params.mismatch_nmax_total(read_seq.len()); let mut transcripts = Vec::new(); // STAR's `trAll`: every window's transcripts, before any cross-window // dedup or filtering, which is what its chimeric detection reads. @@ -577,7 +578,7 @@ fn align_read_inner( || n_match < (params.out_filter_match_nmin_over_lread * lread_m1) as u32 { Some(UnmappedReason::TooShort) - } else if best.n_mismatch > params.out_filter_mismatch_nmax + } else if best.n_mismatch > scorer.n_mm_max || f64::from(best.n_mismatch) / f64::from(r_length.max(1)) > params.out_filter_mismatch_nover_lmax { @@ -709,7 +710,8 @@ pub fn align_paired_read( let len2 = mate2_seq.len(); let debug_pe = !params.read_name_filter.is_empty() && read_name == params.read_name_filter; - let scorer = AlignmentScorer::from_params(params); + let mut scorer = AlignmentScorer::from_params(params); + scorer.n_mm_max = params.mismatch_nmax_total(len1 + len2); let junction_db = if index.junction_db.is_empty() { None } else { @@ -1081,6 +1083,9 @@ pub fn align_paired_read( merge_params.out_filter_match_nmin = 0; merge_params.out_filter_match_nmin_over_lread = 0.0; merge_params.out_filter_score_min_over_lread = 0.0; + // STAR's copyRead hands the merged read the pair's mismatch cap. + merge_params.out_filter_mismatch_nmax = scorer.n_mm_max; + merge_params.out_filter_mismatch_nover_read_lmax = 1.0; let (merged_transcripts, _merged_chim, _n_mapq, _unmapped, _best) = align_read_inner(&merge.merged, read_name, index, &merge_params, false)?; let mut converted: Vec = Vec::new(); @@ -1559,7 +1564,7 @@ fn filter_paired_transcripts( return Some(UnmappedReason::TooShort); } - if combined_nm > params.out_filter_mismatch_nmax + if combined_nm > params.mismatch_nmax_total(mate1_seq.len() + mate2_seq.len()) || f64::from(combined_nm) / f64::from(combined_r_length.max(1)) > params.out_filter_mismatch_nover_lmax { diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 48f9673..579e929 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -45,6 +45,20 @@ impl UnmappedFastqWriter { } } +/// `Y` when the first word of a FASTQ description has `:Y:` from its second byte, +/// STAR's Illumina filter flag (`ReadAlignChunk_processChunks.cpp:121-127`). +fn filter_flag(description: &[u8]) -> u8 { + let field = description + .split(u8::is_ascii_whitespace) + .next() + .unwrap_or(&[]); + if field.get(1..4) == Some(b":Y:") { + b'Y' + } else { + b'N' + } +} + /// A read from a FASTQ file with encoded bases #[derive(Debug, Clone)] pub struct EncodedRead { @@ -58,6 +72,8 @@ pub struct EncodedRead { /// Index of the `--readFilesIn` file this read came from (STAR /// `readFilesIndex`); selects the read group. pub file_index: u32, + /// `Y` when the second field of a FASTQ name line reads `n:Y:...`, else `N`. + pub filter_flag: u8, } /// A paired-end read from two FASTQ files @@ -235,6 +251,7 @@ impl FastqReader { sequence, quality, file_index: self.file_index, + filter_flag: filter_flag(record.description()), })) } Some(Err(e)) => Err(Error::from(e)), @@ -496,6 +513,26 @@ mod tests { assert_eq!(encode_base(b'S'), 4); } + /// Expected flags are STAR's on the same name lines. + #[test] + fn filter_flag_needs_y_as_the_second_colon_field_of_a_one_byte_first_field() { + for (description, flag) in [ + (&b"1:Y:0:ACGT"[..], b'Y'), + (b"2:Y:18:GG", b'Y'), + (b"1:Y:0 2:N:0", b'Y'), + (b"1:N:0:ACGT", b'N'), + (b"Y", b'N'), + (b":Y:", b'N'), + (b"1:Y", b'N'), + (b"ab:Y:", b'N'), + (b"1:y:0", b'N'), + (b"x\ty:Y:z", b'N'), + (b"", b'N'), + ] { + assert_eq!(filter_flag(description), flag, "{description:?}"); + } + } + #[test] fn test_decode_base() { assert_eq!(decode_base(0), b'A'); diff --git a/src/io/sam.rs b/src/io/sam.rs index 40a03e5..9cd64d2 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -213,7 +213,7 @@ impl SamWriter { read_qual, genome, mapq, - max_output, // NH = number of reported alignments + n_alignments, // NH counts every locus, not only the written ones (STAR nTr) hit_index + 1, // 1-based params.out_sam_attr_ih_start, attrs, @@ -363,7 +363,7 @@ impl SamWriter { aligned_qual, genome, mapq, - max_output, // NH = number of reported alignments + n_alignments, // NH counts every locus, not only the written ones (STAR nTr) hit_index + 1, // 1-based params.out_sam_attr_ih_start, attrs, @@ -482,7 +482,7 @@ impl SamWriter { paired_aln.insert_size, params, ), - max_output, // NH = number of reported alignments + n_alignments, // NH counts every locus, not only the written ones (STAR nTr) hit_index, params.out_sam_attr_ih_start, combined_score, @@ -513,7 +513,7 @@ impl SamWriter { -paired_aln.insert_size, // Negative for mate2 (mode 1 default) params, ), - max_output, // NH = number of reported alignments + n_alignments, // NH counts every locus, not only the written ones (STAR nTr) hit_index, params.out_sam_attr_ih_start, combined_score, @@ -3861,13 +3861,13 @@ mod tests { // Only 3 records output despite 5 transcripts assert_eq!(records.len(), 3); - // NH should be 3 (number of reported alignments) + // NH counts all 5 loci, as STAR's does (`outputTranscriptSAM(.., nTrOutSAM, ..)`) for rec in &records { let data = rec.data(); assert_eq!( data.get(&Tag::ALIGNMENT_HIT_COUNT), - Some(&Value::from(3_i32)), - "NH should be 3 (capped by outSAMmultNmax)" + Some(&Value::from(5_i32)), + "NH is not capped by outSAMmultNmax" ); } diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..138c871 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -833,6 +833,51 @@ fn apply_rg_per_file( } } +/// Flag every record of a read whose FASTQ name has the Illumina filter flag +/// `n:Y:` as not passing quality control (0x200), mapped or not, as STAR does +/// (`ReadAlign_outputTranscriptSAM.cpp:27,91`, `ReadAlign_alignBAM.cpp:146,198`). +/// STAR applies `--outSAMflagAND/OR` after setting the bit, on mapped records +/// only, so they are applied again here. +fn apply_qc_fail( + records: &mut [noodles::sam::alignment::record_buf::RecordBuf], + filter_flag: u8, + params: &Parameters, + and_or: bool, +) { + if filter_flag != b'Y' { + return; + } + for rec in records { + rec.flags_mut() + .insert(noodles::sam::alignment::record::Flags::QC_FAIL); + if and_or && !rec.flags().is_unmapped() { + crate::io::sam::apply_sam_flag_or_and(rec, params); + } + } +} + +/// [`apply_rg_per_file`] and [`apply_qc_fail`] over a batch's genomic and +/// transcriptome records; `read(i)` is the batch's `i`-th read (mate 1 for pairs). +fn apply_per_read_batch( + params: &Parameters, + out: &mut BatchOut, + read: impl Fn(usize) -> (u32, u8), +) { + for (i, res) in out.iter_mut().enumerate() { + if let Ok(o) = res { + let filter_flag = read(i).1; + apply_qc_fail(&mut o.sam_records.records, filter_flag, params, true); + apply_qc_fail( + &mut o.transcriptome_records.records, + filter_flag, + params, + false, + ); + } + } + apply_rg_per_file_batch(params, out, |i| read(i).0); +} + /// [`apply_rg_per_file`] over a batch's genomic and transcriptome records; /// `file_index(i)` is the file of the batch's `i`-th read. fn apply_rg_per_file_batch( @@ -2252,7 +2297,9 @@ fn align_reads_single_end( }) }) .collect(); - apply_rg_per_file_batch(params, &mut out, |i| batch[i].file_index); + apply_per_read_batch(params, &mut out, |i| { + (batch[i].file_index, batch[i].filter_flag) + }); out }; Arc::new(process) @@ -2876,12 +2923,11 @@ fn align_reads_solo( }) }) .collect(); - if !rg_ids.is_empty() { - for (i, res) in out.iter_mut().enumerate() { - if let Ok(o) = res { - let fi = batch[i].cdna.file_index; - apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi); - } + for (i, res) in out.iter_mut().enumerate() { + if let Ok(o) = res { + let r = &batch[i].cdna; + apply_qc_fail(&mut o.sam_records.records, r.filter_flag, params, true); + apply_rg_per_file(&mut o.sam_records.records, &rg_ids, r.file_index); } } out @@ -3401,12 +3447,11 @@ fn align_reads_solo_pe( }) }) .collect(); - if !rg_ids.is_empty() { - for (i, res) in out.iter_mut().enumerate() { - if let Ok(o) = res { - let fi = batch[i].mate1.file_index; - apply_rg_per_file(&mut o.sam_records.records, &rg_ids, fi); - } + for (i, res) in out.iter_mut().enumerate() { + if let Ok(o) = res { + let r = &batch[i].mate1; + apply_qc_fail(&mut o.sam_records.records, r.filter_flag, params, true); + apply_rg_per_file(&mut o.sam_records.records, &rg_ids, r.file_index); } } out @@ -3994,7 +4039,9 @@ fn align_reads_paired_end( }) }) .collect(); - apply_rg_per_file_batch(params, &mut out, |i| batch[i].mate1.file_index); + apply_per_read_batch(params, &mut out, |i| { + (batch[i].mate1.file_index, batch[i].mate1.filter_flag) + }); out }; Arc::new(process) diff --git a/src/params/mod.rs b/src/params/mod.rs index 47e4ecf..2e98eb4 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -861,6 +861,10 @@ pub struct Parameters { #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)] pub out_filter_mismatch_nover_lmax: f64, + /// Max ratio of mismatches to read length (both mates) + #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)] + pub out_filter_mismatch_nover_read_lmax: f64, + /// Min alignment score (absolute) #[arg(long = "outFilterScoreMin", default_value_t = 0)] pub out_filter_score_min: i32, @@ -1670,6 +1674,14 @@ impl Parameters { .collect() } + /// STAR's per-read `outFilterMismatchNmaxTotal` (`ReadAlign_oneRead.cpp:78`): + /// `outFilterMismatchNmax` capped by `outFilterMismatchNoverReadLmax` times the + /// read length, both mates after clipping. + pub fn mismatch_nmax_total(&self, read_len: usize) -> u32 { + let by_len = (self.out_filter_mismatch_nover_read_lmax * read_len as f64) as u32; + self.out_filter_mismatch_nmax.min(by_len) + } + /// Read group ID emitted on SAM records (the first RG line's `ID:` value). /// Returns `None` when no RG line is configured. pub fn primary_rg_id(&self) -> Result, crate::error::Error> { diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 0a02059..995b89f 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -1592,9 +1592,7 @@ fn prepare_for_projection( /// STAR's `(nMM + nMM1) > min(outFilterMismatchNmaxTotal, /// outFilterMismatchNoverLmax*(Lread-1))` test; `lread` is STAR's `Lread`. fn within_mismatch_budget(n_mm: u32, read_len_sum: u32, lread: u32, params: &Parameters) -> bool { - // `outFilterMismatchNmaxTotal = min(outFilterMismatchNmax, - // outFilterMismatchNoverReadLmax*(readLength[0]+readLength[1]))`, ratio 1. - let nmax_total = params.out_filter_mismatch_nmax.min(read_len_sum); + let nmax_total = params.mismatch_nmax_total(read_len_sum as usize); let rel = (params.out_filter_mismatch_nover_lmax * f64::from(lread.saturating_sub(1))).floor() as u32; n_mm <= nmax_total.min(rel) diff --git a/src/solo/mod.rs b/src/solo/mod.rs index b91d1fb..19b8ea1 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -2249,6 +2249,7 @@ mod tests { sequence: seq.bytes().map(encode_base).collect(), quality: qual.bytes().collect(), file_index: 0, + filter_flag: b'N', } } @@ -2750,6 +2751,7 @@ mod tests { sequence: seq.bytes().map(encode_base).collect(), quality: vec![b'I'; seq.len()], file_index: 0, + filter_flag: b'N', }; // Default: exactly CB+UMI = 26 bases. let checked = BarcodeReadLength { @@ -2890,6 +2892,7 @@ mod tests { sequence: vec![2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], quality: vec![b'I'; 12], file_index: 0, + filter_flag: b'N', }; let bc = layout.extract(&read).expect("adapter present"); assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); @@ -2901,6 +2904,7 @@ mod tests { sequence: vec![1, 1, 2, 2, 2, 2, 0, 1, 2, 3, 3, 3, 3, 3], quality: vec![b'I'; 14], file_index: 0, + filter_flag: b'N', }; let bc = layout.extract(&read).expect("adapter present, shifted"); assert_eq!(bc.cb_seq, vec![2, 2, 2, 2]); @@ -2912,6 +2916,7 @@ mod tests { sequence: vec![2, 2, 2, 2, 0, 1, 2, 0, 3, 3, 3, 3], quality: vec![b'I'; 12], file_index: 0, + filter_flag: b'N', }; assert!(layout.extract(&read).is_ok()); @@ -2922,6 +2927,7 @@ mod tests { sequence: vec![2; 12], quality: vec![b'I'; 12], file_index: 0, + filter_flag: b'N', }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } @@ -2999,6 +3005,7 @@ mod tests { sequence: vec![2; 12], quality: vec![b'I'; 12], file_index: 0, + filter_flag: b'N', }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoAdapter)); } @@ -3029,6 +3036,7 @@ mod tests { sequence: vec![0, 1, 2, 3, 3, 3, 3], quality: vec![b'I'; 7], file_index: 0, + filter_flag: b'N', }; assert_eq!(layout.extract(&read), Err(BarcodeMiss::NoCb)); } diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 7e71c04..3dde486 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -3673,3 +3673,119 @@ fn solo_out_files(output_dir: &Path) -> Vec { out.sort(); out } + +/// Map `mates` (one FASTQ text per mate) against the genome with a 300 bp +/// segment planted three times, and return `(QNAME, FLAG, NH)` per record. +fn map_on_repeat_genome(mates: &[String], extra: &[&str]) -> Vec<(String, u16, String)> { + let tmpdir = TempDir::new().unwrap(); + let mut genome = build_genome(); + let repeat = lcg_seq(4444, 300); + for start in [2000, 6000, 14000] { + genome[start..start + 300].copy_from_slice(&repeat); + } + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + let mut fastqs = Vec::new(); + for (i, text) in mates.iter().enumerate() { + let path = tmpdir.path().join(format!("r{i}.fq")); + fs::write(&path, text).unwrap(); + fastqs.push(path.to_str().unwrap().to_string()); + } + let prefix = format!("{}/", tmpdir.path().display()); + let mut cmd = cargo_bin_cmd!("rustar-aligner"); + cmd.args(["--genomeDir", genome_dir.to_str().unwrap(), "--readFilesIn"]) + .args(&fastqs) + .args(["--outFileNamePrefix", &prefix, "--outSAMunmapped", "Within"]) + .args(extra) + .assert() + .success(); + let sam = fs::read_to_string(tmpdir.path().join("Aligned.out.sam")).unwrap(); + sam.lines() + .filter(|l| !l.starts_with('@')) + .map(|l| { + let f: Vec<&str> = l.split('\t').collect(); + let nh = f[11..].iter().find(|t| t.starts_with("NH:")).unwrap(); + (f[0].to_string(), f[1].parse().unwrap(), nh.to_string()) + }) + .collect() +} + +fn fastq(records: &[(&str, &[u8])]) -> String { + let mut out = String::new(); + for (name, seq) in records { + let seq = String::from_utf8_lossy(seq); + out.push_str(&format!("@{name}\n{seq}\n+\n{}\n", "I".repeat(seq.len()))); + } + out +} + +/// Three things STAR 2.7.11b does, checked against its output on the same kind +/// of reads: NH counts all loci when `--outSAMmultNmax` writes fewer +/// (`ReadAlign_outputAlignments.cpp:168-177`), a FASTQ name with `1:Y:` sets 0x200 +/// on mapped and unmapped records (`ReadAlign_outputTranscriptSAM.cpp:27,91`), and +/// `--outFilterMismatchNoverReadLmax` caps the mismatches by read length +/// (`ReadAlign_oneRead.cpp:78`). +#[test] +fn test_nh_qc_fail_and_mismatch_over_read_length_match_star() { + let genome = build_genome(); + let repeat = lcg_seq(4444, 300); + let noise = lcg_seq(777, 100); + // 6 mismatches in 100 bases: over 0.05 * 100, under the default 10. + let mut mm6 = genome[16000..16100].to_vec(); + for p in [10, 25, 40, 55, 70, 85] { + mm6[p] = if mm6[p] == b'A' { b'C' } else { b'A' }; + } + let se = fastq(&[ + ("multi 1:Y:0:ACGT", &repeat[0..100]), + ("noise 1:Y:0:ACGT", &noise), + ("unique 1:N:0:ACGT", &genome[1000..1100]), + ("mm6", &mm6), + ]); + let args = [ + "--outSAMmultNmax", + "2", + "--alignEndsType", + "EndToEnd", + "--outFilterMismatchNoverReadLmax", + "0.05", + ]; + let mut got = map_on_repeat_genome(&[se], &args); + got.sort(); + // STAR 2.7.11b's FLAG and NH on these reads. + let want: Vec<(String, u16, String)> = [ + ("mm6", 4, "NH:i:0"), + ("multi", 512, "NH:i:3"), + ("multi", 768, "NH:i:3"), + ("noise", 516, "NH:i:0"), + ("unique", 0, "NH:i:1"), + ] + .iter() + .map(|(n, f, t)| (n.to_string(), *f, t.to_string())) + .collect(); + assert_eq!(got, want); + + // Pairs: NH counts all 6 loci (3 copies of each mate pair up across copies), + // and mate 1's flag marks both mates. STAR 2.7.11b's FLAG and NH. + let m1 = fastq(&[ + ("multi 1:Y:0:A", &repeat[0..100]), + ("noise 1:Y:0:A", &noise), + ]); + let m2 = fastq(&[ + ("multi 2:N:0:A", &rc(&repeat[200..300])), + ("noise 2:N:0:A", &rc(&noise)), + ]); + let got = map_on_repeat_genome(&[m1, m2], &["--outSAMmultNmax", "2"]); + let want: Vec<(String, u16, String)> = [ + ("multi", 611, "NH:i:6"), + ("multi", 659, "NH:i:6"), + ("multi", 867, "NH:i:6"), + ("multi", 915, "NH:i:6"), + ("noise", 589, "NH:i:0"), + ("noise", 653, "NH:i:0"), + ] + .iter() + .map(|(n, f, t)| (n.to_string(), *f, t.to_string())) + .collect(); + assert_eq!(got, want); +} diff --git a/tests/parameter_surface.rs b/tests/parameter_surface.rs index cdd9cc6..50d060e 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -96,7 +96,6 @@ const NOT_YET_ACCEPTED: &[&str] = &[ "alignInsertionFlush", "alignSoftClipAtReferenceEnds", "alignTranscriptsPerReadNmax", - "outFilterMismatchNoverReadLmax", "seedNoneLociPerWindow", // Long reads. "winReadCoverageBasesMin",