From d0896783192e89599d3dc2da6415042b643df92a Mon Sep 17 00:00:00 2001 From: Madhava Jay Date: Tue, 19 May 2026 13:39:36 +1000 Subject: [PATCH 1/2] adding support for illumina files and assembly detection --- rust/bioscript-formats/src/genotype.rs | 5 +- .../src/genotype/delimited.rs | 2 + .../src/genotype/delimited/gsgt.rs | 294 ++++++++++++++ .../src/genotype/delimited/scan.rs | 224 +++++++++-- .../bioscript-formats/src/genotype/loaders.rs | 152 ++++++- rust/bioscript-formats/src/inspect.rs | 18 +- .../src/inspect/assembly_anchors.rs | 371 ++++++++++++++++++ .../src/inspect/heuristics.rs | 68 +++- rust/bioscript-formats/src/inspect/sex.rs | 75 +++- .../tests/file_formats/delimited.rs | 228 +++++++++++ .../carigenetics_ddna_concordance_sample.txt | 8 + .../fixtures/carigenetics_gsgt_sample.txt | 31 ++ rust/bioscript-formats/tests/inspect.rs | 131 +++++++ 13 files changed, 1560 insertions(+), 47 deletions(-) create mode 100644 rust/bioscript-formats/src/genotype/delimited/gsgt.rs create mode 100644 rust/bioscript-formats/src/inspect/assembly_anchors.rs create mode 100644 rust/bioscript-formats/tests/fixtures/carigenetics_ddna_concordance_sample.txt create mode 100644 rust/bioscript-formats/tests/fixtures/carigenetics_gsgt_sample.txt diff --git a/rust/bioscript-formats/src/genotype.rs b/rust/bioscript-formats/src/genotype.rs index ee6f420..ccad4a4 100644 --- a/rust/bioscript-formats/src/genotype.rs +++ b/rust/bioscript-formats/src/genotype.rs @@ -30,9 +30,10 @@ pub use cram_backend::{ observe_cram_deletion_with_reader, observe_cram_indel_with_reader, observe_cram_snp_with_reader, }; pub(crate) use delimited::{ - COMMENT_PREFIXES, DelimitedColumnIndexes, Delimiter, detect_delimiter, parse_streaming_row, + COMMENT_PREFIXES, DelimitedColumnIndexes, Delimiter, GsgtParser, detect_delimiter, + lines_look_like_gsgt, parse_streaming_row, }; -use delimited::{RowParser, scan_delimited_variants}; +use delimited::{RowParser, is_no_call as gsgt_is_no_call, scan_delimited_variants}; use io::{ detect_source_format, is_bgzf_path, looks_like_vcf_lines, read_lines_from_reader, select_zip_entry, diff --git a/rust/bioscript-formats/src/genotype/delimited.rs b/rust/bioscript-formats/src/genotype/delimited.rs index 02b1da3..947da2d 100644 --- a/rust/bioscript-formats/src/genotype/delimited.rs +++ b/rust/bioscript-formats/src/genotype/delimited.rs @@ -4,8 +4,10 @@ use bioscript_core::RuntimeError; use super::normalize_genotype; +mod gsgt; mod scan; +pub(crate) use gsgt::{GsgtParser, is_no_call, lines_look_like_gsgt}; pub(crate) use scan::scan_delimited_variants; pub(crate) const COMMENT_PREFIXES: [&str; 2] = ["#", "//"]; diff --git a/rust/bioscript-formats/src/genotype/delimited/gsgt.rs b/rust/bioscript-formats/src/genotype/delimited/gsgt.rs new file mode 100644 index 0000000..66a4069 --- /dev/null +++ b/rust/bioscript-formats/src/genotype/delimited/gsgt.rs @@ -0,0 +1,294 @@ +//! Illumina `GenomeStudio` GSGT Final Report ("Carigenetics") support. +//! +//! Maps GSGT rows into the existing `ParsedDelimitedRow` so the standard +//! rsid-first / `(chrom,pos)` locus matching engine handles the rest. See +//! `illumina.md` for the format spec and the quirks this module encodes. + +use bioscript_core::RuntimeError; + +use super::{ParsedDelimitedRow, normalize_name, sanitize_evidence_line}; +use crate::genotype::normalize_genotype; + +/// True when the sampled lines are a GSGT Final Report: the first non-empty +/// line is `[Header]` (case-insensitive). +pub(crate) fn lines_look_like_gsgt(lines: &[String]) -> bool { + for line in lines { + let trimmed = line.trim(); + if trimmed.is_empty() { + continue; + } + return trimmed.eq_ignore_ascii_case("[header]"); + } + false +} + +/// Extract the first `rs\d+` substring from a GSGT `SNP Name` probe id. +/// +/// `BOT-rs1135675` -> `rs1135675`, `rs111647200_ilmndup1` -> `rs111647200`, +/// `1:103380393` -> `None` (the embedded coordinate is *not* an rsid). +pub(crate) fn extract_rsid(snp_name: &str) -> Option { + let bytes = snp_name.as_bytes(); + let mut i = 0; + while i + 2 < bytes.len() { + let is_rs = (bytes[i] == b'r' || bytes[i] == b'R') + && (bytes[i + 1] == b's' || bytes[i + 1] == b'S'); + if is_rs && bytes[i + 2].is_ascii_digit() { + let mut j = i + 2; + while j < bytes.len() && bytes[j].is_ascii_digit() { + j += 1; + } + return Some(format!("rs{}", &snp_name[i + 2..j])); + } + i += 1; + } + None +} + +#[derive(Debug, Clone, Copy)] +struct GsgtColumns { + snp_name: usize, + chr: usize, + position: usize, + allele1_plus: usize, + allele2_plus: usize, +} + +#[derive(Debug, Clone, Copy, PartialEq)] +enum Phase { + /// Inside the `[Header]` metadata block, before `[Data]`. + Meta, + /// `[Data]` seen; the next non-empty line is the column header. + ExpectHeader, + /// Column header parsed; remaining lines are data rows. + Body, +} + +/// Streaming GSGT row parser. Feed it every line of the file in order; it +/// skips the metadata block, learns the columns, then yields one +/// `ParsedDelimitedRow` per usable data row. +#[derive(Debug)] +pub(crate) struct GsgtParser { + phase: Phase, + cols: Option, +} + +impl GsgtParser { + pub(crate) fn new() -> Self { + Self { + phase: Phase::Meta, + cols: None, + } + } + + pub(crate) fn consume( + &mut self, + line: &str, + ) -> Result, RuntimeError> { + let trimmed = line.trim(); + match self.phase { + Phase::Meta => { + if trimmed.eq_ignore_ascii_case("[data]") { + self.phase = Phase::ExpectHeader; + } + Ok(None) + } + Phase::ExpectHeader => { + if trimmed.is_empty() { + return Ok(None); + } + self.cols = Some(resolve_columns(&split_tab(line))?); + self.phase = Phase::Body; + Ok(None) + } + Phase::Body => { + if trimmed.is_empty() { + return Ok(None); + } + let cols = self.cols.expect("columns resolved before body"); + let fields = split_tab(line); + Ok(extract_row(line, &fields, &cols)) + } + } + } +} + +fn split_tab(line: &str) -> Vec<&str> { + line.trim_end_matches(['\n', '\r']) + .split('\t') + .map(str::trim) + .collect() +} + +fn resolve_columns(header: &[&str]) -> Result { + let find = |want: &str| header.iter().position(|h| normalize_name(h) == want); + let snp_name = find("snpname"); + let chr = find("chr").or_else(|| find("chromosome")).or_else(|| find("chrom")); + let position = find("position").or_else(|| find("pos")); + let allele1_plus = find("allele1plus"); + let allele2_plus = find("allele2plus"); + match (snp_name, chr, position, allele1_plus, allele2_plus) { + (Some(snp_name), Some(chr), Some(position), Some(allele1_plus), Some(allele2_plus)) => { + Ok(GsgtColumns { + snp_name, + chr, + position, + allele1_plus, + allele2_plus, + }) + } + _ => Err(RuntimeError::Unsupported( + "GSGT Final Report missing required columns (SNP Name, Chr, Position, Allele1 - Plus, Allele2 - Plus)" + .to_owned(), + )), + } +} + +fn extract_row( + line: &str, + fields: &[&str], + cols: &GsgtColumns, +) -> Option { + let snp_name = fields.get(cols.snp_name).copied().unwrap_or(""); + let rsid = extract_rsid(snp_name); + + let chrom_raw = fields.get(cols.chr).copied().unwrap_or("").trim(); + // Skip unplaced markers (Chr/Position == 0). + if chrom_raw.is_empty() || chrom_raw == "0" { + return None; + } + let position = fields + .get(cols.position) + .and_then(|value| value.trim().parse::().ok()) + .filter(|pos| *pos != 0)?; + + if rsid.is_none() { + // No rsid and we still have a locus; keep it for the locus fallback. + // (chrom/position validated above.) + } + + let a1 = fields.get(cols.allele1_plus).copied().unwrap_or("").trim(); + let a2 = fields.get(cols.allele2_plus).copied().unwrap_or("").trim(); + // GSGT no-call is a single `-` per Plus allele column. + let genotype = if (a1 == "-" || a1.is_empty()) && (a2 == "-" || a2.is_empty()) { + "--".to_owned() + } else { + normalize_genotype(&format!("{a1}{a2}")) + }; + + Some(ParsedDelimitedRow { + rsid, + chrom: Some(chrom_raw.to_owned()), + position: Some(position), + genotype, + raw_line: sanitize_evidence_line(line), + }) +} + +/// Whether a normalized genotype string represents a no-call. +pub(crate) fn is_no_call(genotype: &str) -> bool { + genotype == "--" || genotype.is_empty() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn lines_look_like_gsgt_detects_header_marker() { + assert!(lines_look_like_gsgt(&[ + String::new(), + "[Header]".to_owned() + ])); + assert!(lines_look_like_gsgt(&["[HEADER]".to_owned()])); + assert!(!lines_look_like_gsgt(&[ + "# Dynamic DNA".to_owned(), + "[Header]".to_owned() + ])); + assert!(!lines_look_like_gsgt(&[ + "rsid\tchromosome\tposition\tgenotype".to_owned() + ])); + } + + #[test] + fn extract_rsid_recovers_rs_from_every_probe_form() { + assert_eq!(extract_rsid("rs11466023").as_deref(), Some("rs11466023")); + assert_eq!(extract_rsid("BOT-rs1135675").as_deref(), Some("rs1135675")); + assert_eq!( + extract_rsid("rs111647200_ilmndup1").as_deref(), + Some("rs111647200") + ); + assert_eq!(extract_rsid("GSA-rs61660502").as_deref(), Some("rs61660502")); + assert_eq!(extract_rsid("seq-rs786202193").as_deref(), Some("rs786202193")); + // No rs id: chr:pos, CNV, MNV, vendor/clinical. + assert_eq!(extract_rsid("1:103380393"), None); + assert_eq!(extract_rsid("1:159174749-C-T"), None); + assert_eq!(extract_rsid("1:110228436_CNV_GSTM1"), None); + assert_eq!(extract_rsid("1:45332290_MNV"), None); + assert_eq!(extract_rsid("DICER1-chr14-95596479"), None); + assert_eq!(extract_rsid("GALC:c.2002A>C"), None); + } + + fn body_parser() -> GsgtParser { + let mut p = GsgtParser::new(); + for line in [ + "[Header]", + "GSGT Version\t2.0.5", + "[Data]", + "Sample ID\tSample Name\tSNP Name\tSNP\tChr\tPosition\tAllele1 - Top\tAllele2 - Top\tAllele1 - Plus\tAllele2 - Plus\tPlus/Minus Strand", + ] { + assert!(p.consume(line).unwrap().is_none()); + } + p + } + + fn row(p: &mut GsgtParser, snp_name: &str, snp: &str, chr: &str, pos: &str, a1: &str, a2: &str) -> Option { + let line = format!("S1\t\t{snp_name}\t{snp}\t{chr}\t{pos}\tX\tX\t{a1}\t{a2}\t+"); + p.consume(&line).unwrap() + } + + #[test] + fn parses_plus_columns_and_snp_name_not_design_column() { + let mut p = body_parser(); + let r = row(&mut p, "rs9000001", "[A/G]", "1", "1000", "G", "G").unwrap(); + assert_eq!(r.rsid.as_deref(), Some("rs9000001")); + assert_eq!(r.chrom.as_deref(), Some("1")); + assert_eq!(r.position, Some(1000)); + assert_eq!(r.genotype, "GG"); + } + + #[test] + fn design_snp_column_is_never_used_as_rsid() { + let mut p = body_parser(); + // SNP Name has no rs id; SNP design column is [A/G]. rsid must be None. + let r = row(&mut p, "1:7000_CNV_GENEX", "[A/G]", "1", "7000", "A", "A").unwrap(); + assert_eq!(r.rsid, None); + assert_eq!(r.position, Some(7000)); + assert_eq!(r.genotype, "AA"); + } + + #[test] + fn no_call_dash_maps_to_double_dash_and_indels_pass_through() { + let mut p = body_parser(); + let nc = row(&mut p, "rs9000011", "[A/G]", "4", "9000", "-", "-").unwrap(); + assert_eq!(nc.genotype, "--"); + let indel = row(&mut p, "rs9000012", "[I/D]", "4", "9100", "I", "D").unwrap(); + assert!(indel.genotype != "--"); + } + + #[test] + fn skips_unplaced_chr_or_position_zero() { + let mut p = body_parser(); + assert!(row(&mut p, "rs9000009", "[A/G]", "0", "9000", "A", "G").is_none()); + assert!(row(&mut p, "rs9000010", "[A/G]", "4", "0", "A", "G").is_none()); + } + + #[test] + fn empty_rsid_kept_when_locus_present() { + let mut p = body_parser(); + let r = row(&mut p, "2:6000", "[C/T]", "2", "6000", "C", "T").unwrap(); + assert_eq!(r.rsid, None); + assert_eq!(r.chrom.as_deref(), Some("2")); + assert_eq!(r.position, Some(6000)); + assert_eq!(r.genotype, "CT"); + } +} diff --git a/rust/bioscript-formats/src/genotype/delimited/scan.rs b/rust/bioscript-formats/src/genotype/delimited/scan.rs index a9248dd..039ad36 100644 --- a/rust/bioscript-formats/src/genotype/delimited/scan.rs +++ b/rust/bioscript-formats/src/genotype/delimited/scan.rs @@ -6,7 +6,7 @@ use std::{ use zip::ZipArchive; -use crate::inspect::detect_assembly; +use crate::inspect::{AssemblyAnchorScorer, detect_assembly}; use bioscript_core::{RuntimeError, VariantObservation, VariantSpec}; use super::{ @@ -14,8 +14,158 @@ use super::{ GenotypeSourceFormat, backends::delimited_locus_for_assembly, describe_query, types::DelimitedBackend, variant_sort_key, }, - DelimitedColumnIndexes, detect_delimiter, parse_streaming_row, + DelimitedColumnIndexes, GsgtParser, detect_delimiter, is_no_call as gsgt_is_no_call, + lines_look_like_gsgt, parse_streaming_row, }; +use bioscript_core::Assembly; + +/// Apply a candidate genotype to a result slot. +/// +/// Non-GSGT keeps the historical first-match-wins behaviour. GSGT merges +/// replicate probes in a single pass (spec §5): a real call replaces an +/// earlier no-call; two disagreeing real calls collapse to a no-call; equal +/// calls and no-call-after-call are ignored. +#[allow(clippy::too_many_arguments)] +fn apply_match( + slot: &mut VariantObservation, + unresolved: &mut usize, + gsgt: bool, + matched_rsid: Option, + genotype: &str, + evidence_head: &str, + raw_line: &str, + backend_name: &str, +) { + let make = |g: &str, rsid: Option| VariantObservation { + backend: backend_name.to_owned(), + matched_rsid: rsid, + genotype: Some(g.to_owned()), + evidence: vec![evidence_head.to_owned(), format!("source line: {raw_line}")], + ..VariantObservation::default() + }; + + if slot.genotype.is_none() { + *slot = make(genotype, matched_rsid); + *unresolved = unresolved.saturating_sub(1); + return; + } + if !gsgt { + return; + } + let existing = slot.genotype.clone().unwrap_or_default(); + let prev_rsid = slot.matched_rsid.clone(); + let existing_nc = gsgt_is_no_call(&existing); + let new_nc = gsgt_is_no_call(genotype); + if existing_nc && !new_nc { + *slot = make(genotype, matched_rsid); + } else if !existing_nc && !new_nc && existing != genotype { + *slot = make("--", prev_rsid); + } +} + +/// One bounded pass over the file feeding the rsID/locus anchor vote, used +/// only when no build metadata was found. Mirrors the main scan's open +/// logic; reuses the GSGT / delimited row parsers. +fn prescan_assembly_anchors( + backend: &DelimitedBackend, + is_gsgt: bool, +) -> Result, RuntimeError> { + fn vote(mut reader: R, is_gsgt: bool) -> Result, RuntimeError> { + let mut scorer = AssemblyAnchorScorer::new(); + let mut probe = Vec::new(); + let mut buf = String::new(); + for _ in 0..8 { + buf.clear(); + if reader.read_line(&mut buf).unwrap_or(0) == 0 { + break; + } + probe.push(buf.trim_end_matches(['\n', '\r']).to_owned()); + } + let delimiter = detect_delimiter(&probe); + let mut ci: Option = None; + let mut ch: Option> = None; + let mut gsgt = GsgtParser::new(); + let feed = |line: &str, + scorer: &mut AssemblyAnchorScorer, + gsgt: &mut GsgtParser, + ci: &mut Option, + ch: &mut Option>| + -> Result<(), RuntimeError> { + let row = if is_gsgt { + gsgt.consume(line)? + } else { + parse_streaming_row(line, delimiter, ci, ch)? + }; + if let Some(row) = row + && let (Some(c), Some(p)) = (row.chrom.as_ref(), row.position) + { + scorer.observe(row.rsid.as_deref().unwrap_or(""), c, p, &row.genotype); + } + Ok(()) + }; + for line in &probe { + feed(line, &mut scorer, &mut gsgt, &mut ci, &mut ch)?; + } + loop { + buf.clear(); + if reader + .read_line(&mut buf) + .map_err(|err| RuntimeError::Io(format!("anchor prescan read failed: {err}")))? + == 0 + { + break; + } + feed( + buf.trim_end_matches(['\n', '\r']), + &mut scorer, + &mut gsgt, + &mut ci, + &mut ch, + )?; + } + Ok(scorer.decide()) + } + + match backend.format { + GenotypeSourceFormat::Text => { + let file = File::open(&backend.path).map_err(|err| { + RuntimeError::Io(format!( + "failed to open genotype file {}: {err}", + backend.path.display() + )) + })?; + vote(BufReader::new(file), is_gsgt) + } + GenotypeSourceFormat::Zip => { + let entry_name = backend.zip_entry_name.as_ref().ok_or_else(|| { + RuntimeError::Unsupported(format!( + "zip backend missing selected entry for {}", + backend.path.display() + )) + })?; + let file = File::open(&backend.path).map_err(|err| { + RuntimeError::Io(format!( + "failed to open genotype zip {}: {err}", + backend.path.display() + )) + })?; + let mut archive = ZipArchive::new(file).map_err(|err| { + RuntimeError::Io(format!( + "failed to read genotype zip {}: {err}", + backend.path.display() + )) + })?; + let entry = archive.by_name(entry_name).map_err(|err| { + RuntimeError::Io(format!( + "failed to open genotype entry {entry_name} in {}: {err}", + backend.path.display() + )) + })?; + vote(BufReader::new(entry), is_gsgt) + } + _ => Ok(None), + } +} pub(crate) fn scan_delimited_variants( backend: &DelimitedBackend, @@ -67,6 +217,7 @@ pub(crate) fn scan_delimited_variants( } let delimiter = detect_delimiter(&probe_lines); + let is_gsgt = lines_look_like_gsgt(&probe_lines); if detected_assembly.is_none() { let mut label = backend.path.to_string_lossy().to_ascii_lowercase(); if let Some(entry_name) = backend.zip_entry_name.as_ref() { @@ -74,6 +225,12 @@ pub(crate) fn scan_delimited_variants( label.push_str(&entry_name.to_ascii_lowercase()); } detected_assembly = detect_assembly(&label, &probe_lines); + // No declared build (e.g. a GSGT Final Report): resolve it from + // the rsID/locus anchor vote over the whole file instead of + // assuming. Only runs when metadata gave us nothing. + if detected_assembly.is_none() { + detected_assembly = prescan_assembly_anchors(backend, is_gsgt)?; + } if let Some(assembly) = detected_assembly { for (idx, variant) in &indexed { if let Some(locus) = delimited_locus_for_assembly(variant, Some(assembly)) { @@ -95,35 +252,39 @@ pub(crate) fn scan_delimited_variants( } let mut column_indexes: Option = None; let mut comment_header: Option> = None; + let mut gsgt_parser = GsgtParser::new(); + let backend_name = backend.backend_name(); let mut process_line = |line: &str| -> Result { - let Some(row) = + let row = if is_gsgt { + gsgt_parser.consume(line)? + } else { parse_streaming_row(line, delimiter, &mut column_indexes, &mut comment_header)? - else { - return Ok(unresolved == 0); + }; + let Some(row) = row else { + // GSGT must scan the whole file so later replicate probes can + // be merged; non-GSGT may stop once everything resolved. + return Ok(!is_gsgt && unresolved == 0); }; if let Some(rsid) = row.rsid.as_ref() && let Some(target_indexes) = rsid_targets.get(rsid) { for &target_idx in target_indexes { - if results[target_idx].genotype.is_none() { - results[target_idx] = VariantObservation { - backend: backend.backend_name().to_owned(), - matched_rsid: Some(rsid.clone()), - genotype: Some(row.genotype.clone()), - evidence: vec![ - format!("resolved by rsid {rsid}"), - format!("source line: {}", row.raw_line), - ], - ..VariantObservation::default() - }; - unresolved = unresolved.saturating_sub(1); - } + apply_match( + &mut results[target_idx], + &mut unresolved, + is_gsgt, + Some(rsid.clone()), + &row.genotype, + &format!("resolved by rsid {rsid}"), + &row.raw_line, + backend_name, + ); } } - if unresolved == 0 { + if !is_gsgt && unresolved == 0 { return Ok(true); } @@ -134,23 +295,20 @@ pub(crate) fn scan_delimited_variants( ); if let Some(target_indexes) = coord_targets.get(&key) { for &target_idx in target_indexes { - if results[target_idx].genotype.is_none() { - results[target_idx] = VariantObservation { - backend: backend.backend_name().to_owned(), - matched_rsid: row.rsid.clone(), - genotype: Some(row.genotype.clone()), - evidence: vec![ - format!("resolved by locus {}:{}", chrom, position), - format!("source line: {}", row.raw_line), - ], - ..VariantObservation::default() - }; - unresolved = unresolved.saturating_sub(1); - } + apply_match( + &mut results[target_idx], + &mut unresolved, + is_gsgt, + row.rsid.clone(), + &row.genotype, + &format!("resolved by locus {chrom}:{position}"), + &row.raw_line, + backend_name, + ); } } } - Ok(unresolved == 0) + Ok(!is_gsgt && unresolved == 0) }; for line in &probe_lines { diff --git a/rust/bioscript-formats/src/genotype/loaders.rs b/rust/bioscript-formats/src/genotype/loaders.rs index feed1db..bb71808 100644 --- a/rust/bioscript-formats/src/genotype/loaders.rs +++ b/rust/bioscript-formats/src/genotype/loaders.rs @@ -2,11 +2,12 @@ use std::{collections::HashMap, io::BufRead}; use bioscript_core::{Assembly, RuntimeError}; -use crate::inspect::detect_assembly; +use crate::inspect::{AssemblyAnchorScorer, detect_assembly}; use super::{ - COMMENT_PREFIXES, GenotypeSourceFormat, GenotypeStore, QueryBackend, RowParser, RsidMapBackend, - delimited::sanitize_evidence_line, detect_delimiter, vcf_tokens::genotype_from_vcf_gt, + COMMENT_PREFIXES, GenotypeSourceFormat, GenotypeStore, GsgtParser, QueryBackend, RowParser, + RsidMapBackend, delimited::sanitize_evidence_line, detect_delimiter, gsgt_is_no_call, + lines_look_like_gsgt, vcf_tokens::genotype_from_vcf_gt, }; pub(crate) fn from_vcf_reader( @@ -66,8 +67,15 @@ pub(crate) fn from_delimited_reader( } } + if lines_look_like_gsgt(&prelude) { + return from_gsgt_reader(format, &prelude, reader, &mut buf, label); + } + let mut parser = RowParser::new(delimiter.unwrap_or(super::Delimiter::Tab)); - let assembly = detect_assembly(&label.to_ascii_lowercase(), &prelude); + // Build from declared metadata first; only fall back to the rsID/locus + // anchor vote when the file gives us no idea (per spec / no blind guess). + let meta_assembly = detect_assembly(&label.to_ascii_lowercase(), &prelude); + let mut scorer = meta_assembly.is_none().then(AssemblyAnchorScorer::new); let mut values = HashMap::new(); let mut locus_values = HashMap::new(); let mut source_lines = HashMap::new(); @@ -78,6 +86,7 @@ pub(crate) fn from_delimited_reader( &mut values, &mut locus_values, &mut source_lines, + scorer.as_mut(), )?; } loop { @@ -94,9 +103,118 @@ pub(crate) fn from_delimited_reader( &mut values, &mut locus_values, &mut source_lines, + scorer.as_mut(), )?; } + let assembly = meta_assembly.or_else(|| scorer.as_ref().and_then(AssemblyAnchorScorer::decide)); + Ok(from_rsid_map( + format, + values, + locus_values, + assembly, + source_lines, + )) +} + +struct GsgtMergeGroup { + chrom: String, + position: i64, + rsid: Option, + genotypes: Vec, + source_line: String, +} + +/// GSGT loader: parse the Final Report, merge replicate probes keyed by +/// `(chrom, pos, rsid)` (spec §5), then build the in-memory rsid/locus maps. +fn from_gsgt_reader( + format: GenotypeSourceFormat, + prelude: &[String], + mut reader: R, + buf: &mut String, + label: &str, +) -> Result { + let mut parser = GsgtParser::new(); + // Insertion-ordered groups so output is deterministic. + let mut order: Vec<(String, i64, String)> = Vec::new(); + let mut groups: HashMap<(String, i64, String), GsgtMergeGroup> = HashMap::new(); + + let ingest = |row: super::delimited::ParsedDelimitedRow, + order: &mut Vec<(String, i64, String)>, + groups: &mut HashMap<(String, i64, String), GsgtMergeGroup>| { + let (Some(chrom), Some(position)) = (row.chrom.clone(), row.position) else { + return; + }; + let chrom_norm = chrom.trim_start_matches("chr").to_ascii_lowercase(); + let rsid_key = row.rsid.clone().unwrap_or_default(); + let key = (chrom_norm, position, rsid_key); + let entry = groups.entry(key.clone()).or_insert_with(|| { + order.push(key.clone()); + GsgtMergeGroup { + chrom, + position, + rsid: row.rsid.clone(), + genotypes: Vec::new(), + source_line: row.raw_line.clone(), + } + }); + entry.genotypes.push(row.genotype); + }; + + for line in prelude { + if let Some(row) = parser.consume(line)? { + ingest(row, &mut order, &mut groups); + } + } + loop { + buf.clear(); + let bytes = reader + .read_line(buf) + .map_err(|err| RuntimeError::Io(format!("failed to read {label}: {err}")))?; + if bytes == 0 { + break; + } + if let Some(row) = parser.consume(buf.trim_end_matches(['\n', '\r']))? { + ingest(row, &mut order, &mut groups); + } + } + + // GSGT Final Reports declare no build; try real metadata first (none + // here), then resolve the assembly from the rsID/locus anchor vote over + // the merged calls instead of blindly assuming GRCh38. + let meta_assembly = detect_assembly(&label.to_ascii_lowercase(), prelude); + let mut scorer = meta_assembly.is_none().then(AssemblyAnchorScorer::new); + + let mut values = HashMap::new(); + let mut locus_values = HashMap::new(); + let mut source_lines = HashMap::new(); + for key in order { + let group = groups.remove(&key).expect("group recorded in order exists"); + let genotype = merge_genotypes(&group.genotypes); + let chrom_norm = group.chrom.trim_start_matches("chr").to_ascii_lowercase(); + if let Some(scorer) = scorer.as_mut() { + scorer.observe( + group.rsid.as_deref().unwrap_or(""), + &group.chrom, + group.position, + &genotype, + ); + } + locus_values.insert( + (chrom_norm, group.position), + ( + genotype.clone(), + group.rsid.clone(), + group.source_line.clone(), + ), + ); + if let Some(rsid) = group.rsid { + values.insert(rsid.clone(), genotype); + source_lines.insert(rsid, group.source_line); + } + } + + let assembly = meta_assembly.or_else(|| scorer.as_ref().and_then(AssemblyAnchorScorer::decide)); Ok(from_rsid_map( format, values, @@ -106,6 +224,21 @@ pub(crate) fn from_delimited_reader( )) } +/// Collapse a replicate-probe group to one genotype: drop no-calls; if none +/// remain emit a no-call; if the remaining calls all agree emit that; +/// otherwise emit a no-call (genuine disagreement — never auto-pick). +fn merge_genotypes(genotypes: &[String]) -> String { + let mut calls = genotypes.iter().filter(|g| !gsgt_is_no_call(g)); + let Some(first) = calls.next() else { + return "--".to_owned(); + }; + if calls.all(|g| g == first) { + first.clone() + } else { + "--".to_owned() + } +} + pub(crate) fn from_vcf_lines(lines: Vec) -> Result { let mut values = HashMap::new(); for line in lines { @@ -158,9 +291,20 @@ fn consume_delimited_line( values: &mut HashMap, locus_values: &mut HashMap<(String, i64), (String, Option, String)>, source_lines: &mut HashMap, + scorer: Option<&mut AssemblyAnchorScorer>, ) -> Result<(), RuntimeError> { if let Some(row) = parser.consume_record(line)? { let source_line = sanitize_evidence_line(line); + if let (Some(scorer), Some(chrom), Some(position)) = + (scorer, row.chrom.as_ref(), row.position) + { + scorer.observe( + row.rsid.as_deref().unwrap_or(""), + chrom, + position, + &row.genotype, + ); + } if let (Some(chrom), Some(position)) = (row.chrom.as_ref(), row.position) { locus_values.insert( ( diff --git a/rust/bioscript-formats/src/inspect.rs b/rust/bioscript-formats/src/inspect.rs index b94fe05..a687312 100644 --- a/rust/bioscript-formats/src/inspect.rs +++ b/rust/bioscript-formats/src/inspect.rs @@ -31,11 +31,13 @@ impl StubDuration { use bioscript_core::{Assembly, RuntimeError}; +mod assembly_anchors; mod heuristics; mod io; mod render; mod sex; +pub(crate) use assembly_anchors::{AssemblyAnchorScorer, assembly_from_text_bytes}; pub(crate) use heuristics::*; pub(crate) use io::*; #[cfg(test)] @@ -179,7 +181,13 @@ pub fn inspect_bytes( }; let inspection_context = inspect_context_name(&lower, options); let source = detect_source(&inspection_context, &sample_lines, detected_kind); - let assembly = detect_assembly(&inspection_context, &sample_lines); + // Declared metadata first; fall back to the rsID/locus anchor vote over + // the buffer only when the file declares no build (e.g. a GSGT report). + let assembly = detect_assembly(&inspection_context, &sample_lines).or_else(|| { + matches!(detected_kind, DetectedKind::GenotypeText | DetectedKind::Vcf) + .then(|| assembly_from_text_bytes(bytes)) + .flatten() + }); let phased = (detected_kind == DetectedKind::Vcf) .then(|| detect_vcf_phasing(&sample_lines)) .flatten(); @@ -280,7 +288,13 @@ pub fn inspect_file(path: &Path, options: &InspectOptions) -> Result Option { + match base.to_ascii_uppercase() { + b'A' => Some(b'T'), + b'T' => Some(b'A'), + b'C' => Some(b'G'), + b'G' => Some(b'C'), + _ => None, + } +} + +/// Normalize a chromosome label to bare form (`chr7` -> `7`, `23` -> `X`). +fn norm_chrom(raw: &str) -> String { + let c = raw.trim(); + let c = c.strip_prefix("chr").or_else(|| c.strip_prefix("CHR")).unwrap_or(c); + match c.to_ascii_uppercase().as_str() { + "23" => "X".to_owned(), + "24" => "Y".to_owned(), + "25" => "XY".to_owned(), + "26" | "M" => "MT".to_owned(), + _ => c.to_ascii_uppercase(), + } +} + +/// Observed ACGT letters in a genotype must be a subset of the anchor's +/// allele set or its complement (strand-agnostic). Empty/no-call passes. +fn allele_ok(genotype: &str, anchor_alleles: &str) -> bool { + let observed: Vec = genotype + .bytes() + .map(|b| b.to_ascii_uppercase()) + .filter(|b| matches!(b, b'A' | b'C' | b'G' | b'T')) + .collect(); + if observed.is_empty() { + return true; + } + let fwd: Vec = anchor_alleles.bytes().map(|b| b.to_ascii_uppercase()).collect(); + let rev: Vec = fwd.iter().filter_map(|b| complement(*b)).collect(); + let subset = |set: &[u8]| observed.iter().all(|o| set.contains(o)); + subset(&fwd) || subset(&rev) +} + +/// Accumulates anchor votes from a streamed genotype/VCF file. +#[derive(Debug, Default)] +pub(crate) struct AssemblyAnchorScorer { + rs: [u32; 3], + loc: [u32; 3], +} + +impl AssemblyAnchorScorer { + pub(crate) fn new() -> Self { + Self::default() + } + + /// Observe one parsed row. `rsid` may be empty (locus-only vote). + pub(crate) fn observe(&mut self, rsid: &str, chrom: &str, pos: i64, genotype: &str) { + let chrom = norm_chrom(chrom); + for (a_rsid, a_chrom, p36, p37, p38, alleles) in ANCHORS { + if chrom != *a_chrom { + continue; + } + let build = if pos == *p36 { + 0 + } else if pos == *p37 { + 1 + } else if pos == *p38 { + 2 + } else { + continue; + }; + if rsid == *a_rsid { + self.rs[build] += 1; + } else if allele_ok(genotype, alleles) { + self.loc[build] += 1; + } + return; + } + } + + fn scores(&self) -> [u32; 3] { + [ + self.rs[0] * RSID_WEIGHT + self.loc[0] * LOCUS_WEIGHT, + self.rs[1] * RSID_WEIGHT + self.loc[1] * LOCUS_WEIGHT, + self.rs[2] * RSID_WEIGHT + self.loc[2] * LOCUS_WEIGHT, + ] + } + + /// Resolve to a supported assembly, or `None` when there is too little + /// signal, the vote is split, or a build-36 file dominates (unsupported). + pub(crate) fn decide(&self) -> Option { + let anchors: u32 = self.rs.iter().chain(self.loc.iter()).sum(); + if anchors < MIN_ANCHORS { + return None; + } + let s = self.scores(); + let total: u32 = s.iter().sum(); + if total == 0 { + return None; + } + let (best_idx, &best) = s + .iter() + .enumerate() + .max_by_key(|(_, v)| **v) + .expect("three scores"); + if f64::from(best) / f64::from(total) < MIN_DOMINANCE { + return None; + } + match best_idx { + 1 => Some(Assembly::Grch37), + 2 => Some(Assembly::Grch38), + // best_idx == 0 -> GRCh36, which we do not model: honest unknown. + _ => None, + } + } +} + +/// Vote an assembly from an in-memory text buffer (genotype text or VCF), +/// vendor-agnostic. Handles GSGT `[Header]/[Data]`, VCF, and flat +/// rsid/chrom/pos/genotype delimited rows. Returns `None` unless the anchor +/// vote is confident — used only as the metadata-absent fallback. +#[allow(clippy::items_after_statements, clippy::many_single_char_names)] +pub(crate) fn assembly_from_text_bytes(bytes: &[u8]) -> Option { + enum Mode { + Auto, + GsgtMeta, + GsgtHdr, + GsgtBody(usize, usize, usize, usize, usize), + Vcf, + } + let text = String::from_utf8_lossy(bytes); + let mut scorer = AssemblyAnchorScorer::new(); + let mut mode = Mode::Auto; + for raw in text.lines() { + let line = raw.trim_end_matches('\r'); + let t = line.trim(); + if t.is_empty() { + continue; + } + match mode { + Mode::Auto if t.eq_ignore_ascii_case("[header]") => { + mode = Mode::GsgtMeta; + } + Mode::GsgtMeta => { + if t.eq_ignore_ascii_case("[data]") { + mode = Mode::GsgtHdr; + } + } + Mode::GsgtHdr => { + let norm = |s: &str| { + s.trim() + .to_ascii_lowercase() + .replace([' ', '-', '_'], "") + }; + let h: Vec = line.split('\t').map(norm).collect(); + let idx = |w: &str| h.iter().position(|x| x == w); + match ( + idx("snpname"), + idx("chr"), + idx("position"), + idx("allele1plus"), + idx("allele2plus"), + ) { + (Some(s), Some(c), Some(p), Some(a1), Some(a2)) => { + mode = Mode::GsgtBody(s, c, p, a1, a2); + } + _ => return None, + } + } + Mode::GsgtBody(s, c, p, a1, a2) => { + let f: Vec<&str> = line.split('\t').collect(); + if f.len() <= a2.max(p).max(c).max(s) { + continue; + } + if let Ok(pos) = f[p].trim().parse::() { + let rsid = extract_rs(f[s]); + scorer.observe( + rsid.as_deref().unwrap_or(""), + f[c].trim(), + pos, + &format!("{}{}", f[a1].trim(), f[a2].trim()), + ); + } + } + Mode::Vcf => { + let f: Vec<&str> = t.split('\t').collect(); + if f.len() >= 5 + && let Ok(pos) = f[1].trim().parse::() + { + let rsid = if f[2].starts_with("rs") { f[2] } else { "" }; + scorer.observe(rsid, f[0], pos, &format!("{}{}", f[3], f[4])); + } + } + Mode::Auto => { + if t.starts_with("##") { + continue; + } + if t.starts_with("#CHROM\t") { + mode = Mode::Vcf; + continue; + } + if t.starts_with('#') || t.starts_with("//") { + continue; + } + let f: Vec<&str> = if line.contains('\t') { + line.split('\t').collect() + } else { + line.split(',').collect() + }; + if f.len() < 3 { + continue; + } + let c0 = f[0].trim().trim_matches('"').to_ascii_lowercase(); + if matches!(c0.as_str(), "rsid" | "rs id" | "snp" | "name" | "snpname") { + continue; // column header row + } + if let Ok(pos) = f[2].trim().trim_matches('"').parse::() { + let id = f[0].trim().trim_matches('"'); + let rsid = if id.starts_with("rs") || id.starts_with("RS") { + id + } else { + "" + }; + let gt: String = f[3..].iter().map(|x| x.trim().trim_matches('"')).collect(); + scorer.observe(rsid, f[1].trim().trim_matches('"'), pos, >); + } + } + } + } + scorer.decide() +} + +fn extract_rs(name: &str) -> Option { + let b = name.as_bytes(); + let mut i = 0; + while i + 2 < b.len() { + if (b[i] | 32) == b'r' && (b[i + 1] | 32) == b's' && b[i + 2].is_ascii_digit() { + let mut j = i + 2; + while j < b.len() && b[j].is_ascii_digit() { + j += 1; + } + return Some(format!("rs{}", &name[i + 2..j])); + } + i += 1; + } + None +} + +#[cfg(test)] +mod tests { + use super::*; + + fn feed(scorer: &mut AssemblyAnchorScorer, build: usize, with_rsid: bool) { + for (rsid, chrom, p36, p37, p38, alleles) in ANCHORS { + let pos = [*p36, *p37, *p38][build]; + let id = if with_rsid { *rsid } else { "" }; + scorer.observe(id, chrom, pos, alleles); + } + } + + #[test] + fn rsid_votes_resolve_grch37_and_grch38() { + let mut g37 = AssemblyAnchorScorer::new(); + feed(&mut g37, 1, true); + assert_eq!(g37.decide(), Some(Assembly::Grch37)); + + let mut g38 = AssemblyAnchorScorer::new(); + feed(&mut g38, 2, true); + assert_eq!(g38.decide(), Some(Assembly::Grch38)); + } + + #[test] + fn locus_only_no_rsid_still_resolves() { + // VCF-style: no usable rsIDs, alleles corroborate the position. + let mut g38 = AssemblyAnchorScorer::new(); + feed(&mut g38, 2, false); + assert_eq!(g38.decide(), Some(Assembly::Grch38)); + } + + #[test] + fn build36_file_is_honest_unknown_not_misread() { + let mut g36 = AssemblyAnchorScorer::new(); + feed(&mut g36, 0, true); + assert_eq!(g36.decide(), None); + } + + #[test] + fn too_few_anchors_is_unknown() { + let mut s = AssemblyAnchorScorer::new(); + s.observe("rs429358", "19", 45_411_941, "CT"); + s.observe("rs7412", "19", 45_412_079, "C"); + assert_eq!(s.decide(), None); + } + + #[test] + fn split_vote_is_unknown_not_a_guess() { + let mut s = AssemblyAnchorScorer::new(); + for (i, (rsid, chrom, p36, p37, p38, al)) in ANCHORS.iter().enumerate() { + let pos = if i % 2 == 0 { *p37 } else { *p38 }; + let _ = (p36,); + s.observe(rsid, chrom, pos, al); + } + assert_eq!(s.decide(), None); + } + + #[test] + fn chrom_prefix_and_strand_complement_are_handled() { + let mut s = AssemblyAnchorScorer::new(); + // "chr19" prefix + complement-strand alleles for rs429358 (CT -> GA). + for (rsid, chrom, _p36, p37, _p38, _al) in ANCHORS { + s.observe(rsid, &format!("chr{chrom}"), *p37, "GA"); + } + assert_eq!(s.decide(), Some(Assembly::Grch37)); + } +} diff --git a/rust/bioscript-formats/src/inspect/heuristics.rs b/rust/bioscript-formats/src/inspect/heuristics.rs index 2dd9490..37d4ab2 100644 --- a/rust/bioscript-formats/src/inspect/heuristics.rs +++ b/rust/bioscript-formats/src/inspect/heuristics.rs @@ -12,6 +12,11 @@ pub(crate) fn looks_like_vcf_lines(lines: &[String]) -> bool { } pub(crate) fn looks_like_genotype_text(lines: &[String]) -> bool { + // Illumina GSGT Final Report: a [Header]/[Data] block whose data rows + // don't match the flat rsid/chrom/pos/genotype shape but are genotypes. + if detect_gsgt_final_report(lines).is_some() { + return true; + } let mut checked = 0usize; let mut valid = 0usize; for line in lines { @@ -121,6 +126,7 @@ pub(crate) fn detect_source( let mut vendor = None; let mut platform_version = None; let mut confidence = DetectionConfidence::Unknown; + let gsgt_final_report = detect_gsgt_final_report(sample_lines); if normalized.contains("genes for good") || normalized.contains("geneforgood") { vendor = Some("Genes for Good".to_owned()); @@ -189,10 +195,22 @@ pub(crate) fn detect_source( DetectionConfidence::StrongHeuristic }; evidence.push("sequencing.com path/header text".to_owned()); - } else if normalized.contains("carigenetics") || normalized.contains("cari genetics") { + } else if normalized.contains("carigenetics") + || normalized.contains("cari genetics") + || gsgt_final_report.is_some() + { vendor = Some("CariGenetics".to_owned()); confidence = DetectionConfidence::StrongHeuristic; - evidence.push("CariGenetics path/header text".to_owned()); + if normalized.contains("carigenetics") || normalized.contains("cari genetics") { + evidence.push("CariGenetics path/header text".to_owned()); + } + if let Some(report) = gsgt_final_report { + evidence.push("GSGT [Header]/[Data] block".to_owned()); + if report.platform_version.is_some() { + evidence.push("GSGT Version / .bpm manifest".to_owned()); + } + platform_version = report.platform_version; + } } vendor.map(|vendor| SourceMetadata { @@ -203,6 +221,50 @@ pub(crate) fn detect_source( }) } +/// A detected Illumina `GenomeStudio` GSGT Final Report. +struct GsgtReport { + /// `.bpm` manifest stem or `GSGT Version`, when present. + platform_version: Option, +} + +/// Detect an Illumina `GenomeStudio` GSGT Final Report by its +/// `[Header]` … `[Data]` block. +fn detect_gsgt_final_report(sample_lines: &[String]) -> Option { + let mut saw_header = false; + let mut saw_data = false; + let mut version: Option = None; + for line in sample_lines { + let trimmed = line.trim(); + if trimmed.eq_ignore_ascii_case("[header]") { + saw_header = true; + } else if trimmed.eq_ignore_ascii_case("[data]") { + saw_data = true; + } else if version.is_none() { + let lower = trimmed.to_ascii_lowercase(); + if let Some(rest) = lower.strip_prefix("content") { + if let Some(bpm) = rest + .split_whitespace() + .find(|tok| tok.ends_with(".bpm")) + { + version = Some(bpm.trim_end_matches(".bpm").to_owned()); + } + } else if lower.starts_with("gsgt version") { + version = trimmed + .split_whitespace() + .last() + .map(|v| format!("GSGT-{v}")); + } + } + } + if saw_header && saw_data { + Some(GsgtReport { + platform_version: version, + }) + } else { + None + } +} + fn extract_after_marker(text: &str, marker: &str) -> Option { text.lines().find_map(|line| { let trimmed = line.trim(); @@ -251,6 +313,8 @@ pub(crate) fn detect_assembly(lower_name: &str, sample_lines: &[String]) -> Opti || combined.contains("##contig= Result { let mut stats = SexStats::default(); let delimiter = detect_delimiter(lines); + let is_gsgt = lines_look_like_gsgt(lines); + let mut gsgt = None; let mut column_indexes = None; let mut comment_header = None; for line in lines { @@ -189,6 +194,8 @@ pub fn infer_sex_from_text_lines( line, kind, delimiter, + is_gsgt, + &mut gsgt, &mut column_indexes, &mut comment_header, )?; @@ -212,6 +219,8 @@ fn infer_sex_from_reader( let bytes = reader.read_line(&mut line).unwrap_or_default(); if bytes == 0 { let delimiter = detect_delimiter(&probe_lines); + let is_gsgt = lines_look_like_gsgt(&probe_lines); + let mut gsgt = None; let mut column_indexes = None; let mut comment_header = None; for probe_line in &probe_lines { @@ -220,6 +229,8 @@ fn infer_sex_from_reader( probe_line, kind, delimiter, + is_gsgt, + &mut gsgt, &mut column_indexes, &mut comment_header, )?; @@ -229,6 +240,8 @@ fn infer_sex_from_reader( probe_lines.push(line.trim_end_matches(['\n', '\r']).to_owned()); } let delimiter = detect_delimiter(&probe_lines); + let is_gsgt = lines_look_like_gsgt(&probe_lines); + let mut gsgt = None; let mut column_indexes = None; let mut comment_header = None; for probe_line in &probe_lines { @@ -237,6 +250,8 @@ fn infer_sex_from_reader( probe_line, kind, delimiter, + is_gsgt, + &mut gsgt, &mut column_indexes, &mut comment_header, )?; @@ -252,6 +267,8 @@ fn infer_sex_from_reader( line.trim_end_matches(['\n', '\r']), kind, delimiter, + is_gsgt, + &mut gsgt, &mut column_indexes, &mut comment_header, )?; @@ -259,11 +276,14 @@ fn infer_sex_from_reader( Ok(classify_stats(&stats, kind)) } +#[allow(clippy::too_many_arguments)] fn update_stats_from_line( stats: &mut SexStats, line: &str, kind: DetectedKind, delimiter: Delimiter, + is_gsgt: bool, + gsgt: &mut Option, column_indexes: &mut Option, comment_header: &mut Option>, ) -> Result<(), RuntimeError> { @@ -274,22 +294,39 @@ fn update_stats_from_line( match kind { DetectedKind::Vcf => update_vcf_stats(stats, trimmed), DetectedKind::GenotypeText | DetectedKind::Unknown => { - update_genotype_text_stats(stats, trimmed, delimiter, column_indexes, comment_header)?; + update_genotype_text_stats( + stats, trimmed, delimiter, is_gsgt, gsgt, column_indexes, comment_header, + )?; } _ => {} } Ok(()) } +#[allow(clippy::too_many_arguments)] fn update_genotype_text_stats( stats: &mut SexStats, line: &str, delimiter: Delimiter, + is_gsgt: bool, + gsgt: &mut Option, column_indexes: &mut Option, comment_header: &mut Option>, ) -> Result<(), RuntimeError> { - let Some(row) = parse_streaming_row(line, delimiter, column_indexes, comment_header)? else { - return Ok(()); + // GSGT SNP-array exports must normalize into the same + // rsid/chrom/pos/genotype row as every other SNP text export so the + // identical X/Y fingerprint heuristic applies regardless of source. + let row = if is_gsgt { + match gsgt.get_or_insert_with(GsgtParser::new).consume(line)? { + Some(row) => row, + None => return Ok(()), + } + } else { + let Some(row) = parse_streaming_row(line, delimiter, column_indexes, comment_header)? + else { + return Ok(()); + }; + row }; let rsid = row.rsid.as_deref().unwrap_or_default(); let chrom = normalize_chrom(row.chrom.as_deref().unwrap_or_default()); @@ -502,6 +539,36 @@ mod tests { assert_eq!(result.confidence, SexDetectionConfidence::High); } + #[test] + fn gsgt_snp_array_goes_through_the_same_y_fingerprint() { + // A male SNP array wrapped in the Illumina GSGT [Header]/[Data] + // layout must reach the identical X/Y fingerprint and classify the + // same as any other text export — source format is irrelevant. + let mut gsgt = vec![ + "[Header]".to_owned(), + "GSGT Version\t2.0.5".to_owned(), + "[Data]".to_owned(), + "SNP Name\tSNP\tChr\tPosition\tAllele1 - Plus\tAllele2 - Plus".to_owned(), + ]; + let mut flat = Vec::new(); + for i in 0..600 { + let rsid = format!("rs{}", 700_000 + i); + gsgt.push(format!("{rsid}\t[A/G]\tY\t{}\tG\tG", i + 1)); + flat.push(format!("{rsid}\tY\t{}\tG", i + 1)); + } + + let gsgt_result = + infer_sex_from_text_lines(&gsgt, DetectedKind::GenotypeText).unwrap(); + let flat_result = + infer_sex_from_text_lines(&flat, DetectedKind::GenotypeText).unwrap(); + + assert_eq!(gsgt_result.sex, InferredSex::Male); + assert_eq!(gsgt_result.method, "snp_array_x_y_fingerprint"); + // Identical heuristic outcome regardless of which export format. + assert_eq!(gsgt_result.sex, flat_result.sex); + assert_eq!(gsgt_result.confidence, flat_result.confidence); + } + #[test] fn y_fingerprint_uses_genotype_column_not_array_metrics() { let lines: Vec = (0..1001) diff --git a/rust/bioscript-formats/tests/file_formats/delimited.rs b/rust/bioscript-formats/tests/file_formats/delimited.rs index 15b450b..bcfa363 100644 --- a/rust/bioscript-formats/tests/file_formats/delimited.rs +++ b/rust/bioscript-formats/tests/file_formats/delimited.rs @@ -132,3 +132,231 @@ fn delimited_parser_handles_space_delimited_rows_without_headers_and_inline_comm observation.evidence ); } + +fn fixtures_dir() -> PathBuf { + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures") +} + +fn gsgt_sample_bytes() -> Vec { + fs::read(fixtures_dir().join("carigenetics_gsgt_sample.txt")).unwrap() +} + +#[test] +fn gsgt_in_memory_path_parses_merges_and_skips() { + let store = + GenotypeStore::from_bytes("carigenetics_gsgt_sample.txt", &gsgt_sample_bytes()).unwrap(); + + // Plain rs + Plus columns. + assert_eq!(store.get("rs9000001").unwrap().as_deref(), Some("GG")); + // Prefixed rs recovered from SNP Name. + assert_eq!(store.get("rs9000002").unwrap().as_deref(), Some("AG")); + // Replicate probes: a real call + a no-call replicate -> the call. + assert_eq!(store.get("rs9000003").unwrap().as_deref(), Some("CC")); + // Replicate disagreement -> no-call, never auto-picked. + assert_eq!(store.get("rs9000004").unwrap().as_deref(), Some("--")); + // Co-located distinct SNP stays addressable by its own rsid. + assert_eq!(store.get("rs9000007").unwrap().as_deref(), Some("AT")); + // No-call only. + assert_eq!(store.get("rs9000011").unwrap().as_deref(), Some("--")); + // Chr == 0 and Position == 0 rows are skipped entirely. + assert_eq!(store.get("rs9000009").unwrap(), None); + assert_eq!(store.get("rs9000010").unwrap(), None); + + // Non-rs probe resolves only by locus. + let obs = store + .lookup_variant(&VariantSpec { + grch38: Some(bioscript_core::GenomicLocus { + chrom: "2".to_owned(), + start: 6000, + end: 6000, + }), + ..VariantSpec::default() + }) + .unwrap(); + assert_eq!(obs.genotype.as_deref(), Some("CT")); +} + +#[test] +fn gsgt_streaming_scan_path_merges_replicates() { + let dir = temp_dir("gsgt-scan"); + let path = dir.join("carigenetics_gsgt_sample.txt"); + fs::write(&path, gsgt_sample_bytes()).unwrap(); + + let store = + GenotypeStore::from_file_with_options(&path, &GenotypeLoadOptions::default()).unwrap(); + + let by_rsid = |rsid: &str| { + store + .lookup_variant(&VariantSpec { + rsids: vec![rsid.to_owned()], + ..VariantSpec::default() + }) + .unwrap() + .genotype + }; + assert_eq!(by_rsid("rs9000001").as_deref(), Some("GG")); + // call + no-call replicate merged in a single streaming pass. + assert_eq!(by_rsid("rs9000003").as_deref(), Some("CC")); + // disagreeing replicates collapse to a no-call. + assert_eq!(by_rsid("rs9000004").as_deref(), Some("--")); +} + +#[test] +fn gsgt_matches_ddna_concordance_on_shared_fixture() { + let gsgt = + GenotypeStore::from_bytes("carigenetics_gsgt_sample.txt", &gsgt_sample_bytes()).unwrap(); + let ddna = GenotypeStore::from_bytes( + "carigenetics_ddna_concordance_sample.txt", + &fs::read(fixtures_dir().join("carigenetics_ddna_concordance_sample.txt")).unwrap(), + ) + .unwrap(); + + let mut compared = 0; + for rsid in ["rs9000001", "rs9000002", "rs9000003", "rs9000007"] { + let g = gsgt.get(rsid).unwrap().unwrap(); + let d = ddna.get(rsid).unwrap().unwrap(); + assert_eq!(g, d, "genotype mismatch for {rsid}: gsgt={g} ddna={d}"); + compared += 1; + } + assert_eq!(compared, 4); +} + +#[test] +fn ddna_path_unchanged_placeholder_rsids_not_collapsed_by_rsid() { + // Regression guard (illumina.md Phase 3): keying merges by rsid alone + // silently collapses every placeholder rsid into one row. DDNA must not + // go through the GSGT path and must keep distinct loci distinct. + let dir = temp_dir("ddna-placeholder"); + let path = dir.join("ddna.txt"); + fs::write( + &path, + "# rsid\tchromosome\tposition\tgenotype\n\ + .\t1\t100\tAA\n\ + .\t1\t200\tGG\n", + ) + .unwrap(); + + let store = GenotypeStore::from_file_with_options( + &path, + &GenotypeLoadOptions { + assembly: Some(Assembly::Grch38), + ..GenotypeLoadOptions::default() + }, + ) + .unwrap(); + + let at = |start: i64| { + store + .lookup_variant(&VariantSpec { + grch38: Some(bioscript_core::GenomicLocus { + chrom: "1".to_owned(), + start, + end: start, + }), + ..VariantSpec::default() + }) + .unwrap() + .genotype + }; + assert_eq!(at(100).as_deref(), Some("AA")); + assert_eq!(at(200).as_deref(), Some("GG")); +} + +fn carika_dir() -> Option { + if let Some(p) = std::env::var_os("BIOVAULT_CARIKA_DIR") { + let d = PathBuf::from(p); + if d.exists() { + return Some(d); + } + } + let default = PathBuf::from("/Users/madhavajay/dev/my_private_data/carika"); + default.exists().then_some(default) +} + +fn is_acgt_call(g: &str) -> bool { + g.len() == 2 && g.chars().all(|c| matches!(c, 'A' | 'C' | 'G' | 'T')) +} + +fn is_palindromic(g: &str) -> bool { + let mut s: Vec = g.chars().collect(); + s.sort_unstable(); + matches!(s.as_slice(), ['A', 'T'] | ['C', 'G']) +} + +// Phase 4 (illumina.md): prove the real PC0001 GSGT export and the +// already-supported DDNA export are genotype-concordant. Skips when the +// private data dir is unavailable so CI stays green. +#[test] +#[allow(clippy::cast_precision_loss)] +fn pc0001_gsgt_matches_ddna_real_files() { + let Some(dir) = carika_dir() else { + eprintln!("skipping pc0001_gsgt_matches_ddna_real_files: BIOVAULT_CARIKA_DIR unset"); + return; + }; + let gsgt_path = dir.join("PC0001_Raw Data_Carigenetics.txt"); + let ddna_path = dir.join("PC0001_X_X_GSAv3-DTC_GRCh38-07-29-2025.txt"); + if !gsgt_path.exists() || !ddna_path.exists() { + eprintln!("skipping pc0001_gsgt_matches_ddna_real_files: sample files missing"); + return; + } + + let gsgt = GenotypeStore::from_bytes( + "PC0001_Raw Data_Carigenetics.txt", + &fs::read(&gsgt_path).unwrap(), + ) + .unwrap(); + let ddna = GenotypeStore::from_bytes( + "PC0001_X_X_GSAv3-DTC_GRCh38-07-29-2025.txt", + &fs::read(&ddna_path).unwrap(), + ) + .unwrap(); + + // rsids come from the DDNA ground-truth file (col 0 == rsNNN). + let ddna_raw = fs::read_to_string(&ddna_path).unwrap(); + let mut compared = 0u64; + let mut concordant = 0u64; + let mut palindromic = 0u64; + let mut mismatch = 0u64; + for line in ddna_raw.lines() { + let line = line.trim(); + if line.is_empty() || line.starts_with('#') { + continue; + } + let Some(rsid) = line.split('\t').next() else { + continue; + }; + if !rsid.starts_with("rs") { + continue; + } + let (Some(g), Some(d)) = ( + gsgt.get(rsid).unwrap(), + ddna.get(rsid).unwrap(), + ) else { + continue; + }; + if !is_acgt_call(&g) || !is_acgt_call(&d) { + continue; + } + compared += 1; + if g == d { + concordant += 1; + } else if is_palindromic(&g) && is_palindromic(&d) { + palindromic += 1; + } else { + mismatch += 1; + } + } + + let pct = 100.0 * concordant as f64 / compared as f64; + eprintln!( + "PC0001 concordance: compared={compared} concordant={concordant} \ + palindromic={palindromic} mismatch={mismatch} ({pct:.4}%)", + ); + assert!(compared > 100_000, "too few comparable SNPs: {compared}"); + let non_palindromic_mismatch_rate = mismatch as f64 / compared as f64; + assert!( + non_palindromic_mismatch_rate < 0.001, + "non-palindromic mismatch rate {non_palindromic_mismatch_rate:.6} exceeds 0.1% \ + (compared={compared} mismatch={mismatch})", + ); +} diff --git a/rust/bioscript-formats/tests/fixtures/carigenetics_ddna_concordance_sample.txt b/rust/bioscript-formats/tests/fixtures/carigenetics_ddna_concordance_sample.txt new file mode 100644 index 0000000..4a412dd --- /dev/null +++ b/rust/bioscript-formats/tests/fixtures/carigenetics_ddna_concordance_sample.txt @@ -0,0 +1,8 @@ +# Synthetic Dynamic DNA (DDNA) export for the same fabricated individual as +# carigenetics_gsgt_sample.txt. Build 38, plus strand. Values are fake. +# rsid chromosome position genotype gs baf lrr +rs9000001 1 1000 GG 0.99 1.0 0.01 +rs9000002 1 2000 AG 0.98 0.5 0.02 +rs9000003 1 3000 CC 0.97 0.0 0.00 +rs9000007 3 8000 AT 0.99 0.5 0.01 +rs9000011 4 9000 -- 0.10 0.5 0.01 diff --git a/rust/bioscript-formats/tests/fixtures/carigenetics_gsgt_sample.txt b/rust/bioscript-formats/tests/fixtures/carigenetics_gsgt_sample.txt new file mode 100644 index 0000000..83e6124 --- /dev/null +++ b/rust/bioscript-formats/tests/fixtures/carigenetics_gsgt_sample.txt @@ -0,0 +1,31 @@ +[Header] +GSGT Version 2.0.5 +Processing Date 01/02/2030 1:23 PM +Content SyntheticPanel_00000000_A1.bpm +Num SNPs 14 +Total SNPs 14 +Num Samples 1 +Total Samples 1 +File 1 of 1 +[Data] +Sample ID Sample Name SNP Name SNP Chr Position Allele1 - Top Allele2 - Top Allele1 - Forward Allele2 - Forward Allele1 - Plus Allele2 - Plus Allele1 - Design Allele2 - Design Allele1 - AB Allele2 - AB Plus/Minus Strand +FAKE001 rs9000001 [A/G] 1 1000 G G G G G G G G B B + +FAKE001 BOT-rs9000002 [A/G] 1 2000 A G A G A G A G A B - +FAKE001 rs9000003_ilmndup1 [C/T] 1 3000 C C C C C C C C A A + +FAKE001 GSA-rs9000003 [C/T] 1 3000 - - - - - - C C - - + +FAKE001 rs9000004 [A/G] 2 5000 A A A A A A A A A A + +FAKE001 BOT-rs9000004 [A/G] 2 5000 G G G G G G A G B B - +FAKE001 2:6000 [C/T] 2 6000 C T C T C T C T A B + +FAKE001 1:7000_CNV_GENEX [A/G] 1 7000 A A A A A A A A A A + +FAKE001 seq-rs9000007 [A/T] 3 8000 A T A T A T A T A B + +FAKE001 1:8000-I-D_MNV [I/D] 3 8000 I D I D I D I D A B + +FAKE001 rs9000009 [A/G] 0 0 A G A G A G A G A B + +FAKE001 rs9000010 [A/G] 4 0 A G A G A G A G A B + +FAKE001 rs9000011 [A/G] 4 9000 - - - - - - A G - - + +FAKE001 DICER1-chr14-95000000 [A/G] 14 95000000 A A A A A A A A A A + +FAKE001 rs429358 [C/T] 19 44908684 C C C C C C C C A A + +FAKE001 rs7412 [C/T] 19 44908822 C C C C C C C C A A + +FAKE001 rs1801133 [A/G] 1 11796321 A G A G A G A G A B + +FAKE001 rs53576 [A/G] 3 8762685 A A A A A A A A A A + +FAKE001 rs9939609 [A/T] 16 53786615 A T A T A T A T A B + +FAKE001 rs4680 [A/G] 22 19963748 G G G G G G G G A A + diff --git a/rust/bioscript-formats/tests/inspect.rs b/rust/bioscript-formats/tests/inspect.rs index fd10d36..a1c6c3a 100644 --- a/rust/bioscript-formats/tests/inspect.rs +++ b/rust/bioscript-formats/tests/inspect.rs @@ -804,3 +804,134 @@ chr1\t100\trs1\tA\tG\t.\tPASS\t.\tGT\t0/1\n", assert_eq!(inspection.detected_kind, DetectedKind::Vcf); assert_eq!(inspection.assembly, Some(Assembly::Grch38)); } + +#[test] +fn inspect_bytes_detects_gsgt_carigenetics_final_report() { + let bytes = std::fs::read(fixtures_dir().join("carigenetics_gsgt_sample.txt")).unwrap(); + let inspection = + inspect_bytes("export.txt", &bytes, &InspectOptions::default()).unwrap(); + + assert_eq!(inspection.detected_kind, DetectedKind::GenotypeText); + // GSGT carries no build line but is GRCh38 (spec §2.1). + assert_eq!(inspection.assembly, Some(Assembly::Grch38)); + let source = inspection.source.expect("vendor detected"); + assert_eq!(source.vendor.as_deref(), Some("CariGenetics")); + assert!( + source + .evidence + .iter() + .any(|e| e.contains("GSGT [Header]/[Data] block")), + "evidence: {:?}", + source.evidence + ); + // Detection works on content alone, even with a non-CariGenetics name. +} + +#[test] +fn inspect_bytes_still_detects_dynamic_dna_no_regression() { + let bytes = + std::fs::read(fixtures_dir().join("carigenetics_ddna_concordance_sample.txt")).unwrap(); + let inspection = inspect_bytes( + "PC0001_X_X_GSAv3-DTC_GRCh38.txt", + &bytes, + &InspectOptions::default(), + ) + .unwrap(); + assert_eq!(inspection.detected_kind, DetectedKind::GenotypeText); + let source = inspection.source.expect("vendor detected"); + assert_eq!(source.vendor.as_deref(), Some("Dynamic DNA")); +} + +// Real-file proof: the same individual's GSGT and DDNA exports must infer +// the same sex through the identical fingerprint. Skips if private data absent. +#[test] +fn pc0001_gsgt_and_ddna_infer_same_sex_real_files() { + let dir = std::env::var_os("BIOVAULT_CARIKA_DIR") + .map(std::path::PathBuf::from) + .filter(|d| d.exists()) + .or_else(|| { + let d = std::path::PathBuf::from("/Users/madhavajay/dev/my_private_data/carika"); + d.exists().then_some(d) + }); + let Some(dir) = dir else { + eprintln!("skipping pc0001 sex parity: BIOVAULT_CARIKA_DIR unset"); + return; + }; + let gsgt_path = dir.join("PC0001_Raw Data_Carigenetics.txt"); + let ddna_path = dir.join("PC0001_X_X_GSAv3-DTC_GRCh38-07-29-2025.txt"); + if !gsgt_path.exists() || !ddna_path.exists() { + eprintln!("skipping pc0001 sex parity: sample files missing"); + return; + } + + let sex = |p: &std::path::Path| { + let f = std::fs::File::open(p).unwrap(); + bioscript_formats::infer_sex_from_named_reader( + &p.file_name().unwrap().to_string_lossy(), + f, + DetectedKind::GenotypeText, + ) + .unwrap() + }; + let g = sex(&gsgt_path); + let d = sex(&ddna_path); + eprintln!( + "PC0001 sex: gsgt={:?}/{:?} ddna={:?}/{:?}", + g.sex, g.confidence, d.sex, d.confidence + ); + assert_eq!(g.method, "snp_array_x_y_fingerprint"); + assert_eq!(g.sex, d.sex, "GSGT vs DDNA inferred different sex"); + assert_eq!(g.sex, InferredSex::Male); +} + +// Gated real-file matrix: every cached vendor export must resolve to its +// known ground-truth build (or honest unknown for build-36). Skips when the +// local test-data cache is absent. Proves metadata->anchor priority. +#[test] +fn assembly_matrix_matches_ground_truth_across_vendors() { + let cache = std::path::PathBuf::from(std::env::var_os("HOME").unwrap()) + .join(".bioscript/cache/test-data"); + let carika = std::path::PathBuf::from("/Users/madhavajay/dev/my_private_data/carika"); + if !cache.exists() { + eprintln!("skipping assembly_matrix: no test-data cache"); + return; + } + use bioscript_core::Assembly::{Grch37, Grch38}; + let cases: &[(&str, Option)] = &[ + ("23andme/v4/huE18D82/genome__v4_Full_2016.txt", Some(Grch37)), + ("23andme/v5/hu50B3F5/genome_Lisa_Fauman_v5_Full_20180326062517.txt", Some(Grch37)), + ("ancestrydna/huE922FC/AncestryDNA.txt", Some(Grch37)), + ("myheritage/hu33515F/MyHeritage_raw_dna_data.csv", Some(Grch37)), + ("genesforgood/hu80B047/GFG0_filtered_imputed_genotypes_noY_noMT_23andMe.txt", Some(Grch37)), + ("dynamicdna/100001-synthetic/100001_X_X_GSAv3-DTC_GRCh38-07-12-2025.txt", Some(Grch38)), + // build-36 files must NOT be misread as 37/38. + ("23andme/v2/hu0199C8/23data20100526.txt", None), + ("23andme/v3/huE4DAE4/huE4DAE4_20120522224129.txt", None), + ]; + for (rel, want) in cases { + let path = cache.join(rel); + if !path.exists() { + eprintln!(" skip missing {rel}"); + continue; + } + let bytes = std::fs::read(&path).unwrap(); + let got = inspect_bytes(rel, &bytes, &InspectOptions::default()) + .unwrap() + .assembly; + assert_eq!(got, *want, "vendor {rel}: got {got:?} want {want:?}"); + } + + // The GSGT file declares no build; it must resolve GRCh38 via the + // anchor vote (not a hardcode). + for f in ["PC0001_Raw Data_Carigenetics.txt", "PC0159_Raw Data_Carigenetics.txt"] { + let p = carika.join(f); + if !p.exists() { + continue; + } + let bytes = std::fs::read(&p).unwrap(); + let got = inspect_bytes(f, &bytes, &InspectOptions::default()) + .unwrap() + .assembly; + assert_eq!(got, Some(Grch38), "GSGT {f} should anchor-vote GRCh38"); + } +} From a6b9e0782fd286fac6096b71127ec293bb6ed605 Mon Sep 17 00:00:00 2001 From: Madhava Jay Date: Tue, 19 May 2026 13:51:28 +1000 Subject: [PATCH 2/2] linting --- .../src/genotype/delimited/gsgt.rs | 30 +++++-- .../src/genotype/delimited/scan.rs | 8 +- .../bioscript-formats/src/genotype/loaders.rs | 4 +- rust/bioscript-formats/src/inspect.rs | 22 +++-- .../src/inspect/assembly_anchors.rs | 34 +++++-- .../src/inspect/heuristics.rs | 5 +- rust/bioscript-formats/src/inspect/sex.rs | 88 +++---------------- .../src/inspect/sex/calls.rs | 68 ++++++++++++++ .../tests/file_formats/delimited.rs | 5 +- rust/bioscript-formats/tests/inspect.rs | 30 +++++-- 10 files changed, 174 insertions(+), 120 deletions(-) create mode 100644 rust/bioscript-formats/src/inspect/sex/calls.rs diff --git a/rust/bioscript-formats/src/genotype/delimited/gsgt.rs b/rust/bioscript-formats/src/genotype/delimited/gsgt.rs index 66a4069..bc67f4f 100644 --- a/rust/bioscript-formats/src/genotype/delimited/gsgt.rs +++ b/rust/bioscript-formats/src/genotype/delimited/gsgt.rs @@ -122,7 +122,9 @@ fn split_tab(line: &str) -> Vec<&str> { fn resolve_columns(header: &[&str]) -> Result { let find = |want: &str| header.iter().position(|h| normalize_name(h) == want); let snp_name = find("snpname"); - let chr = find("chr").or_else(|| find("chromosome")).or_else(|| find("chrom")); + let chr = find("chr") + .or_else(|| find("chromosome")) + .or_else(|| find("chrom")); let position = find("position").or_else(|| find("pos")); let allele1_plus = find("allele1plus"); let allele2_plus = find("allele2plus"); @@ -143,11 +145,7 @@ fn resolve_columns(header: &[&str]) -> Result { } } -fn extract_row( - line: &str, - fields: &[&str], - cols: &GsgtColumns, -) -> Option { +fn extract_row(line: &str, fields: &[&str], cols: &GsgtColumns) -> Option { let snp_name = fields.get(cols.snp_name).copied().unwrap_or(""); let rsid = extract_rsid(snp_name); @@ -217,8 +215,14 @@ mod tests { extract_rsid("rs111647200_ilmndup1").as_deref(), Some("rs111647200") ); - assert_eq!(extract_rsid("GSA-rs61660502").as_deref(), Some("rs61660502")); - assert_eq!(extract_rsid("seq-rs786202193").as_deref(), Some("rs786202193")); + assert_eq!( + extract_rsid("GSA-rs61660502").as_deref(), + Some("rs61660502") + ); + assert_eq!( + extract_rsid("seq-rs786202193").as_deref(), + Some("rs786202193") + ); // No rs id: chr:pos, CNV, MNV, vendor/clinical. assert_eq!(extract_rsid("1:103380393"), None); assert_eq!(extract_rsid("1:159174749-C-T"), None); @@ -241,7 +245,15 @@ mod tests { p } - fn row(p: &mut GsgtParser, snp_name: &str, snp: &str, chr: &str, pos: &str, a1: &str, a2: &str) -> Option { + fn row( + p: &mut GsgtParser, + snp_name: &str, + snp: &str, + chr: &str, + pos: &str, + a1: &str, + a2: &str, + ) -> Option { let line = format!("S1\t\t{snp_name}\t{snp}\t{chr}\t{pos}\tX\tX\t{a1}\t{a2}\t+"); p.consume(&line).unwrap() } diff --git a/rust/bioscript-formats/src/genotype/delimited/scan.rs b/rust/bioscript-formats/src/genotype/delimited/scan.rs index 039ad36..bb7f3a6 100644 --- a/rust/bioscript-formats/src/genotype/delimited/scan.rs +++ b/rust/bioscript-formats/src/genotype/delimited/scan.rs @@ -86,10 +86,10 @@ fn prescan_assembly_anchors( let mut ch: Option> = None; let mut gsgt = GsgtParser::new(); let feed = |line: &str, - scorer: &mut AssemblyAnchorScorer, - gsgt: &mut GsgtParser, - ci: &mut Option, - ch: &mut Option>| + scorer: &mut AssemblyAnchorScorer, + gsgt: &mut GsgtParser, + ci: &mut Option, + ch: &mut Option>| -> Result<(), RuntimeError> { let row = if is_gsgt { gsgt.consume(line)? diff --git a/rust/bioscript-formats/src/genotype/loaders.rs b/rust/bioscript-formats/src/genotype/loaders.rs index bb71808..a8759a5 100644 --- a/rust/bioscript-formats/src/genotype/loaders.rs +++ b/rust/bioscript-formats/src/genotype/loaders.rs @@ -140,8 +140,8 @@ fn from_gsgt_reader( let mut groups: HashMap<(String, i64, String), GsgtMergeGroup> = HashMap::new(); let ingest = |row: super::delimited::ParsedDelimitedRow, - order: &mut Vec<(String, i64, String)>, - groups: &mut HashMap<(String, i64, String), GsgtMergeGroup>| { + order: &mut Vec<(String, i64, String)>, + groups: &mut HashMap<(String, i64, String), GsgtMergeGroup>| { let (Some(chrom), Some(position)) = (row.chrom.clone(), row.position) else { return; }; diff --git a/rust/bioscript-formats/src/inspect.rs b/rust/bioscript-formats/src/inspect.rs index a687312..9cc7520 100644 --- a/rust/bioscript-formats/src/inspect.rs +++ b/rust/bioscript-formats/src/inspect.rs @@ -184,9 +184,12 @@ pub fn inspect_bytes( // Declared metadata first; fall back to the rsID/locus anchor vote over // the buffer only when the file declares no build (e.g. a GSGT report). let assembly = detect_assembly(&inspection_context, &sample_lines).or_else(|| { - matches!(detected_kind, DetectedKind::GenotypeText | DetectedKind::Vcf) - .then(|| assembly_from_text_bytes(bytes)) - .flatten() + matches!( + detected_kind, + DetectedKind::GenotypeText | DetectedKind::Vcf + ) + .then(|| assembly_from_text_bytes(bytes)) + .flatten() }); let phased = (detected_kind == DetectedKind::Vcf) .then(|| detect_vcf_phasing(&sample_lines)) @@ -291,9 +294,16 @@ pub fn inspect_file(path: &Path, options: &InspectOptions) -> Result Option { /// Normalize a chromosome label to bare form (`chr7` -> `7`, `23` -> `X`). fn norm_chrom(raw: &str) -> String { let c = raw.trim(); - let c = c.strip_prefix("chr").or_else(|| c.strip_prefix("CHR")).unwrap_or(c); + let c = c + .strip_prefix("chr") + .or_else(|| c.strip_prefix("CHR")) + .unwrap_or(c); match c.to_ascii_uppercase().as_str() { "23" => "X".to_owned(), "24" => "Y".to_owned(), @@ -91,7 +108,10 @@ fn allele_ok(genotype: &str, anchor_alleles: &str) -> bool { if observed.is_empty() { return true; } - let fwd: Vec = anchor_alleles.bytes().map(|b| b.to_ascii_uppercase()).collect(); + let fwd: Vec = anchor_alleles + .bytes() + .map(|b| b.to_ascii_uppercase()) + .collect(); let rev: Vec = fwd.iter().filter_map(|b| complement(*b)).collect(); let subset = |set: &[u8]| observed.iter().all(|o| set.contains(o)); subset(&fwd) || subset(&rev) @@ -203,11 +223,7 @@ pub(crate) fn assembly_from_text_bytes(bytes: &[u8]) -> Option { } } Mode::GsgtHdr => { - let norm = |s: &str| { - s.trim() - .to_ascii_lowercase() - .replace([' ', '-', '_'], "") - }; + let norm = |s: &str| s.trim().to_ascii_lowercase().replace([' ', '-', '_'], ""); let h: Vec = line.split('\t').map(norm).collect(); let idx = |w: &str| h.iter().position(|x| x == w); match ( diff --git a/rust/bioscript-formats/src/inspect/heuristics.rs b/rust/bioscript-formats/src/inspect/heuristics.rs index 37d4ab2..1c03a2b 100644 --- a/rust/bioscript-formats/src/inspect/heuristics.rs +++ b/rust/bioscript-formats/src/inspect/heuristics.rs @@ -242,10 +242,7 @@ fn detect_gsgt_final_report(sample_lines: &[String]) -> Option { } else if version.is_none() { let lower = trimmed.to_ascii_lowercase(); if let Some(rest) = lower.strip_prefix("content") { - if let Some(bpm) = rest - .split_whitespace() - .find(|tok| tok.ends_with(".bpm")) - { + if let Some(bpm) = rest.split_whitespace().find(|tok| tok.ends_with(".bpm")) { version = Some(bpm.trim_end_matches(".bpm").to_owned()); } } else if lower.starts_with("gsgt version") { diff --git a/rust/bioscript-formats/src/inspect/sex.rs b/rust/bioscript-formats/src/inspect/sex.rs index 43f76d5..2623bfd 100644 --- a/rust/bioscript-formats/src/inspect/sex.rs +++ b/rust/bioscript-formats/src/inspect/sex.rs @@ -15,11 +15,16 @@ use crate::genotype::{ use super::{DetectedKind, InspectOptions}; mod alignment_depth; +mod calls; mod classify; pub use alignment_depth::infer_sex_from_alignment_reader; pub(crate) use alignment_depth::infer_sex_from_alignment_path; +use calls::{ + genotype_allele_count, is_called_genotype_text, is_called_vcf_gt, is_genotype_text_het, + is_non_par_x, is_vcf_gt_het, normalize_chrom, vcf_gt_allele_count, +}; use classify::{classify_stats, supports_sex_detection, unsupported_sex_inference}; const MAX_SEX_DETECTION_LINES: usize = 50_000_000; @@ -295,7 +300,13 @@ fn update_stats_from_line( DetectedKind::Vcf => update_vcf_stats(stats, trimmed), DetectedKind::GenotypeText | DetectedKind::Unknown => { update_genotype_text_stats( - stats, trimmed, delimiter, is_gsgt, gsgt, column_indexes, comment_header, + stats, + trimmed, + delimiter, + is_gsgt, + gsgt, + column_indexes, + comment_header, )?; } _ => {} @@ -402,75 +413,6 @@ fn update_vcf_stats(stats: &mut SexStats, line: &str) { } } -fn normalize_chrom(value: &str) -> String { - let normalized = value - .trim() - .trim_start_matches("chr") - .trim_start_matches("CHR") - .to_ascii_uppercase(); - match normalized.as_str() { - "23" => "X".to_owned(), - "24" => "Y".to_owned(), - "25" => "XY".to_owned(), - "26" | "M" => "MT".to_owned(), - _ => normalized, - } -} - -fn is_called_genotype_text(value: &str) -> bool { - let value = value.trim(); - if value.is_empty() || matches!(value, "--" | "00" | "." | "./." | ".|.") { - return false; - } - value - .chars() - .all(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) -} - -fn genotype_allele_count(value: &str) -> usize { - value - .chars() - .filter(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) - .count() -} - -fn is_genotype_text_het(value: &str) -> bool { - let alleles: Vec = value - .chars() - .filter(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) - .map(|ch| ch.to_ascii_uppercase()) - .collect(); - alleles.len() == 2 && alleles[0] != alleles[1] -} - -fn is_called_vcf_gt(value: &str) -> bool { - let value = value.trim(); - !value.is_empty() - && !value.contains('.') - && (value != "0" || matches!(value, "0" | "1" | "2" | "3")) -} - -fn vcf_gt_allele_count(gt: &str) -> usize { - gt.split(['/', '|']) - .filter(|part| !part.is_empty() && *part != ".") - .count() -} - -fn is_vcf_gt_het(gt: &str) -> bool { - let alleles: Vec<&str> = gt - .split(['/', '|']) - .filter(|part| !part.is_empty() && *part != ".") - .collect(); - alleles.len() == 2 && alleles[0] != alleles[1] -} - -fn is_non_par_x(pos: u32) -> bool { - // Human GRCh38 non-PAR X used by bcftools +guess-ploidy. - // GRCh37 differs slightly, but these bounds cover the common non-PAR body - // and avoid both pseudoautosomal ends for this QC heuristic. - (2_781_480..=154_931_043).contains(&pos) -} - fn select_sex_detection_zip_entry( archive: &mut ZipArchive, ) -> Result { @@ -557,10 +499,8 @@ mod tests { flat.push(format!("{rsid}\tY\t{}\tG", i + 1)); } - let gsgt_result = - infer_sex_from_text_lines(&gsgt, DetectedKind::GenotypeText).unwrap(); - let flat_result = - infer_sex_from_text_lines(&flat, DetectedKind::GenotypeText).unwrap(); + let gsgt_result = infer_sex_from_text_lines(&gsgt, DetectedKind::GenotypeText).unwrap(); + let flat_result = infer_sex_from_text_lines(&flat, DetectedKind::GenotypeText).unwrap(); assert_eq!(gsgt_result.sex, InferredSex::Male); assert_eq!(gsgt_result.method, "snp_array_x_y_fingerprint"); diff --git a/rust/bioscript-formats/src/inspect/sex/calls.rs b/rust/bioscript-formats/src/inspect/sex/calls.rs new file mode 100644 index 0000000..b757677 --- /dev/null +++ b/rust/bioscript-formats/src/inspect/sex/calls.rs @@ -0,0 +1,68 @@ +pub(super) fn normalize_chrom(value: &str) -> String { + let normalized = value + .trim() + .trim_start_matches("chr") + .trim_start_matches("CHR") + .to_ascii_uppercase(); + match normalized.as_str() { + "23" => "X".to_owned(), + "24" => "Y".to_owned(), + "25" => "XY".to_owned(), + "26" | "M" => "MT".to_owned(), + _ => normalized, + } +} + +pub(super) fn is_called_genotype_text(value: &str) -> bool { + let value = value.trim(); + if value.is_empty() || matches!(value, "--" | "00" | "." | "./." | ".|.") { + return false; + } + value + .chars() + .all(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) +} + +pub(super) fn genotype_allele_count(value: &str) -> usize { + value + .chars() + .filter(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) + .count() +} + +pub(super) fn is_genotype_text_het(value: &str) -> bool { + let alleles: Vec = value + .chars() + .filter(|ch| matches!(ch.to_ascii_uppercase(), 'A' | 'C' | 'G' | 'T')) + .map(|ch| ch.to_ascii_uppercase()) + .collect(); + alleles.len() == 2 && alleles[0] != alleles[1] +} + +pub(super) fn is_called_vcf_gt(value: &str) -> bool { + let value = value.trim(); + !value.is_empty() + && !value.contains('.') + && (value != "0" || matches!(value, "0" | "1" | "2" | "3")) +} + +pub(super) fn vcf_gt_allele_count(gt: &str) -> usize { + gt.split(['/', '|']) + .filter(|part| !part.is_empty() && *part != ".") + .count() +} + +pub(super) fn is_vcf_gt_het(gt: &str) -> bool { + let alleles: Vec<&str> = gt + .split(['/', '|']) + .filter(|part| !part.is_empty() && *part != ".") + .collect(); + alleles.len() == 2 && alleles[0] != alleles[1] +} + +pub(super) fn is_non_par_x(pos: u32) -> bool { + // Human GRCh38 non-PAR X used by bcftools +guess-ploidy. + // GRCh37 differs slightly, but these bounds cover the common non-PAR body + // and avoid both pseudoautosomal ends for this QC heuristic. + (2_781_480..=154_931_043).contains(&pos) +} diff --git a/rust/bioscript-formats/tests/file_formats/delimited.rs b/rust/bioscript-formats/tests/file_formats/delimited.rs index bcfa363..ef7bca0 100644 --- a/rust/bioscript-formats/tests/file_formats/delimited.rs +++ b/rust/bioscript-formats/tests/file_formats/delimited.rs @@ -328,10 +328,7 @@ fn pc0001_gsgt_matches_ddna_real_files() { if !rsid.starts_with("rs") { continue; } - let (Some(g), Some(d)) = ( - gsgt.get(rsid).unwrap(), - ddna.get(rsid).unwrap(), - ) else { + let (Some(g), Some(d)) = (gsgt.get(rsid).unwrap(), ddna.get(rsid).unwrap()) else { continue; }; if !is_acgt_call(&g) || !is_acgt_call(&d) { diff --git a/rust/bioscript-formats/tests/inspect.rs b/rust/bioscript-formats/tests/inspect.rs index a1c6c3a..ef7a8c4 100644 --- a/rust/bioscript-formats/tests/inspect.rs +++ b/rust/bioscript-formats/tests/inspect.rs @@ -808,8 +808,7 @@ chr1\t100\trs1\tA\tG\t.\tPASS\t.\tGT\t0/1\n", #[test] fn inspect_bytes_detects_gsgt_carigenetics_final_report() { let bytes = std::fs::read(fixtures_dir().join("carigenetics_gsgt_sample.txt")).unwrap(); - let inspection = - inspect_bytes("export.txt", &bytes, &InspectOptions::default()).unwrap(); + let inspection = inspect_bytes("export.txt", &bytes, &InspectOptions::default()).unwrap(); assert_eq!(inspection.detected_kind, DetectedKind::GenotypeText); // GSGT carries no build line but is GRCh38 (spec §2.1). @@ -889,6 +888,7 @@ fn pc0001_gsgt_and_ddna_infer_same_sex_real_files() { // local test-data cache is absent. Proves metadata->anchor priority. #[test] fn assembly_matrix_matches_ground_truth_across_vendors() { + use bioscript_core::Assembly::{Grch37, Grch38}; let cache = std::path::PathBuf::from(std::env::var_os("HOME").unwrap()) .join(".bioscript/cache/test-data"); let carika = std::path::PathBuf::from("/Users/madhavajay/dev/my_private_data/carika"); @@ -896,14 +896,25 @@ fn assembly_matrix_matches_ground_truth_across_vendors() { eprintln!("skipping assembly_matrix: no test-data cache"); return; } - use bioscript_core::Assembly::{Grch37, Grch38}; let cases: &[(&str, Option)] = &[ ("23andme/v4/huE18D82/genome__v4_Full_2016.txt", Some(Grch37)), - ("23andme/v5/hu50B3F5/genome_Lisa_Fauman_v5_Full_20180326062517.txt", Some(Grch37)), + ( + "23andme/v5/hu50B3F5/genome_Lisa_Fauman_v5_Full_20180326062517.txt", + Some(Grch37), + ), ("ancestrydna/huE922FC/AncestryDNA.txt", Some(Grch37)), - ("myheritage/hu33515F/MyHeritage_raw_dna_data.csv", Some(Grch37)), - ("genesforgood/hu80B047/GFG0_filtered_imputed_genotypes_noY_noMT_23andMe.txt", Some(Grch37)), - ("dynamicdna/100001-synthetic/100001_X_X_GSAv3-DTC_GRCh38-07-12-2025.txt", Some(Grch38)), + ( + "myheritage/hu33515F/MyHeritage_raw_dna_data.csv", + Some(Grch37), + ), + ( + "genesforgood/hu80B047/GFG0_filtered_imputed_genotypes_noY_noMT_23andMe.txt", + Some(Grch37), + ), + ( + "dynamicdna/100001-synthetic/100001_X_X_GSAv3-DTC_GRCh38-07-12-2025.txt", + Some(Grch38), + ), // build-36 files must NOT be misread as 37/38. ("23andme/v2/hu0199C8/23data20100526.txt", None), ("23andme/v3/huE4DAE4/huE4DAE4_20120522224129.txt", None), @@ -923,7 +934,10 @@ fn assembly_matrix_matches_ground_truth_across_vendors() { // The GSGT file declares no build; it must resolve GRCh38 via the // anchor vote (not a hardcode). - for f in ["PC0001_Raw Data_Carigenetics.txt", "PC0159_Raw Data_Carigenetics.txt"] { + for f in [ + "PC0001_Raw Data_Carigenetics.txt", + "PC0159_Raw Data_Carigenetics.txt", + ] { let p = carika.join(f); if !p.exists() { continue;