solo: --soloFeatures Transcript3p, with --soloClusterCBfile

CHANGELOG.md

@@ -97,6 +97,20 @@ 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.
 
+- **`--soloFeatures Transcript3p`** quantifies transcripts rather than
+  genes, using how far each read's 3' end sits from each transcript's.
+  In a 3'-biased assay that distance discriminates between isoforms: a
+  read 200 bases from the end of one and 4000 from the end of another
+  is evidence for the first. The distance distribution is estimated
+  from the data, then used as the likelihood in an EM over UMIs. Output
+  is per *cluster* rather than per cell — `--soloClusterCBfile` (new,
+  and required for this feature) says which cell is in which cluster,
+  because one cell has too few UMIs to resolve isoforms. Reads sharing
+  a UMI contribute the intersection of their transcript sets, not the
+  union: they came from one molecule. Writes `matrix.mtx`,
+  `features.tsv` and `transcriptEndDistanceDistribution.txt` under
+  `Solo.out/Transcript3p/raw/`.
+
 ### Bug fixes
 
 - **STARsolo `Gene` assignment now requires exon concordance**, matching

src/params/mod.rs

@@ -1113,6 +1113,11 @@ pub struct Parameters {
     #[arg(long = "soloCBmatchWLtype", default_value = "1MM_multi")]
     pub solo_cb_match_wl_type: String,
 
+    /// Two-column `CB cluster` file assigning cells to clusters, for
+    /// `--soloFeatures Transcript3p`.
+    #[arg(long = "soloClusterCBfile")]
+    pub solo_cluster_cb_file: Option<PathBuf>,
+
     /// Cell-calling / matrix filtering: None, CellRanger2.2, EmptyDrops_CR, TopCells.
     #[arg(long = "soloCellFilter", num_args = 1.., default_values_t = vec!["CellRanger2.2".to_string(), "3000".to_string(), "0.99".to_string(), "10".to_string()])]
     pub solo_cell_filter: Vec<String>,
@@ -1615,15 +1620,15 @@ impl Parameters {
                     ));
                 }
             }
-            // Gene / GeneFull / SJ / Velocyto are implemented.
+            // Gene / GeneFull / SJ / Velocyto / Transcript3p are implemented.
             for f in &params.solo_features {
-                if !matches!(f.as_str(), "SJ" | "Velocyto")
+                if !matches!(f.as_str(), "SJ" | "Velocyto" | "Transcript3p")
                     && f.parse::<crate::solo::SoloFeature>().is_err()
                 {
                     return Err(command.error(
                         ErrorKind::InvalidValue,
                         format!(
-                            "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto"
+                            "unsupported --soloFeatures '{f}'; supported: Gene, GeneFull, SJ, Velocyto, Transcript3p"
                         ),
                     ));
                 }
@@ -1710,6 +1715,16 @@ impl Parameters {
                     ));
                 }
             }
+            // Transcript3p quantifies per cluster, so it needs the clustering.
+            if params.solo_features.iter().any(|f| f == "Transcript3p")
+                && params.solo_cluster_cb_file.is_none()
+            {
+                return Err(command.error(
+                    ErrorKind::MissingRequiredArgument,
+                    "--soloFeatures Transcript3p requires --soloClusterCBfile: the EM runs \
+                     per cluster of cells, since one cell has too few UMIs to resolve isoforms",
+                ));
+            }
             // Validate --clipAdapterType.
             if !matches!(
                 params.clip_adapter_type.as_str(),

src/solo/count.rs

@@ -1323,6 +1323,47 @@ pub fn write_gene_matrix(
         }
     }
 
+    // Transcript3p: rows are transcripts, columns are clusters rather than
+    // cells, since the isoform EM needs more UMIs than one cell provides.
+    if let (Some(acc), Some(tx)) = (&ctx.transcript3p, &ctx.transcriptome) {
+        let dir = params.output_path(&format!("{solo_dir}Transcript3p/raw/"));
+        std::fs::create_dir_all(&dir).map_err(|e| Error::io(e, &dir))?;
+
+        let cluster_cb = match &params.solo_cluster_cb_file {
+            Some(path) => {
+                let text = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?;
+                crate::solo::transcript3p::load_cluster_cb(&text, |cb| {
+                    ctx.whitelist.index_of_barcode(cb.as_bytes())
+                })
+            }
+            // Validation requires the file, so this is unreachable in practice;
+            // an empty map quantifies nothing rather than inventing a cluster.
+            None => std::collections::BTreeMap::new(),
+        };
+
+        let acc = acc.lock().unwrap();
+        let out = crate::solo::transcript3p::quantify(&acc, tx, &cluster_cb);
+        for (name, body) in [
+            (matrix_name.as_str(), &out.matrix),
+            (features_name.as_str(), &out.features),
+            (
+                "transcriptEndDistanceDistribution.txt",
+                &out.distance_distribution,
+            ),
+        ] {
+            let path = dir.join(name);
+            std::fs::write(&path, body).map_err(|e| Error::io(e, &path))?;
+        }
+        log::info!(
+            "STARsolo: wrote Transcript3p/raw ({} transcripts × {} clusters)",
+            tx.n_transcripts(),
+            cluster_cb
+                .values()
+                .collect::<std::collections::BTreeSet<_>>()
+                .len(),
+        );
+    }
+
     // SJ (splice-junction) feature: rows are the SJ.out.tab junctions.
     if ctx.sj_enabled
         && let Some(sjs) = sj_stats

src/solo/mod.rs

@@ -12,6 +12,7 @@
 pub mod count;
 pub mod gene;
 pub mod smartseq;
+pub mod transcript3p;
 pub mod whitelist;
 
 pub use count::{UmiDedup, UmiFiltering, write_gene_matrix};
@@ -564,6 +565,11 @@ pub struct SoloContext {
     /// `--soloMultiMappers` includes a non-`Unique` method → capture gene-
     /// ambiguous reads for distribution into `UniqueAndMult-*.mtx`.
     pub want_multi: bool,
+    /// `--soloFeatures Transcript3p`: the transcriptome to assign reads to, and
+    /// the accumulated per-read records. Both `None`/absent unless asked for,
+    /// since building the transcriptome costs a GTF pass.
+    pub transcriptome: Option<crate::quant::transcriptome::TranscriptomeIndex>,
+    pub transcript3p: Option<Mutex<crate::solo::transcript3p::Transcript3pAcc>>,
 }
 
 /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped
@@ -677,6 +683,7 @@ impl SoloContext {
         let feature_reads = features.iter().map(|_| AtomicU64::new(0)).collect();
         let sj_enabled = params.solo_features.iter().any(|f| f == "SJ");
         let velocyto_enabled = params.solo_features.iter().any(|f| f == "Velocyto");
+        let transcript3p = params.solo_features.iter().any(|f| f == "Transcript3p");
         let want_multi = params.solo_multi_mappers.iter().any(|m| m != "Unique");
 
         Ok(Self {
@@ -695,6 +702,20 @@ impl SoloContext {
             velocyto_enabled,
             velocyto_records: Mutex::new(Vec::new()),
             want_multi,
+            transcriptome: transcript3p
+                .then(|| {
+                    crate::quant::transcriptome::TranscriptomeIndex::from_gtf_exons_configured(
+                        &exons,
+                        genome,
+                        &params.sjdb_gtf_tag_exon_parent_transcript,
+                        &params.sjdb_gtf_tag_exon_parent_gene,
+                        &params.sjdb_gtf_tag_exon_parent_gene_name,
+                        &params.sjdb_gtf_tag_exon_parent_gene_type,
+                    )
+                })
+                .transpose()?,
+            transcript3p: transcript3p
+                .then(|| Mutex::new(crate::solo::transcript3p::Transcript3pAcc::new())),
         })
     }
 
@@ -871,6 +892,25 @@ impl SoloContext {
                 fo
             })
             .collect();
+
+        // Transcript3p: every transcript this read is concordant with, and how
+        // far its 3' end sits from each transcript's. Uniquely-mapped reads
+        // only — a read at several genomic loci says nothing about isoforms.
+        if let (Some(acc), Some(tx), Some(cb)) =
+            (&self.transcript3p, &self.transcriptome, cb_resolved)
+            && n_loci == 1
+            && let Some(align) = cdna_transcripts.first()
+        {
+            let hits = crate::solo::transcript3p::concordant_transcripts(
+                align,
+                tx,
+                align.read_length() as u32,
+            );
+            if !hits.is_empty() {
+                acc.lock().unwrap().add(cb, umi, hits);
+            }
+        }
+
         out
     }