diff --git a/cli/src/commands/genotype_to_vcf.rs b/cli/src/commands/genotype_to_vcf.rs index 05ba1e5..0b520ef 100644 --- a/cli/src/commands/genotype_to_vcf.rs +++ b/cli/src/commands/genotype_to_vcf.rs @@ -1,8 +1,7 @@ use std::collections::HashMap; use std::fs::File; -use std::io::{BufRead, BufWriter, Write}; +use std::io::{BufRead, BufReader, BufWriter, Write}; use std::path::{Path, PathBuf}; -use std::process::Command; use anyhow::{bail, Context, Result}; use rayon::prelude::*; @@ -50,6 +49,10 @@ pub fn run_genotype_to_vcf(args: GenotypeToVcfArgs) -> Result<()> { }; let reference_map = load_reference_map(&store)?; let position_map = load_position_map(&store)?; + let coordinate_map = match args.coordinate_map.as_deref() { + Some(path) => load_coordinate_map(path)?, + None => HashMap::new(), + }; let cache_path = args.cache.clone().unwrap_or_else(default_cache_path); let cache = if cache_path.exists() { @@ -81,6 +84,7 @@ pub fn run_genotype_to_vcf(args: GenotypeToVcfArgs) -> Result<()> { &output_paths, &reference_map, &position_map, + &coordinate_map, &cache, &sample_name, &missing_log_path, @@ -128,6 +132,7 @@ fn convert_file( output_paths: &OutputPaths, reference_map: &HashMap, position_map: &HashMap<(String, i64), ReferenceVariant>, + coordinate_map: &HashMap, cache: &RsidCache, sample_name: &str, missing_log_path: &Path, @@ -179,7 +184,16 @@ fn convert_file( let results = pool.install(|| { rows.par_iter() - .map(|row| convert_row(reference_map, position_map, cache, row, include_metrics)) + .map(|row| { + convert_row( + reference_map, + position_map, + coordinate_map, + cache, + row, + include_metrics, + ) + }) .collect::>() }); @@ -195,6 +209,9 @@ fn convert_file( #[derive(Clone)] struct Emit { key: String, + chrom_rank: u16, + chrom: String, + pos: i64, rsid: String, gt: String, line: String, @@ -213,6 +230,9 @@ fn convert_file( let (chrom, pos) = result.locus.clone().unwrap_or_default(); emits.push(Emit { key: format!("{}\t{}\t{}", chrom, pos, result.rsid), + chrom_rank: contig_rank(&chrom), + chrom, + pos, rsid: result.rsid, gt: result.gt.unwrap_or_else(|| "./.".to_string()), line, @@ -232,6 +252,11 @@ fn convert_file( } groups.entry(e.key.clone()).or_default().push(e); } + order.sort_by(|left, right| { + let a = &groups[left][0]; + let b = &groups[right][0]; + (a.chrom_rank, &a.chrom, a.pos, &a.rsid).cmp(&(b.chrom_rank, &b.chrom, b.pos, &b.rsid)) + }); let mut written = 0usize; for key in order { @@ -296,6 +321,7 @@ fn collect_row( fn convert_row( reference_map: &HashMap, position_map: &HashMap<(String, i64), ReferenceVariant>, + coordinate_map: &HashMap, cache: &RsidCache, row: &crate::genotype_reader::GenotypeRow, include_metrics: bool, @@ -332,11 +358,20 @@ fn convert_row( } let rsid_label = rsid_label.as_str(); + let rsid_int = normalize_rsid(rsid_label) + .trim_start_matches("rs") + .parse::() + .ok(); + let coordinate_override = rsid_int.and_then(|rsid| coordinate_map.get(&rsid)); + let mut unresolved = false; let mut reference = "N".to_string(); let mut alternates: Vec = Vec::new(); - if let Some(resolved_ref) = resolved.as_ref() { + if let Some(override_ref) = coordinate_override { + reference = override_ref.reference.clone(); + alternates = override_ref.alternates.clone(); + } else if let Some(resolved_ref) = resolved.as_ref() { reference = resolved_ref.reference.clone(); alternates = resolved_ref.alternates.clone(); } else { @@ -400,12 +435,15 @@ fn convert_row( result.stats.written_rows += 1; result.rsid = rsid_label.to_string(); - result.locus = Some((row.chrom.clone(), row.pos)); + let (emit_chrom, emit_pos) = coordinate_override + .map(|entry| (entry.chrom.clone(), entry.pos)) + .unwrap_or_else(|| (row.chrom.clone(), row.pos)); + result.locus = Some((emit_chrom.clone(), emit_pos)); result.gt = Some(gt_call.clone()); result.line = Some(format!( "{}\t{}\t{}\t{}\t{}\t.\tPASS\t{}\t{}\t{}\n", - row.chrom, - row.pos, + emit_chrom, + emit_pos, rsid_label, reference, alt_field, @@ -438,6 +476,16 @@ fn force_nocall(line: &str) -> String { out } +fn contig_rank(chrom: &str) -> u16 { + let normalized = chrom.trim().trim_start_matches("chr"); + match normalized { + "X" => 23, + "Y" => 24, + "M" | "MT" => 25, + _ => normalized.parse::().unwrap_or(10_000), + } +} + fn resolve_reference( rsid: &str, reference_map: &HashMap, @@ -468,6 +516,67 @@ fn resolve_reference( None } +fn load_coordinate_map(path: &Path) -> Result> { + let reader = BufReader::new( + File::open(path).with_context(|| format!("Open coordinate map {:?}", path))?, + ); + let mut map = HashMap::new(); + for (line_no, line) in reader.lines().enumerate() { + let line = line.with_context(|| format!("Read coordinate map {:?}", path))?; + let trimmed = line.trim(); + if trimmed.is_empty() || trimmed.starts_with('#') { + continue; + } + let fields: Vec<&str> = trimmed.split('\t').collect(); + if fields.len() < 5 { + bail!( + "Coordinate map {:?}: line {} has fewer than 5 tab-separated fields", + path, + line_no + 1 + ); + } + if line_no == 0 && fields[0].eq_ignore_ascii_case("rsid") { + continue; + } + let rsid = normalize_rsid(fields[0]) + .trim_start_matches("rs") + .parse::() + .with_context(|| { + format!( + "Coordinate map {:?}: invalid rsid '{}' on line {}", + path, + fields[0], + line_no + 1 + ) + })?; + let chrom = fields[1].trim().trim_start_matches("chr").to_string(); + let pos = fields[2].trim().parse::().with_context(|| { + format!( + "Coordinate map {:?}: invalid position '{}' on line {}", + path, + fields[2], + line_no + 1 + ) + })?; + let reference = normalize_sequence(fields[3]); + let alternates = parse_alternates(fields[4]); + if reference.is_empty() || alternates.is_empty() { + bail!( + "Coordinate map {:?}: empty ref/alt on line {}", + path, + line_no + 1 + ); + } + map.entry(rsid).or_insert(CoordinateOverride { + chrom, + pos, + reference, + alternates, + }); + } + Ok(map) +} + fn write_vcf_header( writer: &mut BufWriter, sample_name: &str, @@ -492,6 +601,10 @@ fn write_vcf_header( writer, "##INFO=" )?; + writeln!( + writer, + "##INFO=" + )?; writeln!( writer, "##FILTER=" @@ -604,14 +717,17 @@ fn default_sample_name(input: &Path) -> String { } fn gzip_output(path: &Path) -> Result<()> { - let status = Command::new("gzip") - .arg("-f") - .arg(path) - .status() - .with_context(|| "Failed to execute gzip")?; - if !status.success() { - bail!("gzip exited with status {}", status); - } + let mut final_name = path.as_os_str().to_os_string(); + final_name.push(".gz"); + let final_path = PathBuf::from(final_name); + let input = File::open(path).with_context(|| format!("Open {:?}", path))?; + let mut reader = BufReader::new(input); + let output = File::create(&final_path).with_context(|| format!("Create {:?}", final_path))?; + let mut writer = noodles::bgzf::io::Writer::new(output); + std::io::copy(&mut reader, &mut writer) + .with_context(|| format!("BGZF-compress {:?} to {:?}", path, final_path))?; + writer.finish().context("Finish BGZF stream")?; + std::fs::remove_file(path).with_context(|| format!("Remove temporary VCF {:?}", path))?; Ok(()) } @@ -695,6 +811,14 @@ struct ResolvedReference { alternates: Vec, } +#[derive(Clone, Debug)] +struct CoordinateOverride { + chrom: String, + pos: i64, + reference: String, + alternates: Vec, +} + fn build_genotype_call( genotype: Option<&str>, reference: &str, diff --git a/cli/src/main.rs b/cli/src/main.rs index e238197..1f177c6 100644 --- a/cli/src/main.rs +++ b/cli/src/main.rs @@ -201,6 +201,10 @@ pub struct GenotypeToVcfArgs { /// Optional rsid cache file for fallback lookups. #[arg(long)] pub cache: Option, + /// Optional TSV map overriding emitted VCF loci/ref/alt by rsID. + /// Expected columns: rsid, chrom, pos, ref, alt. + #[arg(long)] + pub coordinate_map: Option, } #[derive(Args, Clone)]