diff --git a/CHANGELOG.md b/CHANGELOG.md
index b6e4893..2c40a09 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -97,6 +97,32 @@ Sections commonly used: Features, Bug fixes, Other changes.
   glibc's malloc and per-thread heaps that return whole segments to
   the OS when abandoned, so allocator cache size stays bounded.
 
+- **`--genomeTransformOutput SAM`** reports alignments in the original
+  genome's coordinates. `--genomeTransformType Haploid` bakes a VCF's
+  variants into the genome, so reads carrying those alleles align
+  without mismatches, but every reported coordinate then refers to a
+  genome nobody else has. This maps each alignment back through the
+  conversion blocks written at build time: an indel baked into the
+  sequence reappears as an `I`/`D` CIGAR operation at the original
+  position, and junction motifs are reclassified against the original
+  genome rather than the transformed one. The SAM header comes from
+  the original genome too. `SJ` and `Quant` remain unimplemented and
+  are refused, as are the flag combinations whose other outputs would
+  stay in transformed coordinates.
+
+- **`--genomeTransformOutput SAM`** reports alignments in the original
+  genome's coordinates. `--genomeTransformType Haploid` bakes a VCF's
+  variants into the genome, so reads carrying those alleles align
+  without mismatches, but every reported coordinate then refers to a
+  genome nobody else has. This maps each alignment back through the
+  conversion blocks written at build time: an indel baked into the
+  sequence reappears as an `I`/`D` CIGAR operation at the original
+  position, and junction motifs are reclassified against the original
+  genome. The SAM header comes from the original genome too. `SJ` and
+  `Quant` remain unimplemented and are refused, as are the flag
+  combinations whose other outputs would stay in transformed
+  coordinates.
+
 ### Bug fixes
 
 - **STARsolo `Gene` assignment now requires exon concordance**, matching
diff --git a/src/align/read_align.rs b/src/align/read_align.rs
index b4d2f20..33a2e97 100644
--- a/src/align/read_align.rs
+++ b/src/align/read_align.rs
@@ -1616,6 +1616,7 @@ mod tests {
             transcriptome: None,
             prepared_junctions: Vec::new(),
             sjdb_overhang: 0,
+            transform_out: None,
         }
     }
 
diff --git a/src/align/stitch.rs b/src/align/stitch.rs
index 1c25ecc..97d97e7 100644
--- a/src/align/stitch.rs
+++ b/src/align/stitch.rs
@@ -3191,6 +3191,7 @@ mod tests {
             transcriptome: None,
             prepared_junctions: Vec::new(),
             sjdb_overhang: 0,
+            transform_out: None,
         }
     }
 
@@ -3312,6 +3313,7 @@ mod tests {
             transcriptome: None,
             prepared_junctions: Vec::new(),
             sjdb_overhang: 0,
+            transform_out: None,
         }
     }
 
diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs
index 4997049..af23d35 100644
--- a/src/chimeric/detect.rs
+++ b/src/chimeric/detect.rs
@@ -1060,6 +1060,7 @@ mod tests {
             transcriptome: None,
             prepared_junctions: Vec::new(),
             sjdb_overhang: 0,
+            transform_out: None,
         }
     }
 
diff --git a/src/genome/mod.rs b/src/genome/mod.rs
index 174ba42..3058413 100644
--- a/src/genome/mod.rs
+++ b/src/genome/mod.rs
@@ -1,5 +1,6 @@
 pub mod fasta;
 pub mod transform;
+pub mod transform_align;
 
 use std::path::Path;
 
diff --git a/src/genome/transform.rs b/src/genome/transform.rs
index 328e4f0..59895c5 100644
--- a/src/genome/transform.rs
+++ b/src/genome/transform.rs
@@ -21,6 +21,7 @@
 use std::collections::BTreeMap;
 use std::fmt::Write as _;
 
+use crate::error::Error;
 use crate::io::fastq::encode_base;
 
 use super::compute_chr_starts;
@@ -328,6 +329,67 @@ pub fn blocks_to_tsv(blocks: &[[u64; 3]]) -> String {
     s
 }
 
+/// The stop entry appended to a parsed block map.
+///
+/// The back-transform walks forward over the blocks an exon reaches into and
+/// stops on the first one that starts past its end. A sentinel that no
+/// coordinate can reach makes that a plain comparison instead of a bounds test
+/// on every iteration, which is how STAR writes it (`genomeOutLoad`).
+pub const BLOCK_SENTINEL: [u64; 3] = [u64::MAX, 0, 0];
+
+/// Parse `transformGenomeBlocks.tsv` back into the conversion map used by the
+/// align-time back-transform.
+///
+/// The file is written for reverse conversion, so its columns are already
+/// `[transformed_start, length, original_start]` — the order
+/// [`crate::genome::transform_align::transform_transcript`] wants, and the
+/// reverse of the `[u64; 3]` order the build side keeps in memory. The result
+/// is sorted by transformed start and terminated by [`BLOCK_SENTINEL`].
+pub fn blocks_from_tsv(text: &str) -> Result<Vec<[u64; 3]>, Error> {
+    let mut lines = text.lines();
+    let header = lines
+        .next()
+        .ok_or_else(|| Error::Index("transformGenomeBlocks.tsv is empty".to_string()))?;
+    let declared: usize = header
+        .split_whitespace()
+        .next()
+        .and_then(|n| n.parse().ok())
+        .ok_or_else(|| {
+            Error::Index(format!(
+                "transformGenomeBlocks.tsv: bad header line {header:?}"
+            ))
+        })?;
+
+    let mut blocks = Vec::with_capacity(declared + 1);
+    for line in lines {
+        if line.trim().is_empty() {
+            continue;
+        }
+        let mut f = line.split_whitespace();
+        let mut next = |what: &str| -> Result<u64, Error> {
+            f.next().and_then(|v| v.parse().ok()).ok_or_else(|| {
+                Error::Index(format!("transformGenomeBlocks.tsv: bad {what} in {line:?}"))
+            })
+        };
+        let new_start = next("new_start")?;
+        let length = next("length")?;
+        let orig_start = next("orig_start")?;
+        blocks.push([new_start, length, orig_start]);
+    }
+
+    if blocks.len() != declared {
+        return Err(Error::Index(format!(
+            "transformGenomeBlocks.tsv declares {declared} blocks but has {}",
+            blocks.len()
+        )));
+    }
+    // The map is written in chromosome order, which is already ascending by
+    // transformed start; sorting is a cheap guarantee rather than a fix.
+    blocks.sort_unstable();
+    blocks.push(BLOCK_SENTINEL);
+    Ok(blocks)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -474,4 +536,24 @@ chr1\t6\t.\tG\tA\t.\t.\t.\tGT\t1/1
         assert_eq!(t.blocks[1][0], 0); // same shared orig_start, not offset
         assert_eq!(t.blocks[1][2], offset); // new_start shifted by hap0's padded genome length
     }
+
+    #[test]
+    fn block_map_round_trips_through_the_tsv() {
+        // In-memory order is [orig_start, length, new_start]; the file stores
+        // the reverse, which is also what the back-transform reads.
+        let written = vec![[0u64, 50, 0], [55, 100, 50]];
+        let parsed = blocks_from_tsv(&blocks_to_tsv(&written)).unwrap();
+        assert_eq!(
+            parsed,
+            vec![[0, 50, 0], [50, 100, 55], BLOCK_SENTINEL],
+            "the parsed map is keyed on transformed coordinates, plus the stop entry"
+        );
+    }
+
+    #[test]
+    fn a_block_count_that_disagrees_with_the_body_is_an_error() {
+        let text = "3\t-1\n0\t50\t0\n50\t100\t55\n";
+        let err = blocks_from_tsv(text).unwrap_err().to_string();
+        assert!(err.contains("declares 3 blocks but has 2"), "{err}");
+    }
 }
diff --git a/src/genome/transform_align.rs b/src/genome/transform_align.rs
new file mode 100644
index 0000000..95c3fa2
--- /dev/null
+++ b/src/genome/transform_align.rs
@@ -0,0 +1,499 @@
+//! Align-time back-transform (`--genomeTransformOutput SAM`).
+//!
+//! STAR `Transcript_transformGenome.cpp` + `ReadAlign_transformGenome.cpp`.
+//!
+//! `--genomeTransformType Haploid` bakes a VCF's variants into the genome
+//! sequence, so reads carrying those alleles align without mismatches. The
+//! price is that every coordinate in the output then refers to a genome nobody
+//! else has. `--genomeTransformOutput SAM` pays it back: each alignment is
+//! mapped through the conversion blocks written at build time
+//! (`transformGenomeBlocks.tsv`) onto the original genome, so an indel baked
+//! into the transformed sequence reappears as an `I`/`D` CIGAR operation at the
+//! original coordinates.
+//!
+//! Three things happen, in STAR's order:
+//!
+//! 1. **Remap.** Each exon is looked up in the block map and split wherever it
+//!    straddles a block boundary, since the two halves land at unrelated
+//!    original coordinates.
+//! 2. **Merge.** Splits that turn out to be seamless in the original genome are
+//!    collapsed again, and a block gap that is shorter on one side than the
+//!    other has its shared part folded back into the neighbouring exon. What
+//!    survives as a gap is a real indel.
+//! 3. **Reclassify.** Junction motifs and annotation flags are recomputed
+//!    against the *original* genome. A junction that was canonical in the
+//!    transformed genome need not be canonical in the original one, and STAR
+//!    reports what the original genome says.
+//!
+//! Only the SAM path is implemented here. `--genomeTransformOutput SJ` and
+//! `Quant` are still rejected by parameter validation.
+
+use noodles::sam::alignment::record::cigar::{self, op::Kind};
+
+use crate::align::score::SpliceMotif;
+use crate::align::transcript::{Exon, Transcript};
+use crate::genome::Genome;
+use crate::junction::SpliceJunctionDb;
+
+/// Map a transcript from transformed-genome coordinates back to the original
+/// genome.
+///
+/// `blocks` is the conversion map in `[transformed_start, length,
+/// original_start]` order, ascending by `[0]`, terminated by the sentinel
+/// [`super::transform::BLOCK_SENTINEL`] — the walk over trailing blocks relies
+/// on a stop entry rather than a bounds test, as STAR's does.
+///
+/// Returns `None` when the transcript cannot be converted: no block covers its
+/// start, or it ends up empty.
+pub fn transform_transcript(
+    orig: &Genome,
+    orig_junctions: &SpliceJunctionDb,
+    blocks: &[[u64; 3]],
+    tr: &Transcript,
+    align_intron_min: u64,
+    align_intron_max: u64,
+) -> Option<Transcript> {
+    if tr.exons.is_empty() || blocks.is_empty() {
+        return None;
+    }
+
+    let exons = remap_exons(blocks, &tr.exons)?;
+    let exons = merge_adjacent(exons);
+
+    let (motifs, annotated) = reclassify_junctions(
+        orig,
+        orig_junctions,
+        &exons,
+        tr.is_reverse,
+        align_intron_min,
+    );
+
+    let (chr_idx, _) = orig.position_to_chr(exons[0].genome_start)?;
+    let cigar = rebuild_cigar(tr, &exons, align_intron_min, align_intron_max);
+
+    let genome_start = exons[0].genome_start;
+    let genome_end = exons[exons.len() - 1].genome_end;
+    let n_junction = motifs.len() as u32;
+    // Gaps that are not junctions are indels. Counting them off the rebuilt
+    // CIGAR rather than the exon list keeps the two consistent: the CIGAR is
+    // what decides which gap is `N` and which is `D`.
+    let n_gap = cigar
+        .iter()
+        .filter(|op| matches!(op.kind(), Kind::Insertion | Kind::Deletion))
+        .count() as u32;
+
+    Some(Transcript {
+        chr_idx,
+        genome_start,
+        genome_end,
+        is_reverse: tr.is_reverse,
+        exons,
+        cigar,
+        score: tr.score,
+        n_mismatch: tr.n_mismatch,
+        n_gap,
+        n_junction,
+        junction_motifs: motifs,
+        junction_annotated: annotated,
+        read_seq: tr.read_seq.clone(),
+    })
+}
+
+/// Step 1: per-exon remap, splitting at block boundaries.
+fn remap_exons(blocks: &[[u64; 3]], exons: &[Exon]) -> Option<Vec<Exon>> {
+    let mut out: Vec<Exon> = Vec::with_capacity(exons.len());
+    for exon in exons {
+        let len = exon.genome_end - exon.genome_start;
+        if len == 0 {
+            continue;
+        }
+        let g1 = exon.genome_start;
+        let g2 = exon.genome_end - 1;
+        let read_start = exon.read_start;
+
+        // The last block starting at or before this exon. STAR reaches it with
+        // `upper_bound` then `--`; with nothing at or before it, the exon lies
+        // outside the map entirely and the alignment does not convert.
+        let idx = blocks.partition_point(|b| b[0] <= g1);
+        if idx == 0 {
+            return None;
+        }
+        let mut ci = idx - 1;
+
+        let b_start = blocks[ci][0];
+        let b_end = blocks[ci][0] + blocks[ci][1] - 1;
+        if g1 <= b_end {
+            let piece = if g2 <= b_end { len } else { b_end - g1 + 1 };
+            out.push(Exon {
+                genome_start: blocks[ci][2] + g1 - b_start,
+                genome_end: blocks[ci][2] + g1 - b_start + piece,
+                read_start,
+                read_end: read_start + piece as usize,
+                i_frag: exon.i_frag,
+            });
+        }
+
+        // Any further blocks this exon reaches into, each a separate piece.
+        ci += 1;
+        while ci < blocks.len() && g2 >= blocks[ci][0] {
+            let piece = if g2 < blocks[ci][0] + blocks[ci][1] {
+                g2 - blocks[ci][0] + 1
+            } else {
+                blocks[ci][1]
+            };
+            let r = read_start + (blocks[ci][0] - g1) as usize;
+            out.push(Exon {
+                genome_start: blocks[ci][2],
+                genome_end: blocks[ci][2] + piece,
+                read_start: r,
+                read_end: r + piece as usize,
+                i_frag: exon.i_frag,
+            });
+            ci += 1;
+        }
+    }
+    if out.is_empty() { None } else { Some(out) }
+}
+
+/// Step 2: collapse the splits that are seamless in the original genome, and
+/// fold the shared part of an uneven gap back into the exons around it.
+fn merge_adjacent(exons: Vec<Exon>) -> Vec<Exon> {
+    let mut out: Vec<Exon> = Vec::with_capacity(exons.len());
+    for exon in exons {
+        let Some(prev) = out.last_mut() else {
+            out.push(exon);
+            continue;
+        };
+        if prev.i_frag != exon.i_frag {
+            out.push(exon);
+            continue;
+        }
+        let gap_r = (exon.read_start - prev.read_end) as u64;
+        let gap_g = exon.genome_start - prev.genome_end;
+        if gap_r == gap_g {
+            // Same gap on both sides: not an indel, so the split was an
+            // artefact of the block boundary. Absorb it, gap included.
+            prev.genome_end = exon.genome_end;
+            prev.read_end = exon.read_end;
+        } else {
+            // Uneven gap: the smaller side is aligned sequence, not part of the
+            // indel. Give it to the following exon so the remaining gap is the
+            // indel alone.
+            let shared = gap_r.min(gap_g);
+            let mut exon = exon;
+            if shared > 0 {
+                exon.genome_start -= shared;
+                exon.read_start -= shared as usize;
+            }
+            out.push(exon);
+        }
+    }
+    out
+}
+
+/// Step 3: junction motifs and annotation flags, read off the original genome.
+fn reclassify_junctions(
+    orig: &Genome,
+    orig_junctions: &SpliceJunctionDb,
+    exons: &[Exon],
+    is_reverse: bool,
+    align_intron_min: u64,
+) -> (Vec<SpliceMotif>, Vec<bool>) {
+    let mut motifs = Vec::new();
+    let mut annotated = Vec::new();
+    for pair in exons.windows(2) {
+        let (a, b) = (&pair[0], &pair[1]);
+        if a.i_frag != b.i_frag {
+            continue; // the mate gap is not a junction
+        }
+        let gap_g = b.genome_start.saturating_sub(a.genome_end);
+        let gap_r = b.read_start - a.read_end;
+        if gap_r > 0 || gap_g < align_intron_min {
+            continue; // insertion or deletion, not a junction
+        }
+        let motif = junction_motif(orig, a.genome_end, b.genome_start - 1);
+        let is_annot = orig
+            .position_to_chr(a.genome_end)
+            .is_some_and(|(chr, offset)| {
+                let end = b.genome_start - 1 - orig.chr_start[chr];
+                let strand = if is_reverse { 2 } else { 1 };
+                orig_junctions.is_annotated(chr, offset, end, strand)
+                    || orig_junctions.is_annotated(chr, offset, end, 0)
+            });
+        motifs.push(motif);
+        annotated.push(is_annot);
+    }
+    (motifs, annotated)
+}
+
+/// The donor/acceptor dinucleotides of intron `[j_s, j_e]` in the original
+/// genome (STAR `Transcript_transformGenome.cpp:144-156`).
+fn junction_motif(orig: &Genome, j_s: u64, j_e: u64) -> SpliceMotif {
+    let at = |i: u64| orig.get_base(i).unwrap_or(5);
+    match (at(j_s), at(j_s + 1), at(j_e.wrapping_sub(1)), at(j_e)) {
+        (2, 3, 0, 2) => SpliceMotif::GtAg,
+        (1, 3, 0, 1) => SpliceMotif::CtAc,
+        (2, 1, 0, 2) => SpliceMotif::GcAg,
+        (1, 3, 2, 1) => SpliceMotif::CtGc,
+        (0, 3, 0, 1) => SpliceMotif::AtAc,
+        (2, 3, 0, 3) => SpliceMotif::GtAt,
+        _ => SpliceMotif::NonCanonical,
+    }
+}
+
+/// Rebuild the CIGAR from the remapped exons.
+///
+/// The soft clips are the original transcript's: back-transforming moves where
+/// the read aligns, never how much of it aligns.
+fn rebuild_cigar(
+    tr: &Transcript,
+    exons: &[Exon],
+    align_intron_min: u64,
+    align_intron_max: u64,
+) -> Vec<cigar::Op> {
+    let mut ops: Vec<cigar::Op> = Vec::new();
+    let push_match = |ops: &mut Vec<cigar::Op>, len: usize| {
+        if len == 0 {
+            return;
+        }
+        match ops.last_mut() {
+            Some(op) if op.kind() == Kind::Match => {
+                *op = cigar::Op::new(Kind::Match, op.len() + len);
+            }
+            _ => ops.push(cigar::Op::new(Kind::Match, len)),
+        }
+    };
+
+    let [left_clip, right_clip] = tr.count_soft_clips();
+    if left_clip > 0 {
+        ops.push(cigar::Op::new(Kind::SoftClip, left_clip));
+    }
+    for (i, exon) in exons.iter().enumerate() {
+        if i > 0 {
+            let prev = &exons[i - 1];
+            let gap_r = exon.read_start - prev.read_end;
+            let gap_g = exon.genome_start - prev.genome_end;
+            if gap_r > 0 {
+                ops.push(cigar::Op::new(Kind::Insertion, gap_r));
+            }
+            if gap_g > 0 {
+                let kind = if gap_g >= align_intron_min && gap_g <= align_intron_max {
+                    Kind::Skip
+                } else {
+                    Kind::Deletion
+                };
+                ops.push(cigar::Op::new(kind, gap_g as usize));
+            }
+        }
+        push_match(&mut ops, exon.read_end - exon.read_start);
+    }
+    if right_clip > 0 {
+        ops.push(cigar::Op::new(Kind::SoftClip, right_clip));
+    }
+    ops
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::genome::transform::BLOCK_SENTINEL;
+
+    /// A genome of one chromosome whose bases are given as codes 0..3.
+    fn genome_from_codes(codes: &[u8]) -> Genome {
+        let n = codes.len() as u64;
+        let mut seq = vec![0u8; 2 * n as usize];
+        seq[..codes.len()].copy_from_slice(codes);
+        Genome {
+            transform_blocks: None,
+            sequence: seq.into(),
+            n_genome: n,
+            n_genome_real: n,
+            n_chr_real: 1,
+            chr_name: vec!["chr1".to_string()],
+            chr_length: vec![n],
+            chr_start: vec![0, n],
+        }
+    }
+
+    fn exon(genome_start: u64, len: u64, read_start: usize) -> Exon {
+        Exon {
+            genome_start,
+            genome_end: genome_start + len,
+            read_start,
+            read_end: read_start + len as usize,
+            i_frag: 0,
+        }
+    }
+
+    fn transcript(exons: Vec<Exon>, cigar: Vec<cigar::Op>) -> Transcript {
+        let genome_start = exons[0].genome_start;
+        let genome_end = exons[exons.len() - 1].genome_end;
+        Transcript {
+            chr_idx: 0,
+            genome_start,
+            genome_end,
+            is_reverse: false,
+            exons,
+            cigar,
+            score: 0,
+            n_mismatch: 0,
+            n_gap: 0,
+            n_junction: 0,
+            junction_motifs: vec![],
+            junction_annotated: vec![],
+            read_seq: vec![],
+        }
+    }
+
+    /// Blocks for a genome where 5 original bases were deleted at original
+    /// position 50: transformed [0,50) is original [0,50), and transformed
+    /// [50,..) is original [55,..).
+    fn deletion_blocks() -> Vec<[u64; 3]> {
+        vec![[0, 50, 0], [50, 100, 55], BLOCK_SENTINEL]
+    }
+
+    #[test]
+    fn an_exon_inside_one_block_is_shifted_by_that_block() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        let tr = transcript(vec![exon(60, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]);
+        let out =
+            transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap();
+        // transformed 60 is 10 into the second block, which starts at original 55.
+        assert_eq!(out.genome_start, 65);
+        assert_eq!(out.exons.len(), 1);
+        assert_eq!(out.cigar_string(), "20M");
+    }
+
+    /// The point of the whole module: a variant baked into the genome comes
+    /// back out as a CIGAR operation.
+    #[test]
+    fn an_exon_spanning_a_block_boundary_becomes_a_deletion() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        // One 20-base exon straddling the boundary at transformed 50.
+        let tr = transcript(vec![exon(40, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]);
+        let out =
+            transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap();
+        assert_eq!(out.genome_start, 40);
+        assert_eq!(out.exons.len(), 2, "the split survives as a real gap");
+        // 10 bases, the 5 deleted bases, then 10 more.
+        assert_eq!(out.cigar_string(), "10M5D10M");
+        assert_eq!(out.n_gap, 1);
+    }
+
+    /// A gap the transform did not create must not be turned into one: when the
+    /// two pieces are contiguous in the original genome too, the split
+    /// disappears again.
+    #[test]
+    fn a_seamless_split_is_merged_back() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        // Two blocks that are adjacent in both coordinate systems: the split is
+        // pure bookkeeping.
+        let blocks = vec![[0, 50, 0], [50, 100, 50], BLOCK_SENTINEL];
+        let tr = transcript(vec![exon(40, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]);
+        let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap();
+        assert_eq!(out.exons.len(), 1);
+        assert_eq!(out.cigar_string(), "20M");
+        assert_eq!(out.n_gap, 0);
+    }
+
+    /// Soft clips describe the read, not the genome, so they survive unchanged.
+    #[test]
+    fn soft_clips_are_carried_through() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        let tr = transcript(
+            vec![exon(60, 20, 5)],
+            vec![
+                cigar::Op::new(Kind::SoftClip, 5),
+                cigar::Op::new(Kind::Match, 20),
+                cigar::Op::new(Kind::SoftClip, 3),
+            ],
+        );
+        let out =
+            transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap();
+        assert_eq!(out.cigar_string(), "5S20M3S");
+    }
+
+    /// A junction is reclassified against the original genome, not the
+    /// transformed one. Here the original bases spell GT..AG.
+    #[test]
+    fn a_junction_motif_is_read_from_the_original_genome() {
+        let mut codes = vec![0u8; 200];
+        // intron [80, 119]: GT at the donor, AG at the acceptor.
+        codes[80] = 2;
+        codes[81] = 3;
+        codes[118] = 0;
+        codes[119] = 2;
+        let orig = genome_from_codes(&codes);
+        let jdb = SpliceJunctionDb::empty();
+        // Both exons sit in the identity block, so the coordinates survive.
+        let blocks = vec![[0, 200, 0], BLOCK_SENTINEL];
+        let tr = transcript(
+            vec![exon(60, 20, 0), exon(120, 20, 20)],
+            vec![
+                cigar::Op::new(Kind::Match, 20),
+                cigar::Op::new(Kind::Skip, 40),
+                cigar::Op::new(Kind::Match, 20),
+            ],
+        );
+        let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap();
+        assert_eq!(out.junction_motifs, vec![SpliceMotif::GtAg]);
+        assert_eq!(out.cigar_string(), "20M40N20M");
+        assert_eq!(out.n_junction, 1);
+    }
+
+    /// The same geometry over a genome that does not spell a canonical motif
+    /// there. STAR reports what the original genome says, even when the
+    /// transformed genome said otherwise.
+    #[test]
+    fn a_junction_canonical_only_in_the_transformed_genome_becomes_noncanonical() {
+        let orig = genome_from_codes(&[0u8; 200]); // all A: no motif anywhere
+        let jdb = SpliceJunctionDb::empty();
+        let blocks = vec![[0, 200, 0], BLOCK_SENTINEL];
+        let tr = transcript(
+            vec![exon(60, 20, 0), exon(120, 20, 20)],
+            vec![
+                cigar::Op::new(Kind::Match, 20),
+                cigar::Op::new(Kind::Skip, 40),
+                cigar::Op::new(Kind::Match, 20),
+            ],
+        );
+        let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap();
+        assert_eq!(out.junction_motifs, vec![SpliceMotif::NonCanonical]);
+    }
+
+    /// A gap shorter than `--alignIntronMin` is a deletion, not a junction, and
+    /// carries no motif.
+    #[test]
+    fn a_short_gap_is_a_deletion_not_a_junction() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        let blocks = vec![[0, 200, 0], BLOCK_SENTINEL];
+        let tr = transcript(
+            vec![exon(60, 20, 0), exon(85, 20, 20)],
+            vec![
+                cigar::Op::new(Kind::Match, 20),
+                cigar::Op::new(Kind::Deletion, 5),
+                cigar::Op::new(Kind::Match, 20),
+            ],
+        );
+        let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap();
+        assert!(out.junction_motifs.is_empty());
+        assert_eq!(out.n_junction, 0);
+        assert_eq!(out.cigar_string(), "20M5D20M");
+    }
+
+    /// An alignment starting before every block cannot be converted, and is
+    /// dropped rather than guessed at.
+    #[test]
+    fn an_alignment_outside_the_block_map_does_not_convert() {
+        let orig = genome_from_codes(&[0u8; 200]);
+        let jdb = SpliceJunctionDb::empty();
+        let blocks = vec![[100, 50, 100], BLOCK_SENTINEL];
+        let tr = transcript(vec![exon(10, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]);
+        assert!(transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).is_none());
+    }
+}
diff --git a/src/index/io.rs b/src/index/io.rs
index ce6de64..4edff7f 100644
--- a/src/index/io.rs
+++ b/src/index/io.rs
@@ -6,6 +6,7 @@ use byteorder::{LittleEndian, ReadBytesExt};
 use crate::error::Error;
 use crate::genome::Genome;
 use crate::index::GenomeIndex;
+use crate::index::TransformOut;
 use crate::index::packed_array::PackedArray;
 use crate::index::sa_index::SaIndex;
 use crate::index::suffix_array::SuffixArray;
@@ -138,6 +139,15 @@ impl GenomeIndex {
             );
         }
 
+        // `--genomeTransformOutput SAM`: the original genome and the block map
+        // needed to report alignments in its coordinates rather than the
+        // transformed ones the search runs against.
+        let transform_out = if params.transform_output_sam() {
+            Some(load_transform_out(genome_dir, params)?)
+        } else {
+            None
+        };
+
         Ok(GenomeIndex {
             genome,
             suffix_array,
@@ -146,10 +156,71 @@ impl GenomeIndex {
             transcriptome,
             prepared_junctions,
             sjdb_overhang,
+            transform_out,
         })
     }
 }
 
+/// Load everything `--genomeTransformOutput SAM` needs: the untransformed
+/// genome written to `OriginalGenome/` at build time, its annotated junctions,
+/// and the `transformGenomeBlocks.tsv` conversion map.
+///
+/// Both files exist only in an index built with `--genomeTransformType`, so a
+/// missing one means the flag was asked for against an ordinary index. That is
+/// an error rather than a silent fallback: falling back would report
+/// transformed coordinates while claiming they are original.
+fn load_transform_out(genome_dir: &Path, params: &Parameters) -> Result<TransformOut, Error> {
+    let blocks_path = genome_dir.join("transformGenomeBlocks.tsv");
+    let text = std::fs::read_to_string(&blocks_path).map_err(|e| {
+        if e.kind() == std::io::ErrorKind::NotFound {
+            Error::Index(format!(
+                "--genomeTransformOutput SAM needs {}, which is only written by a \
+                 --genomeTransformType genomeGenerate run",
+                blocks_path.display()
+            ))
+        } else {
+            Error::io(e, &blocks_path)
+        }
+    })?;
+    let blocks = crate::genome::transform::blocks_from_tsv(&text)?;
+
+    let orig_dir = genome_dir.join("OriginalGenome");
+    if !orig_dir.join("chrName.txt").exists() {
+        return Err(Error::Index(format!(
+            "--genomeTransformOutput SAM needs the untransformed index in {}",
+            orig_dir.display()
+        )));
+    }
+    let genome = load_genome(&orig_dir, params)?;
+
+    // The original genome's annotated junctions, for the motif/annotation flags
+    // recomputed after the back-transform. Absent when the index was built
+    // without a GTF, in which case every junction is novel there too.
+    let sjdb_info_path = orig_dir.join("sjdbInfo.txt");
+    let junctions = if sjdb_info_path.exists() {
+        let tab = sjdb_insert::read_sjdb_info_tab(&sjdb_info_path, &genome)?;
+        let raw: Vec<(usize, u64, u64, u8)> = tab
+            .junctions
+            .iter()
+            .map(|j| (j.chr_idx, j.stored_start(), j.stored_end(), j.strand))
+            .collect();
+        SpliceJunctionDb::from_raw_junctions(&raw)
+    } else {
+        SpliceJunctionDb::empty()
+    };
+
+    log::info!(
+        "genomeTransformOutput SAM: original genome {} chromosomes, {} conversion blocks",
+        genome.n_chr_real,
+        blocks.len() - 1, // the sentinel is not a block
+    );
+    Ok(TransformOut {
+        genome,
+        junctions,
+        blocks,
+    })
+}
+
 /// Read `genomeFileSizes\t<n_genome> <sa_size>` from genomeParameters.txt
 /// and return the first field (total genome byte count, including Gsj if
 /// sjdb was baked in). Returns `Ok(None)` if the file or line is absent,
diff --git a/src/index/mod.rs b/src/index/mod.rs
index 88ee62d..bd7296c 100644
--- a/src/index/mod.rs
+++ b/src/index/mod.rs
@@ -42,6 +42,59 @@ pub struct GenomeIndex {
     /// `sjdbOverhang` recorded in `sjdbInfo.txt`. Zero when no sjdb
     /// junctions are present.
     pub sjdb_overhang: u32,
+    /// Populated only for `--genomeTransformOutput SAM`: what the output
+    /// coordinate space is, when it is not the one the search runs against.
+    pub transform_out: Option<TransformOut>,
+}
+
+/// The original genome an alignment is reported against under
+/// `--genomeTransformOutput SAM`.
+///
+/// Reads are searched against the transformed genome, since that is what the
+/// suffix array indexes, and then mapped back through `blocks`. The SAM header
+/// comes from the genome here, not from the searched one: reporting original
+/// coordinates under transformed reference lengths would produce a file no
+/// downstream tool could read correctly.
+#[derive(Clone)]
+pub struct TransformOut {
+    /// The untransformed genome, from `OriginalGenome/`.
+    pub genome: Genome,
+    /// Its annotated junctions, for the motifs recomputed after the transform.
+    pub junctions: SpliceJunctionDb,
+    /// `[transformed_start, length, original_start]`, ascending, sentinel-terminated.
+    pub blocks: Vec<[u64; 3]>,
+}
+
+impl GenomeIndex {
+    /// The genome alignments are *reported* against, which is the searched one
+    /// unless `--genomeTransformOutput SAM` moved output into the original
+    /// coordinate space.
+    pub fn output_genome(&self) -> &Genome {
+        match &self.transform_out {
+            Some(t) => &t.genome,
+            None => &self.genome,
+        }
+    }
+
+    /// Map a transcript into the output coordinate space, dropping it when it
+    /// does not convert. Without a transform this hands back what it was given.
+    pub fn to_output_space(
+        &self,
+        tr: crate::align::transcript::Transcript,
+        params: &Parameters,
+    ) -> Option<crate::align::transcript::Transcript> {
+        let Some(t) = &self.transform_out else {
+            return Some(tr);
+        };
+        crate::genome::transform_align::transform_transcript(
+            &t.genome,
+            &t.junctions,
+            &t.blocks,
+            &tr,
+            params.align_intron_min as u64,
+            params.align_intron_max as u64,
+        )
+    }
 }
 
 /// Output of [`GenomeIndex::build_prep`] — the shared setup
@@ -89,6 +142,9 @@ impl GenomeIndex {
             sa_index,
             junction_db,
             transcriptome,
+            // Built here, not loaded: the back-transform is an align-time
+            // concern and `build` has no index directory to read it from.
+            transform_out: None,
             prepared_junctions,
             sjdb_overhang,
         })
diff --git a/src/lib.rs b/src/lib.rs
index 6fce173..3e2ae9c 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -675,21 +675,21 @@ fn run_single_pass(
         OutStd::Sam => {
             info!("Writing SAM to stdout (--outStd SAM)");
             Box::new(crate::io::sam::SamStdoutWriter::create(
-                &index.genome,
+                index.output_genome(),
                 params,
             )?)
         }
         OutStd::BamUnsorted => {
             info!("Writing unsorted BAM to stdout (--outStd BAM_Unsorted)");
             Box::new(crate::io::bam::BamStdoutWriter::create(
-                &index.genome,
+                index.output_genome(),
                 params,
             )?)
         }
         OutStd::BamSortedByCoordinate => {
             info!("Writing coordinate-sorted BAM to stdout (--outStd BAM_SortedByCoordinate)");
             Box::new(crate::io::bam::SortedBamStdoutWriter::create(
-                &index.genome,
+                index.output_genome(),
                 params,
             )?)
         }
@@ -700,7 +700,11 @@ fn run_single_pass(
                 if let Some(parent) = output_path.parent() {
                     std::fs::create_dir_all(parent)?;
                 }
-                Box::new(SamWriter::create(&output_path, &index.genome, params)?)
+                Box::new(SamWriter::create(
+                    &output_path,
+                    index.output_genome(),
+                    params,
+                )?)
             }
             OutSamFormat::Bam => {
                 let sorted = out_type.sort_order == Some(OutSamSortOrder::SortedByCoordinate);
@@ -716,11 +720,15 @@ fn run_single_pass(
                 if sorted {
                     Box::new(SortedBamWriter::create(
                         &output_path,
-                        &index.genome,
+                        index.output_genome(),
                         params,
                     )?)
                 } else {
-                    Box::new(BamWriter::create(&output_path, &index.genome, params)?)
+                    Box::new(BamWriter::create(
+                        &output_path,
+                        index.output_genome(),
+                        params,
+                    )?)
                 }
             }
             OutSamFormat::None => {
@@ -750,7 +758,7 @@ fn run_single_pass(
         }
         let sj_output_path = params.output_path("SJ.out.tab");
         if !sj_stats.is_empty() {
-            sj_stats.write_output(&sj_output_path, &index.genome, params)?;
+            sj_stats.write_output(&sj_output_path, index.output_genome(), params)?;
         }
         // Per-cell count matrices (raw + filtered), Summary.csv, and the SJ
         // feature matrix — written here where sj_stats is available.
@@ -801,7 +809,7 @@ fn run_single_pass(
             "Writing splice junction statistics to {}",
             sj_output_path.display()
         );
-        sj_stats.write_output(&sj_output_path, &index.genome, params)?;
+        sj_stats.write_output(&sj_output_path, index.output_genome(), params)?;
     }
 
     // 6. Print summary
@@ -829,7 +837,7 @@ fn run_two_pass(
     let pass1_path = pass1_dir.join("SJ.out.tab");
 
     info!("Writing pass 1 junctions to {}", pass1_path.display());
-    sj_stats_pass1.write_output(&pass1_path, &index.genome, params)?;
+    sj_stats_pass1.write_output(&pass1_path, index.output_genome(), params)?;
     info!(
         "Pass 1 discovered {} novel junctions",
         novel_junctions.len()
@@ -1273,8 +1281,11 @@ fn extract_junction_keys(
                 let intron_start = genome_pos;
                 let intron_end = genome_pos + intron_len as u64 - 1;
 
-                let motif =
-                    scorer.detect_splice_motif(genome_pos, intron_len as u32, &index.genome);
+                let motif = scorer.detect_splice_motif(
+                    genome_pos,
+                    intron_len as u32,
+                    index.output_genome(),
+                );
                 let strand = match motif.implied_strand() {
                     Some('+') => 1u8,
                     Some('-') => 2u8,
@@ -1413,7 +1424,8 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
         let write_file = tf
             .reopen()
             .map_err(|e| anyhow::anyhow!("BySJout: temp file reopen error: {e}"))?;
-        let (hdr, w) = crate::io::sam::create_bysj_writer(write_file, &index.genome, params)?;
+        let (hdr, w) =
+            crate::io::sam::create_bysj_writer(write_file, index.output_genome(), params)?;
         (Some(hdr), Some(w))
     } else {
         (None, None)
@@ -1806,6 +1818,21 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                         let (transcripts, chimeric_results, n_for_mapq, unmapped_reason) =
                             align_read(&clipped_seq, &read.name, &index, params)?;
 
+                        // `--genomeTransformOutput SAM`: the search ran against the
+                        // transformed genome; everything downstream — records,
+                        // junction counts, statistics — is in the original genome's
+                        // coordinates. Converting here rather than at the writer
+                        // keeps a single space below this line. An alignment that
+                        // does not convert is dropped, as STAR drops it.
+                        let transcripts: Vec<_> = if index.transform_out.is_some() {
+                            transcripts
+                                .into_iter()
+                                .filter_map(|t| index.to_output_space(t, params))
+                                .collect()
+                        } else {
+                            transcripts
+                        };
+
                         // Collect chimeric alignments if enabled
                         if params.chim_segment_min > 0 {
                             chimeric_alns.extend(chimeric_results);
@@ -1891,7 +1918,7 @@ fn align_reads_single_end<W: AlignmentWriter + ?Sized>(
                                     clip5p,
                                     clip3p,
                                     &transcripts,
-                                    &index.genome,
+                                    index.output_genome(),
                                     params,
                                     n_for_mapq,
                                 )?;
@@ -2637,7 +2664,7 @@ fn align_reads_solo_pe<W: AlignmentWriter + ?Sized>(
                                     clip5p_m2,
                                     clip3p_m2,
                                     &paired_alns,
-                                    &index.genome,
+                                    index.output_genome(),
                                     params,
                                     n_for_mapq,
                                 )?;
@@ -2751,7 +2778,8 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
         let write_file = tf
             .reopen()
             .map_err(|e| anyhow::anyhow!("BySJout: temp file reopen error: {e}"))?;
-        let (hdr, w) = crate::io::sam::create_bysj_writer(write_file, &index.genome, params)?;
+        let (hdr, w) =
+            crate::io::sam::create_bysj_writer(write_file, index.output_genome(), params)?;
         (Some(hdr), Some(w))
     } else {
         (None, None)
@@ -3157,6 +3185,19 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                         let (results, pe_chimeric, n_for_mapq, unmapped_reason) =
                             align_paired_read(&m1_seq, &m2_seq, &paired_read.name, &index, params)?;
 
+                        // `--genomeTransformOutput SAM`: as in the single-end path,
+                        // move to the original genome's coordinates before anything
+                        // reads the transcripts. A pair whose mates do not both
+                        // convert is dropped whole rather than half-reported.
+                        let results: Vec<_> = if index.transform_out.is_some() {
+                            results
+                                .into_iter()
+                                .filter_map(|r| transform_pair_result(r, &index, params))
+                                .collect()
+                        } else {
+                            results
+                        };
+
                         // Classify the result for stats and SAM output
                         let has_half_mapped = results
                             .iter()
@@ -3334,7 +3375,7 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                                     m2_clip3p,
                                     mapped_transcript,
                                     *mate1_is_mapped,
-                                    &index.genome,
+                                    index.output_genome(),
                                     params,
                                     n_for_mapq,
                                 )?;
@@ -3360,7 +3401,7 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                                 m2_clip5p,
                                 m2_clip3p,
                                 &paired_alns,
-                                &index.genome,
+                                index.output_genome(),
                                 params,
                                 n_for_mapq,
                             )?;
@@ -3480,6 +3521,45 @@ struct ReadJunction {
 /// `is_unique` reflects the read's overall mapping multiplicity (n_loci == 1),
 /// not per-locus. Recording per-locus/per-mate instead double-counts junctions
 /// crossed by both mates of a pair, inflating the SJ.out.tab multi counts.
+/// Move one paired result into the output coordinate space under
+/// `--genomeTransformOutput SAM`.
+///
+/// A `BothMapped` pair is kept only if both mates convert: reporting one mate
+/// in original coordinates and dropping the other would turn a pair into a
+/// half-mapped read that never existed. The insert size is recomputed from the
+/// converted mates, since the transform can change the distance between them.
+fn transform_pair_result(
+    result: crate::align::read_align::PairedAlignmentResult,
+    index: &crate::index::GenomeIndex,
+    params: &Parameters,
+) -> Option<crate::align::read_align::PairedAlignmentResult> {
+    use crate::align::read_align::PairedAlignmentResult;
+    match result {
+        PairedAlignmentResult::BothMapped(pa) => {
+            let mut pa = *pa;
+            let m1 = index.to_output_space(pa.mate1_transcript, params)?;
+            let m2 = index.to_output_space(pa.mate2_transcript, params)?;
+            let (lo, hi) = if m1.genome_start <= m2.genome_start {
+                (&m1, &m2)
+            } else {
+                (&m2, &m1)
+            };
+            let span = (hi.genome_end - lo.genome_start) as i32;
+            pa.insert_size = if pa.insert_size < 0 { -span } else { span };
+            pa.mate1_transcript = m1;
+            pa.mate2_transcript = m2;
+            Some(PairedAlignmentResult::BothMapped(Box::new(pa)))
+        }
+        PairedAlignmentResult::HalfMapped {
+            mapped_transcript,
+            mate1_is_mapped,
+        } => Some(PairedAlignmentResult::HalfMapped {
+            mapped_transcript: index.to_output_space(mapped_transcript, params)?,
+            mate1_is_mapped,
+        }),
+    }
+}
+
 fn record_read_junctions<'a>(
     transcripts: impl IntoIterator<Item = &'a crate::align::transcript::Transcript>,
     index: &crate::index::GenomeIndex,
@@ -3561,8 +3641,11 @@ fn extract_transcript_junctions(
                 let intron_end = genome_pos + intron_len as u64 - 1;
 
                 // Detect splice motif
-                let motif =
-                    scorer.detect_splice_motif(genome_pos, intron_len as u32, &index.genome);
+                let motif = scorer.detect_splice_motif(
+                    genome_pos,
+                    intron_len as u32,
+                    index.output_genome(),
+                );
 
                 // Compute overhang: min(left_exon_length, right_exon_length)
                 let left_exon = exon_lengths[junction_idx];
diff --git a/src/params/mod.rs b/src/params/mod.rs
index 0536a85..9774037 100644
--- a/src/params/mod.rs
+++ b/src/params/mod.rs
@@ -508,6 +508,12 @@ pub struct Parameters {
     #[arg(long = "genomeSAsparseD", default_value_t = 1)]
     pub genome_sa_sparse_d: u32,
 
+    /// Coordinate space alignments are reported in when the genome was built
+    /// with `--genomeTransformType`. `None` (the default) reports transformed
+    /// coordinates; `SAM` and `SJ` map them back to the original genome.
+    #[arg(long = "genomeTransformOutput", num_args = 1.., default_values_t = vec!["None".to_string()])]
+    pub genome_transform_output: Vec<String>,
+
     /// Substitute VCF alleles into the genome at genomeGenerate (`None`,
     /// `Haploid`, or `Diploid`). Requires `--genomeTransformVCF`; incompatible
     /// with `--sjdbGTFfile`. `Diploid` is genotype-aware and duplicates the
@@ -1175,6 +1181,12 @@ impl Parameters {
         !matches!(self.out_std, OutStd::None) || self.out_sam_type.format != OutSamFormat::None
     }
 
+    /// Whether `--genomeTransformOutput` asks for SAM records in the original
+    /// genome's coordinates rather than the transformed genome's.
+    pub fn transform_output_sam(&self) -> bool {
+        self.genome_transform_output.iter().any(|o| o == "SAM")
+    }
+
     /// Whether `--chimOutType` includes `Junctions` (write Chimeric.out.junction).
     pub fn chim_out_junctions(&self) -> bool {
         self.chim_out_type.iter().any(|s| s == "Junctions")
@@ -1551,6 +1563,51 @@ impl Parameters {
             ));
         }
 
+        // --genomeTransformOutput: `SAM` maps alignments back to the original
+        // genome's coordinates. `SJ` and `Quant` do the same for the junction
+        // table and the transcriptome BAM and are not implemented, so they are
+        // refused: a silent no-op would report transformed coordinates as if
+        // they were original ones, which is worse than failing.
+        for o in &params.genome_transform_output {
+            if o != "None" && o != "SAM" {
+                return Err(command.error(
+                    ErrorKind::InvalidValue,
+                    format!(
+                        "--genomeTransformOutput {o} is not supported: only SAM output is \
+                         mapped back to original coordinates"
+                    ),
+                ));
+            }
+        }
+        // The back-transform moves the SAM records alone. Everything else that
+        // reports coordinates — the junction table, the transcriptome BAM, the
+        // count matrices, the coverage signal — would still be in transformed
+        // space, and a file that mixes the two spaces is worse than no file.
+        if params.transform_output_sam() {
+            let conflicting = [
+                (!params.quant_mode.iter().all(|m| m == "-"), "--quantMode"),
+                (params.solo_type != SoloType::None, "--soloType"),
+                (
+                    params
+                        .out_wig_type
+                        .iter()
+                        .any(|t| !t.eq_ignore_ascii_case("None")),
+                    "--outWigType",
+                ),
+            ];
+            for (hit, flag) in conflicting {
+                if hit {
+                    return Err(command.error(
+                        ErrorKind::ArgumentConflict,
+                        format!(
+                            "--genomeTransformOutput SAM cannot be combined with {flag}: only \
+                             the SAM records are mapped back to original coordinates"
+                        ),
+                    ));
+                }
+            }
+        }
+
         // ── STARsolo validation ─────────────────────────────────────────
         if params.run_mode == RunMode::AlignReads && params.solo_enabled() {
             // CB_UMI_Complex needs one CB position + whitelist per segment.
@@ -2636,4 +2693,25 @@ mod tests {
         assert_eq!(p.out_sam_strand_field, "None");
         assert!(!p.out_sam_attributes.contains(SamAttributes::XS));
     }
+
+    #[test]
+    fn genome_transform_output_accepts_sam_and_refuses_the_rest() {
+        let gg = |extra: &[&str]| {
+            let mut a = vec!["--runMode", "genomeGenerate", "--genomeFastaFiles", "g.fa"];
+            a.extend_from_slice(extra);
+            try_parse(&a)
+        };
+
+        // Defaults parse, and are the STAR ones.
+        let p = gg(&[]).unwrap();
+        assert_eq!(p.genome_transform_output, vec!["None".to_string()]);
+
+        // SAM output is mapped back to the original genome's coordinates. The
+        // junction table and the transcriptome BAM are not, and a silent no-op
+        // there would report transformed coordinates as if they were original.
+        assert!(gg(&["--genomeTransformOutput", "SAM"]).is_ok());
+        assert!(gg(&["--genomeTransformOutput", "SJ"]).is_err());
+        assert!(gg(&["--genomeTransformOutput", "Quant"]).is_err());
+        assert!(gg(&["--genomeTransformOutput", "None"]).is_ok());
+    }
 }
diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs
index 2db7d39..91a9344 100644
--- a/tests/alignment_features.rs
+++ b/tests/alignment_features.rs
@@ -2024,3 +2024,116 @@ fn test_wasp_samtag() {
         "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)"
     );
 }
+
+// ---------------------------------------------------------------------------
+// --genomeTransformOutput SAM — alignments reported in original coordinates
+// ---------------------------------------------------------------------------
+
+/// A transformed genome carries a deletion the original does not. Reads align
+/// against the transformed genome, where every downstream coordinate has moved
+/// by the deleted length; `--genomeTransformOutput SAM` has to undo that.
+///
+/// Without the back-transform the reported POS is 5 short, which is exactly the
+/// kind of error that looks plausible in a browser and is wrong everywhere.
+#[test]
+fn test_genome_transform_output_sam_reports_original_coordinates() {
+    let tmpdir = TempDir::new().unwrap();
+    let genome = build_genome();
+    let fasta = write_fasta(&tmpdir, &genome);
+
+    // A 5-base deletion at 1-based position 5001: REF spans 6 bases, ALT keeps
+    // the first. Everything after it shifts down by 5 in the transformed genome.
+    let del_pos_1based = 5001usize;
+    let vcf_path = tmpdir.path().join("variants.vcf");
+    {
+        let mut f = fs::File::create(&vcf_path).unwrap();
+        writeln!(f, "##fileformat=VCFv4.2").unwrap();
+        writeln!(f, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO").unwrap();
+        let start = del_pos_1based - 1;
+        let reference = String::from_utf8(genome[start..start + 6].to_vec()).unwrap();
+        let alt = String::from_utf8(genome[start..start + 1].to_vec()).unwrap();
+        writeln!(
+            f,
+            "chr1\t{del_pos_1based}\t.\t{reference}\t{alt}\t.\tPASS\t."
+        )
+        .unwrap();
+    }
+
+    let genome_dir = tmpdir.path().join("genome_transformed");
+    fs::create_dir_all(&genome_dir).unwrap();
+    cargo_bin_cmd!("rustar-aligner")
+        .args([
+            "--runMode",
+            "genomeGenerate",
+            "--genomeDir",
+            genome_dir.to_str().unwrap(),
+            "--genomeFastaFiles",
+            fasta.to_str().unwrap(),
+            "--genomeSAindexNbases",
+            "7",
+            "--genomeTransformType",
+            "Haploid",
+            "--genomeTransformVCF",
+            vcf_path.to_str().unwrap(),
+        ])
+        .assert()
+        .success();
+
+    // One 60 bp read taken from the original genome well downstream of the
+    // deletion. Its sequence is unchanged by the substitution, so it aligns
+    // cleanly in both coordinate systems — only the position differs.
+    let read_start_0based = 6000usize;
+    let fastq_path = tmpdir.path().join("reads.fq");
+    {
+        let mut f = fs::File::create(&fastq_path).unwrap();
+        writeln!(f, "@r1").unwrap();
+        f.write_all(&genome[read_start_0based..read_start_0based + 60])
+            .unwrap();
+        writeln!(f).unwrap();
+        writeln!(f, "+").unwrap();
+        writeln!(f, "{}", "I".repeat(60)).unwrap();
+    }
+
+    let run = |extra: &[&str], out_name: &str| -> Vec<Vec<String>> {
+        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(|l| l.split('\t').map(str::to_string).collect())
+            .collect()
+    };
+
+    let transformed = run(&[], "out_transformed");
+    let original = run(&["--genomeTransformOutput", "SAM"], "out_original");
+
+    assert_eq!(transformed.len(), 1, "expected one alignment");
+    assert_eq!(original.len(), 1, "expected one alignment");
+
+    let pos_transformed: usize = transformed[0][3].parse().unwrap();
+    let pos_original: usize = original[0][3].parse().unwrap();
+
+    // 1-based SAM position in the original genome.
+    assert_eq!(pos_original, read_start_0based + 1);
+    // The transformed genome is 5 bases shorter ahead of this read.
+    assert_eq!(pos_transformed, pos_original - 5);
+    // The read itself spans no variant, so the CIGAR is unaffected.
+    assert_eq!(original[0][5], transformed[0][5]);
+}