Skip to content

Write gene_id_to_name.tsv when indexing a probe set CSV #202

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
an-altosian wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Write gene_id_to_name.tsv when indexing a probe set CSV #202

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 39 additions & 11 deletions src/simpleaf_commands/indexing.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -175,6 +176,17 @@ struct ProbeRow {
probe_id: String,
included: Option<Included>,
region: Option<ProbeRegion>,
// 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<String>,
}

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 {
Expand Down Expand Up @@ -275,7 +287,6 @@ fn parse_csv_record(
has_region: bool,
seq_id_hs: &mut HashSet<String>,
ref_seq_writer: &mut BufWriter<File>,
// id_to_name_writer: &mut BufWriter<File>,
t2g_writer: &mut BufWriter<File>,
) -> anyhow::Result<()> {
if !include {
Expand All @@ -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(())
Expand Down Expand Up @@ -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 {
Expand All @@ -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<String, String> = BTreeMap::new();

match csv_reader {
CsvReader::Feature(mut rdr) => {
Expand All @@ -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,
)?;
}
Expand All @@ -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(),
Expand All @@ -501,20 +522,27 @@ 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,
)?;
}
}
}

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)?;
Expand Down
73 changes: 57 additions & 16 deletions src/utils/probe_utils.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,41 @@ fn get_required_idx(headers: &csv::StringRecord, name: &str) -> anyhow::Result<u
.with_context(|| format!("probe CSV is missing required column `{}`", name))
}

/// Record a `gene_id -> 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<String, String>,
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<String, String>,
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
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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));
Expand Down Expand Up @@ -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");
Expand Down
Loading