Skip to content

Commit ac97387

Browse files
committed
feat(index): write gene_id_to_name.tsv for --probe-csv indexes
When `simpleaf index` builds from a 10x probe set CSV that carries a gene symbol column (`gene_name` or `gene_symbol`), it now emits a `gene_id_to_name.tsv` mapping into the reference and the built index. Previously this mapping was only produced on the auto-build path inside `multiplex-quant`/`quant` (via probe_utils), so a prebuilt probe index passed with `--index` carried no gene names. The GTF/roers path already wrote the file; this brings the probe-csv path to parity, letting downstream `quant` surface gene names for prebuilt probe indexes. The mapping covers every gene in the probe set independent of the `included` flag (it is a complete gene annotation, not matrix contents); conflicting names for the same gene_id are rejected. Adds a unit test.
1 parent 52ccc73 commit ac97387

1 file changed

Lines changed: 78 additions & 12 deletions

File tree

src/simpleaf_commands/indexing.rs

Lines changed: 78 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ use anyhow::{Context, anyhow, bail};
77
use roers;
88
use serde::Deserialize;
99
use serde_json::json;
10-
use std::collections::HashSet;
10+
use std::collections::{BTreeMap, HashSet};
1111
use std::fs::File;
1212
use std::io::{BufWriter, Write};
1313
use std::path::{Path, PathBuf};
@@ -90,7 +90,26 @@ fn write_index_log_stage(
9090

9191
#[cfg(test)]
9292
mod tests {
93-
use super::derive_kmer_and_minimizer;
93+
use super::{derive_kmer_and_minimizer, insert_gene_name};
94+
use std::collections::BTreeMap;
95+
96+
#[test]
97+
fn insert_gene_name_dedups_and_detects_conflicts() {
98+
let mut m = BTreeMap::new();
99+
insert_gene_name(&mut m, "G1", "GeneOne").expect("first insert ok");
100+
insert_gene_name(&mut m, "G1", "GeneOne").expect("identical re-insert ok");
101+
insert_gene_name(&mut m, "G2", "GeneTwo").expect("distinct gene ok");
102+
assert_eq!(m.len(), 2);
103+
assert_eq!(m.get("G1").map(String::as_str), Some("GeneOne"));
104+
105+
let err = insert_gene_name(&mut m, "G1", "Different")
106+
.expect_err("conflicting name for same gene_id should error");
107+
assert!(
108+
format!("{:#}", err).contains("inconsistent gene annotations"),
109+
"unexpected error: {:#}",
110+
err
111+
);
112+
}
94113

95114
#[test]
96115
fn derive_kmer_and_minimizer_fails_for_short_reference() {
@@ -175,6 +194,17 @@ struct ProbeRow {
175194
probe_id: String,
176195
included: Option<Included>,
177196
region: Option<ProbeRegion>,
197+
// optional gene symbol column (10x probe set v2 CSVs include `gene_name`;
198+
// some panels name it `gene_symbol`). Used to emit a gene_id -> name map.
199+
#[serde(default, alias = "gene_symbol")]
200+
gene_name: Option<String>,
201+
}
202+
203+
impl ProbeRow {
204+
/// The gene symbol/name for this probe's gene, if the CSV provided one.
205+
fn gene_name(&self) -> Option<&str> {
206+
self.gene_name.as_deref()
207+
}
178208
}
179209

180210
impl CsvRow<'_> for ProbeRow {
@@ -265,6 +295,26 @@ impl std::fmt::Display for ProbeRegion {
265295
}
266296
}
267297

298+
/// Record a `gene_id -> gene_name` association, erroring if a different name was
299+
/// already seen for the same `gene_id` (an internally inconsistent probe set).
300+
fn insert_gene_name(
301+
map: &mut BTreeMap<String, String>,
302+
gene_id: &str,
303+
gene_name: &str,
304+
) -> anyhow::Result<()> {
305+
if let Some(prev) = map.insert(gene_id.to_string(), gene_name.to_string())
306+
&& prev != gene_name
307+
{
308+
bail!(
309+
"probe CSV contains inconsistent gene annotations for `{}`: saw both `{}` and `{}`.",
310+
gene_id,
311+
prev,
312+
gene_name
313+
);
314+
}
315+
Ok(())
316+
}
317+
268318
#[allow(clippy::too_many_arguments)]
269319
fn parse_csv_record(
270320
ref_id: &str,
@@ -275,7 +325,6 @@ fn parse_csv_record(
275325
has_region: bool,
276326
seq_id_hs: &mut HashSet<String>,
277327
ref_seq_writer: &mut BufWriter<File>,
278-
// id_to_name_writer: &mut BufWriter<File>,
279328
t2g_writer: &mut BufWriter<File>,
280329
) -> anyhow::Result<()> {
281330
if !include {
@@ -301,9 +350,6 @@ fn parse_csv_record(
301350
writeln!(t2g_writer, "{}\t{}", seq_id, ref_id)?;
302351
};
303352

304-
// insert into gene id to name
305-
// writeln!(id_to_name_writer, "{}\t{}", ref_id, ref_id)?;
306-
307353
// insert into ref seq
308354
writeln!(ref_seq_writer, ">{}\n{}", seq_id, sequence)?;
309355
Ok(())
@@ -453,7 +499,6 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu
453499

454500
// define file names
455501
let ref_seq_path = outref.join("ref.fa");
456-
// let id_to_name_path = outref.join("gene_id_to_name.tsv");
457502
let t2g_path = if has_region {
458503
outref.join("t2g_3col.tsv")
459504
} else {
@@ -462,9 +507,11 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu
462507

463508
// define buffer writers
464509
let mut ref_seq_writer = BufWriter::new(File::create(&ref_seq_path)?);
465-
// let mut id_to_name_writer = BufWriter::new(File::create(&id_to_name_path)?);
466510
let mut t2g_writer = BufWriter::new(File::create(&t2g_path)?);
467511
let mut msl = u32::MAX;
512+
// collected gene_id -> gene_name for probe CSVs that carry a gene symbol column;
513+
// written out as gene_id_to_name.tsv so downstream `quant` can surface gene names.
514+
let mut gene_id_to_name_map: BTreeMap<String, String> = BTreeMap::new();
468515

469516
match csv_reader {
470517
CsvReader::Feature(mut rdr) => {
@@ -482,7 +529,6 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu
482529
has_region,
483530
&mut seq_id_hs,
484531
&mut ref_seq_writer,
485-
// &mut id_to_name_writer,
486532
&mut t2g_writer,
487533
)?;
488534
}
@@ -492,6 +538,15 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu
492538
for row in rdr.deserialize() {
493539
let record: ProbeRow = row?;
494540

541+
// record gene_id -> gene_name for every probe that carries a name,
542+
// independent of the `included` flag: the mapping is a complete gene
543+
// annotation, written whenever the probe set provides gene symbols.
544+
if let Some(gene_name) =
545+
record.gene_name().map(str::trim).filter(|s| !s.is_empty())
546+
{
547+
insert_gene_name(&mut gene_id_to_name_map, record.ref_id(), gene_name)?;
548+
}
549+
495550
parse_csv_record(
496551
record.ref_id(),
497552
record.seq_id(),
@@ -501,20 +556,31 @@ pub fn build_ref_and_index(af_home_path: &Path, opts: IndexOpts) -> anyhow::Resu
501556
has_region,
502557
&mut seq_id_hs,
503558
&mut ref_seq_writer,
504-
// &mut id_to_name_writer,
505559
&mut t2g_writer,
506560
)?;
507561
}
508562
}
509563
}
510564

511565
index_info["t2g_file"] = json!(&t2g_path);
512-
// index_info["gene_id_to_name"] = json!(&id_to_name_path);
566+
567+
// If the (probe) CSV carried gene symbols, emit a gene_id -> gene_name map.
568+
// This parallels the GTF/roers path (above) and the multiplex-quant auto-build
569+
// path, so a prebuilt probe index also lets `quant` surface gene names.
570+
if !gene_id_to_name_map.is_empty() {
571+
let id_to_name_path = outref.join("gene_id_to_name.tsv");
572+
let mut id_to_name_writer = BufWriter::new(File::create(&id_to_name_path)?);
573+
for (gene_id, gene_name) in &gene_id_to_name_map {
574+
writeln!(id_to_name_writer, "{}\t{}", gene_id, gene_name)?;
575+
}
576+
id_to_name_writer.flush()?;
577+
index_info["gene_id_to_name"] = json!(&id_to_name_path);
578+
gene_id_to_name = Some(id_to_name_path);
579+
}
513580

514581
min_seq_len = Some(msl);
515582
reference_sequence = Some(ref_seq_path);
516583
t2g = Some(t2g_path);
517-
// _gene_id_to_name = Some(id_to_name_path);
518584
}
519585

520586
io::write_json_pretty(&info_file, &index_info)?;

0 commit comments

Comments
 (0)