Merge feat/transcriptome-sam: add --quantMode TranscriptomeSAM output

# Conflicts:
#	src/io/sam.rs
#	src/params.rs

Author: ewels

src/align/read_align.rs

@@ -17,7 +17,7 @@ use std::hash::{DefaultHasher, Hash, Hasher};
 /// via rayon, so we instead fold the read name into the seed — this keeps tie
 /// breaks reproducible regardless of thread count while still honoring the
 /// user's `--runRNGseed` value.
-fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
+pub(crate) fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
     let mut hasher = DefaultHasher::new();
     read_name.hash(&mut hasher);
     run_rng_seed.wrapping_mul(hasher.finish().wrapping_add(1))

src/io/bam.rs

@@ -2,6 +2,7 @@
 use crate::error::Error;
 use crate::genome::Genome;
 use crate::params::Parameters;
+use crate::quant::transcriptome::TranscriptomeIndex;
 use noodles::bam;
 use noodles::sam;
 use noodles::sam::alignment::io::Write as SamWrite;
@@ -41,26 +42,41 @@ pub struct BamWriter {
 }
 
 impl BamWriter {
-    /// Create a new BAM writer with header from genome index
-    ///
-    /// # Arguments
-    /// * `output_path` - Path to output BAM file
-    /// * `genome` - Genome index with chromosome information
-    /// * `params` - Parameters (for @PG header)
-    pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result<Self, Error> {
-        let file = File::create(output_path)?;
-        let buf_writer = BufWriter::new(file);
-
-        // Reuse SAM header building logic
-        let header = crate::io::sam::build_sam_header(genome, params)?;
-
-        // Create BAM writer (Writer::new handles BGZF compression internally)
+    /// Create a BAM writer at `output_path` with the given prepared header.
+    fn with_header(output_path: &Path, header: sam::Header) -> Result<Self, Error> {
+        let buf_writer = BufWriter::new(File::create(output_path)?);
         let mut writer = bam::io::Writer::new(buf_writer);
         writer.write_header(&header)?;
-
         Ok(Self { writer, header })
     }
 
+    /// Create a new BAM writer with header from genome index.
+    pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result<Self, Error> {
+        Self::with_header(
+            output_path,
+            crate::io::sam::build_sam_header(genome, params)?,
+        )
+    }
+
+    /// Create a BAM writer whose @SQ header lists the transcripts (not
+    /// chromosomes) from `tr_idx`.  Used for
+    /// `Aligned.toTranscriptome.out.bam`.
+    pub fn create_transcriptome(
+        output_path: &Path,
+        tr_idx: &TranscriptomeIndex,
+        params: &Parameters,
+    ) -> Result<Self, Error> {
+        let refs = tr_idx
+            .tr_ids
+            .iter()
+            .zip(tr_idx.tr_length.iter())
+            .map(|(id, len)| (id.as_str(), *len as usize));
+        Self::with_header(
+            output_path,
+            crate::io::sam::build_sam_header_from_refs(refs, params)?,
+        )
+    }
+
     /// Write batch of buffered records (for parallel processing)
     ///
     /// # Arguments
@@ -191,6 +207,41 @@ mod tests {
         assert!(result.is_ok(), "Finishing BAM file should succeed");
     }
 
+    #[test]
+    fn test_bam_transcriptome_writer_creation() {
+        use crate::junction::gtf::GtfRecord;
+        use crate::quant::transcriptome::TranscriptomeIndex;
+        use std::collections::HashMap;
+
+        let genome = create_test_genome();
+        // Stretch genome / exon to fit the tiny test chromosome.
+        let mut attrs = HashMap::new();
+        attrs.insert("gene_id".to_string(), "G1".to_string());
+        attrs.insert("transcript_id".to_string(), "T1".to_string());
+        let exons = vec![GtfRecord {
+            seqname: "chr1".to_string(),
+            feature: "exon".to_string(),
+            start: 1,
+            end: 8,
+            strand: '+',
+            attributes: attrs,
+        }];
+        let tr_idx = TranscriptomeIndex::from_gtf_exons(&exons, &genome).unwrap();
+        assert_eq!(tr_idx.n_transcripts(), 1);
+
+        let params = create_test_params();
+        let temp_file = NamedTempFile::new().unwrap();
+        let writer = BamWriter::create_transcriptome(temp_file.path(), &tr_idx, &params);
+        assert!(
+            writer.is_ok(),
+            "transcriptome BAM writer creation should succeed"
+        );
+
+        // Header should contain exactly 1 @SQ entry (matching n_transcripts).
+        let writer = writer.unwrap();
+        assert_eq!(writer.header.reference_sequences().len(), 1);
+    }
+
     #[test]
     fn test_bam_batch_write() {
         let genome = create_test_genome();

src/io/sam.rs

@@ -547,6 +547,118 @@ impl SamWriter {
         Ok(records)
     }
 
+    /// Build transcriptome-space SAM records for `--quantMode TranscriptomeSAM`.
+    ///
+    /// Each projected `Transcript` is converted to a record where:
+    ///   * `chr_idx` is the transcript index (matches the transcriptome
+    ///     header's @SQ order),
+    ///   * `genome_start` is the 0-based transcript-space position (→ POS =
+    ///     t-space_pos + 1),
+    ///   * splice-aware tags (`jM`, `jI`, `XS`) are not emitted (splices
+    ///     collapse in t-space and have no meaning there),
+    ///   * standard tags (`NH`, `HI`, `AS`, `NM`/`nM`, `MD`) are emitted per
+    ///     the `--outSAMattributes` set.
+    ///
+    /// `primary_hit_idx` (0-based) is the projected alignment selected as
+    /// primary (randomly among ties per STAR's `rngUniformReal0to1`).  All
+    /// other records get the SECONDARY flag (0x100).
+    #[allow(clippy::too_many_arguments)]
+    pub fn build_transcriptome_records(
+        read_name: &str,
+        read_seq: &[u8],
+        read_qual: &[u8],
+        projected: &[Transcript],
+        mapq: u8,
+        params: &Parameters,
+        primary_hit_idx: usize,
+    ) -> Result<Vec<RecordBuf>, Error> {
+        if projected.is_empty() {
+            return Ok(Vec::new());
+        }
+        let mut attrs = params.sam_attribute_set();
+        // Splice tags are meaningless in t-space.
+        attrs.remove("jM");
+        attrs.remove("jI");
+        attrs.remove("XS");
+        // MD-tag would require the transcript's t-space reference which we do
+        // not precompute; drop it to keep this writer simple.  STAR also does
+        // not emit MD for transcriptome SAM.
+        attrs.remove("MD");
+
+        let n_alignments = projected.len();
+        let mut records = Vec::with_capacity(n_alignments);
+
+        for (hit_idx, t) in projected.iter().enumerate() {
+            let mut record = RecordBuf::default();
+            record.name_mut().replace(read_name.into());
+
+            // FLAGS: SECONDARY if not the primary; REVERSE if is_reverse.
+            let mut flags = sam::alignment::record::Flags::empty();
+            if t.is_reverse {
+                flags |= sam::alignment::record::Flags::REVERSE_COMPLEMENTED;
+            }
+            if hit_idx != primary_hit_idx {
+                flags |= sam::alignment::record::Flags::SECONDARY;
+            }
+            *record.flags_mut() = flags;
+
+            // RNAME = transcript index (maps to transcriptome header).
+            *record.reference_sequence_id_mut() = Some(t.chr_idx);
+
+            // POS = t-space position + 1 (1-based).
+            let pos = (t.genome_start + 1) as usize;
+            *record.alignment_start_mut() = Some(pos.try_into().map_err(|e| {
+                Error::Alignment(format!("invalid t-space position {}: {}", pos, e))
+            })?);
+
+            // MAPQ
+            *record.mapping_quality_mut() = MappingQuality::new(mapq);
+
+            // CIGAR (already has N ops stripped by align_to_transcripts)
+            *record.cigar_mut() = convert_cigar(&t.cigar)?;
+
+            // SEQ / QUAL — STAR writes the original-orientation sequence when
+            // FLAG 0x10 is unset (forward alignment in t-space) and RC'd seq
+            // when 0x10 is set.  We follow SAM spec: SEQ matches the CIGAR's
+            // read orientation, so we mirror `transcript_to_record`.
+            if t.is_reverse {
+                let seq_bytes: Vec<u8> = read_seq
+                    .iter()
+                    .rev()
+                    .map(|&b| decode_base(complement_base(b)))
+                    .collect();
+                *record.sequence_mut() = Sequence::from(seq_bytes);
+                let mut qual = read_qual.to_vec();
+                qual.reverse();
+                *record.quality_scores_mut() = QualityScores::from(qual);
+            } else {
+                let seq_bytes: Vec<u8> = read_seq.iter().map(|&b| decode_base(b)).collect();
+                *record.sequence_mut() = Sequence::from(seq_bytes);
+                *record.quality_scores_mut() = QualityScores::from(read_qual.to_vec());
+            }
+
+            // Optional tags
+            let data = record.data_mut();
+            if attrs.contains("NH") {
+                data.insert(Tag::ALIGNMENT_HIT_COUNT, Value::from(n_alignments as i32));
+            }
+            if attrs.contains("HI") {
+                // HI is 1-based; primary = 1, secondaries > 1 in emission order.
+                data.insert(Tag::HIT_INDEX, Value::from((hit_idx + 1) as i32));
+            }
+            if attrs.contains("AS") {
+                data.insert(Tag::ALIGNMENT_SCORE, Value::from(t.score));
+            }
+            if attrs.contains("NM") || attrs.contains("nM") {
+                data.insert(Tag::new(b'n', b'M'), Value::from(t.n_mismatch as i32));
+            }
+
+            records.push(record);
+        }
+
+        Ok(records)
+    }
+
     /// Build unmapped paired records (both mates unmapped)
     pub fn build_paired_unmapped_records(
         read_name: &str,
@@ -600,21 +712,33 @@ impl SamWriter {
 
 /// Build paired SAM header from genome
 pub fn build_sam_header(genome: &Genome, params: &Parameters) -> Result<sam::Header, Error> {
+    build_sam_header_from_refs(
+        (0..genome.n_chr_real)
+            .map(|i| (genome.chr_name[i].as_str(), genome.chr_length[i] as usize)),
+        params,
+    )
+}
+
+/// Build a SAM header from an iterator of (name, length) reference pairs.
+///
+/// Used both for the genome header (chromosomes) and the transcriptome header
+/// (one @SQ per transcript, length = transcript-space length).
+pub fn build_sam_header_from_refs<'a, I>(refs: I, params: &Parameters) -> Result<sam::Header, Error>
+where
+    I: IntoIterator<Item = (&'a str, usize)>,
+{
     let mut builder = sam::Header::builder();
 
     // @HD line (default version and unsorted)
     builder = builder.set_header(Default::default());
 
-    // @SQ lines for each chromosome
-    for i in 0..genome.n_chr_real {
-        let name = &genome.chr_name[i];
-        let length = genome.chr_length[i] as usize;
-
+    // @SQ lines for each reference
+    for (name, length) in refs {
         let length_nz = NonZeroUsize::new(length)
-            .ok_or_else(|| Error::Index(format!("chromosome {} has zero length", name)))?;
+            .ok_or_else(|| Error::Index(format!("reference {} has zero length", name)))?;
 
         builder = builder.add_reference_sequence(
-            name.as_str(),
+            name,
             Map::<sam::header::record::value::map::ReferenceSequence>::new(length_nz),
         );
     }