diff --git a/src/simpleaf_commands/indexing.rs b/src/simpleaf_commands/indexing.rs index 34b5e8c..fdd0e99 100644 --- a/src/simpleaf_commands/indexing.rs +++ b/src/simpleaf_commands/indexing.rs @@ -2,12 +2,13 @@ use crate::core::{context, exec, io, runtime}; use crate::utils::af_utils::create_dir_if_absent; use crate::utils::prog_utils; use crate::utils::prog_utils::ReqProgs; +use crate::utils::probe_utils; use anyhow::{Context, anyhow, bail}; use roers; use serde::Deserialize; use serde_json::json; -use std::collections::HashSet; +use std::collections::{BTreeMap, HashSet}; use std::fs::File; use std::io::{BufWriter, Write}; use std::path::{Path, PathBuf}; @@ -175,6 +176,17 @@ struct ProbeRow { probe_id: String, included: Option, region: Option, + // optional gene symbol column (10x probe set v2 CSVs include `gene_name`; + // some panels name it `gene_symbol`). Used to emit a gene_id -> name map. + #[serde(default, alias = "gene_symbol")] + gene_name: Option, +} + +impl ProbeRow { + /// The gene symbol/name for this probe's gene, if the CSV provided one. + fn gene_name(&self) -> Option<&str> { + self.gene_name.as_deref() + } } impl CsvRow<'_> for ProbeRow { @@ -275,7 +287,6 @@ fn parse_csv_record( has_region: bool, seq_id_hs: &mut HashSet, ref_seq_writer: &mut BufWriter, - // id_to_name_writer: &mut BufWriter, t2g_writer: &mut BufWriter, ) -> anyhow::Result<()> { if !include { @@ -301,9 +312,6 @@ fn parse_csv_record( writeln!(t2g_writer, "{}\t{}", seq_id, ref_id)?; }; - // insert into gene id to name - // writeln!(id_to_name_writer, "{}\t{}", ref_id, ref_id)?; - // insert into ref seq writeln!(ref_seq_writer, ">{}\n{}", seq_id, sequence)?; Ok(()) @@ -453,7 +461,6 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu // define file names let ref_seq_path = outref.join("ref.fa"); - // let id_to_name_path = outref.join("gene_id_to_name.tsv"); let t2g_path = if has_region { outref.join("t2g_3col.tsv") } else { @@ -462,9 +469,11 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu // define buffer writers let mut ref_seq_writer = BufWriter::new(File::create(&ref_seq_path)?); - // let mut id_to_name_writer = BufWriter::new(File::create(&id_to_name_path)?); let mut t2g_writer = BufWriter::new(File::create(&t2g_path)?); let mut msl = u32::MAX; + // collected gene_id -> gene_name for probe CSVs that carry a gene symbol column; + // written out as gene_id_to_name.tsv so downstream `quant` can surface gene names. + let mut gene_id_to_name_map: BTreeMap = BTreeMap::new(); match csv_reader { CsvReader::Feature(mut rdr) => { @@ -482,7 +491,6 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu has_region, &mut seq_id_hs, &mut ref_seq_writer, - // &mut id_to_name_writer, &mut t2g_writer, )?; } @@ -492,6 +500,19 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu for row in rdr.deserialize() { let record: ProbeRow = row?; + // record gene_id -> gene_name for every probe that carries a name, + // independent of the `included` flag: the mapping is a complete gene + // annotation, written whenever the probe set provides gene symbols. + if let Some(gene_name) = + record.gene_name().map(str::trim).filter(|s| !s.is_empty()) + { + probe_utils::insert_gene_name( + &mut gene_id_to_name_map, + record.ref_id(), + gene_name, + )?; + } + parse_csv_record( record.ref_id(), record.seq_id(), @@ -501,7 +522,6 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu has_region, &mut seq_id_hs, &mut ref_seq_writer, - // &mut id_to_name_writer, &mut t2g_writer, )?; } @@ -509,12 +529,20 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu } index_info["t2g_file"] = json!(&t2g_path); - // index_info["gene_id_to_name"] = json!(&id_to_name_path); + + // If the (probe) CSV carried gene symbols, emit a gene_id -> gene_name map. + // This parallels the GTF/roers path (above) and the multiplex-quant auto-build + // path, so a prebuilt probe index also lets `quant` surface gene names. + if !gene_id_to_name_map.is_empty() { + let id_to_name_path = outref.join("gene_id_to_name.tsv"); + probe_utils::write_gene_id_to_name(&gene_id_to_name_map, &id_to_name_path)?; + index_info["gene_id_to_name"] = json!(&id_to_name_path); + gene_id_to_name = Some(id_to_name_path); + } min_seq_len = Some(msl); reference_sequence = Some(ref_seq_path); t2g = Some(t2g_path); - // _gene_id_to_name = Some(id_to_name_path); } io::write_json_pretty(&info_file, &index_info)?; diff --git a/src/utils/probe_utils.rs b/src/utils/probe_utils.rs index 5627181..459a922 100644 --- a/src/utils/probe_utils.rs +++ b/src/utils/probe_utils.rs @@ -52,6 +52,41 @@ fn get_required_idx(headers: &csv::StringRecord, name: &str) -> anyhow::Result gene_name` association into `map`, erroring if a *different* +/// name was already recorded for the same `gene_id` (an internally inconsistent probe +/// set). Shared by the probe-set conversion (auto-build) and `simpleaf index --probe-csv`. +pub fn insert_gene_name( + map: &mut BTreeMap, + gene_id: &str, + gene_name: &str, +) -> anyhow::Result<()> { + if let Some(prev) = map.insert(gene_id.to_string(), gene_name.to_string()) + && prev != gene_name + { + bail!( + "probe CSV contains inconsistent gene annotations for `{}`: saw both `{}` and `{}`.", + gene_id, + prev, + gene_name, + ); + } + Ok(()) +} + +/// Write a `gene_id -> gene_name` map as a 2-column TSV (rows sorted by `gene_id`, +/// since the map is a `BTreeMap`). Shared by both probe-index build paths. +pub fn write_gene_id_to_name( + map: &BTreeMap, + path: &Path, +) -> anyhow::Result<()> { + let mut writer = BufWriter::new(std::fs::File::create(path)?); + for (gene_id, gene_name) in map { + writeln!(writer, "{}\t{}", gene_id, gene_name)?; + } + writer.flush()?; + Ok(()) +} + /// Convert a 10x probe set CSV file to a FASTA file suitable for indexing. /// /// Also generates a collapsed gene-level transcript-to-gene (t2g) map and, when @@ -140,16 +175,7 @@ pub fn convert_probe_csv_to_reference_files( && let Some(gene_name) = record.get(gene_name_i).map(str::trim) && !gene_name.is_empty() { - if let Some(prev) = gene_id_to_name.insert(gene_id.to_string(), gene_name.to_string()) - && prev != gene_name - { - bail!( - "probe CSV contains inconsistent gene annotations for `{}`: saw both `{}` and `{}`.", - gene_id, - prev, - gene_name, - ); - } + insert_gene_name(&mut gene_id_to_name, gene_id, gene_name)?; } if let Some(region_i) = region_idx { @@ -179,11 +205,7 @@ Expected `spliced` or `unspliced`.", writer.flush()?; } if let Some(ref path) = gene_id_to_name_path { - let mut writer = BufWriter::new(std::fs::File::create(path)?); - for (gene_id, gene_name) in &gene_id_to_name { - writeln!(writer, "{}\t{}", gene_id, gene_name)?; - } - writer.flush()?; + write_gene_id_to_name(&gene_id_to_name, path)?; } metadata.insert("num_probes".to_string(), json!(num_probes)); @@ -312,11 +334,30 @@ Provide a probe CSV with a `region` column (`spliced` / `unspliced`), or a pre-b mod tests { use super::{ ProbeT2gMode, collapse_t2g_to_gene, convert_probe_csv_to_reference_files, ensure_t2g_mode, - t2g_has_usa_mapping, + insert_gene_name, t2g_has_usa_mapping, }; + use std::collections::BTreeMap; use std::fs; use tempfile::tempdir; + #[test] + fn insert_gene_name_dedups_and_detects_conflicts() { + let mut m = BTreeMap::new(); + insert_gene_name(&mut m, "G1", "GeneOne").expect("first insert ok"); + insert_gene_name(&mut m, "G1", "GeneOne").expect("identical re-insert ok"); + insert_gene_name(&mut m, "G2", "GeneTwo").expect("distinct gene ok"); + assert_eq!(m.len(), 2); + assert_eq!(m.get("G1").map(String::as_str), Some("GeneOne")); + + let err = insert_gene_name(&mut m, "G1", "Different") + .expect_err("conflicting name for same gene_id should error"); + assert!( + format!("{:#}", err).contains("inconsistent gene annotations"), + "unexpected error: {:#}", + err + ); + } + #[test] fn convert_probe_csv_writes_gene_and_usa_t2g_files() { let td = tempdir().expect("failed to create tempdir");