Unmapped.out.mate* names and sequences as STAR writes them

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.

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");

src/lib.rs

@@ -2016,7 +2016,12 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                         }
                         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<W: AlignmentWriter + ?Sized>(
                     };
                     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<W: AlignmentWriter + ?Sized>(
                     } 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<W: AlignmentWriter + ?Sized>(
                             }
                         }
                         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<W: AlignmentWriter + ?Sized>(
                     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())

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));
     }

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<String> {
+        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();