diff --git a/CHANGELOG.md b/CHANGELOG.md index 857f467..f071515 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- `Unmapped.out.mate*` names are written as STAR writes them: `@q9 0:N: ` for single-end reads, and `@p 0:N: 10` / `@p 1:N: 10` for pairs, with the filter flag read from the FASTQ name line and whether each mate mapped. A FASTQ mate 2 record carries mate 1's name, and a pair's records are no longer clipped. - FASTA reads are accepted, in one or several lines, mixed with FASTQ, as in STAR. Their records have `*` for quality, and they go to `Unmapped.out.mate*` as FASTA. - A blank line at the end of a read file no longer stops the run. A blank line with reads after it is an error. See `DIVERGENCE.md` §4.8. - Mates named `SRR1.5.1` / `SRR1.5.2` are accepted as a pair. `SRR1.1` / `SRR1.2` are still refused, since there the number is the read's. diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 38edd78..4c68f74 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -50,6 +50,27 @@ impl UnmappedFastqWriter { } } +/// The name of an `Unmapped.out.mate*` record, without its `@` or `>`: the read +/// name, the 0-based mate, the filter flag, an empty barcode field, then `status`. +/// That is `""` for single-end reads, and for pairs a space and whether each mate +/// mapped (` 10`). Layout observed in STAR 2.7.11b's output, source not read. +pub fn unmapped_name(name: &str, mate: u8, filter_flag: u8, status: &str) -> String { + format!("{name} {mate}:{}: {status}", filter_flag as char) +} + +/// `Y` when the first word of a FASTQ description has `:Y:` from its second byte. +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 { @@ -63,6 +84,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 @@ -303,6 +326,7 @@ impl FastqReader { sequence, quality, file_index: self.file_index, + filter_flag: filter_flag(record.description()), })) } Some(Err(e)) => Err(Error::io(e, &self.path)), @@ -366,6 +390,7 @@ impl FastqReader { sequence, quality: Vec::new(), file_index: self.file_index, + filter_flag: b'N', }) } @@ -895,6 +920,26 @@ mod tests { FastqReader::open(f.path(), None).unwrap() } + /// 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 fasta_reads_have_empty_quality_and_may_span_lines() { let mut r = reader_over(">f1 desc\nAC\nGT\n>f2/1\nTT\n@q3\nAA\n+\nII\n>f4\nCC"); diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..6389626 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2016,7 +2016,12 @@ fn align_reads_single_end( } let unmapped_m1 = if write_unmapped_fastq { vec![( - out_read_name.clone(), + crate::io::fastq::unmapped_name( + &out_read_name, + 0, + read.filter_flag, + "", + ), read.sequence.clone(), read.quality.clone(), )] @@ -2232,7 +2237,12 @@ fn align_reads_single_end( }; let unmapped_m1 = if write_unmapped_fastq && is_unmapped_se { vec![( - out_read_name.clone(), + crate::io::fastq::unmapped_name( + &out_read_name, + 0, + read.filter_flag, + "", + ), read.sequence.clone(), read.quality.clone(), )] @@ -3583,12 +3593,18 @@ fn align_reads_paired_end( } else { paired_read.name.clone() }; - let (fastx_name1, fastx_name2) = if out_read_id_number { - (out_read_name.clone(), out_read_name.clone()) - } else { + // STAR's mate 2 record carries mate 1's name and flag, except that a + // FASTA mate 2 keeps its own name. + let fastx_names = |status: &str| { + let name2 = if !out_read_id_number && paired_read.mate2.quality.is_empty() { + &paired_read.mate2.name + } else { + &out_read_name + }; + let flag = paired_read.mate1.filter_flag; ( - paired_read.mate1.name.clone(), - paired_read.mate2.name.clone(), + crate::io::fastq::unmapped_name(&out_read_name, 0, flag, status), + crate::io::fastq::unmapped_name(name2, 1, flag, status), ) }; @@ -3643,14 +3659,15 @@ fn align_reads_paired_end( } } let (um1, um2) = if write_unmapped_fastq { + let (fastx_name1, fastx_name2) = fastx_names(" 00"); ( vec![( - fastx_name1.clone(), + fastx_name1, paired_read.mate1.sequence.clone(), paired_read.mate1.quality.clone(), )], vec![( - fastx_name2.clone(), + fastx_name2, paired_read.mate2.sequence.clone(), paired_read.mate2.quality.clone(), )], @@ -3971,9 +3988,26 @@ fn align_reads_paired_end( let (unmapped_mate1, unmapped_mate2) = if write_unmapped_fastq { let pair_unmapped = results.is_empty() || has_half_mapped; if pair_unmapped { + let status = match results.first() { + Some(PairedAlignmentResult::HalfMapped { + mate1_is_mapped: true, + .. + }) => " 10", + Some(PairedAlignmentResult::HalfMapped { .. }) => " 01", + _ => " 00", + }; + let (fastx_name1, fastx_name2) = fastx_names(status); ( - vec![(fastx_name1.clone(), m1_seq.clone(), m1_qual.clone())], - vec![(fastx_name2.clone(), m2_seq.clone(), m2_qual.clone())], + vec![( + fastx_name1, + paired_read.mate1.sequence.clone(), + paired_read.mate1.quality.clone(), + )], + vec![( + fastx_name2, + paired_read.mate2.sequence.clone(), + paired_read.mate2.quality.clone(), + )], ) } else { (Vec::new(), Vec::new()) 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..b269f17 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -964,6 +964,116 @@ fn test_unmapped_reads_output() { // Test 8 — Two-pass mode // --------------------------------------------------------------------------- +/// `Unmapped.out.mate*` byte for byte. Expected text is STAR 2.7.11b's output on +/// this genome and these reads: the name fields, mate 2 carrying mate 1's name and +/// flag in FASTQ but its own name in FASTA, and the mates not clipped. +#[test] +fn test_unmapped_mate_files_match_star() { + let tmpdir = TempDir::new().unwrap(); + let genome = lcg_seq(88888, 20000); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("idx"); + build_index(&fasta, &genome_dir, "7", None); + let seq = |bytes: &[u8]| String::from_utf8(bytes.to_vec()).unwrap(); + let junk = seq(&lcg_seq(424_242, 100)); + let junk2 = seq(&lcg_seq(777_777, 100)); + let mapped = seq(&genome[2000..2100]); + let all_n = "N".repeat(100); + let qual = "I".repeat(100); + let fq = |name: &str, s: &str| format!("@{name}\n{s}\n+\n{qual}\n"); + let fa = |name: &str, s: &str| format!(">{name}\n{s}\n"); + + let run = |tag: &str, files: &[(&str, String)], extra: &[&str]| -> Vec { + let mut cmd = cargo_bin_cmd!("rustar-aligner"); + cmd.arg("--genomeDir").arg(&genome_dir); + cmd.arg("--readFilesIn"); + for (file, content) in files { + let path = tmpdir.path().join(format!("{tag}_{file}")); + fs::write(&path, content).unwrap(); + cmd.arg(path); + } + let prefix = tmpdir.path().join(format!("{tag}_")); + cmd.arg("--outFileNamePrefix").arg(&prefix); + cmd.args(["--outSAMtype", "None", "--outReadsUnmapped", "Fastx"]); + cmd.args(extra); + cmd.assert().success(); + (1..=files.len()) + .map(|m| { + let name = format!("{tag}_Unmapped.out.mate{m}"); + fs::read_to_string(tmpdir.path().join(name)).unwrap() + }) + .collect() + }; + + // Single-end FASTQ: the filter flag comes from the second name field. + let reads = [ + fq("q9", &junk), + fq("m", &mapped), + fq("q10 1:Y:0:ACGT", &junk), + fq("q11/1 2:N:0", &junk), + ] + .concat(); + assert_eq!( + run("se", &[("r.fq", reads)], &[]), + [[ + fq("q9 0:N: ", &junk), + fq("q10 0:Y: ", &junk), + fq("q11 0:N: ", &junk) + ] + .concat()] + ); + + // Single-end FASTA: never flagged. + let reads = [fa("f1 1:Y:0", &junk), fa("f2/1", &junk)].concat(); + assert_eq!( + run("sefa", &[("r.fa", reads)], &[]), + [[fa("f1 0:N: ", &junk), fa("f2 0:N: ", &junk)].concat()] + ); + + // Pair, FASTQ: mate 1's name and flag on both records. + let m1 = fq("p_1 1:Y:0", &junk); + let m2 = fq("p_2 2:N:0", &junk2); + assert_eq!( + run("pe", &[("1.fq", m1), ("2.fq", m2)], &[]), + [fq("p_1 0:Y: 00", &junk), fq("p_1 1:Y: 00", &junk2)] + ); + + // Pair, FASTA: each mate's own name. + let m1 = fa("p_1", &junk); + let m2 = fa("p_2", &junk2); + assert_eq!( + run("pefa", &[("1.fa", m1), ("2.fa", m2)], &[]), + [fa("p_1 0:N: 00", &junk), fa("p_2 1:N: 00", &junk2)] + ); + + // Half-mapped pairs: which mate mapped, and the records are not clipped. + let m1 = [fq("h1", &mapped), fq("h2", &all_n)].concat(); + let m2 = [fq("h1", &all_n), fq("h2", &mapped)].concat(); + let half = run( + "half", + &[("1.fq", m1), ("2.fq", m2)], + &[ + "--outFilterScoreMinOverLread", + "0.3", + "--outFilterMatchNminOverLread", + "0.3", + "--clip5pNbases", + "3", + "3", + "--clip3pNbases", + "2", + "2", + ], + ); + assert_eq!( + half, + [ + [fq("h1 0:N: 10", &mapped), fq("h2 0:N: 01", &all_n)].concat(), + [fq("h1 1:N: 10", &all_n), fq("h2 1:N: 01", &mapped)].concat(), + ] + ); +} + #[test] fn test_two_pass_mode() { let tmpdir = TempDir::new().unwrap();