diff --git a/CHANGELOG.md b/CHANGELOG.md index dc7167ea..b2cd0f58 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n ## Unreleased +- The Rust `og` binary reads the clinical sections of summary.json: `og summary collect --only b`. That covers the ClinVar screen (stars, call floor, gnomAD frequencies and rsIDs from step 13's indexed VEP output, the cDNA check, scan and deep intronic tables), the carrier list, PharmCAT, CPIC, pypgx, Cyrius, CPSR, the clinical variant filter with its ACMG SF v3.3 tier and the two-hit rule, slivar with its phased compound-het pairs, and the CYP2D6 block. Nothing calls it yet: the reports still read the summary that `bin/collect_summary.py` writes. The sections are the same as the Python's on 13 parity cases, compared as JSON with key order; a fourteenth holds the deliberate fix below. Quoted headers in SMN, interval and dosage tables preserve deletion evidence. Non-ASCII consensus contigs leave the collector running. Each input is read once: the ClinVar hits feed four sections, and the clinical VCF is read line by line in one pass for its record count and its HIGH-impact records. On a synthetic folder with 300 ClinVar hits, a 200,000-record clinical VCF and a 200,000-row slivar table, it takes 0.3 s and 51 MB, where the Python collector took 1.3 s and 135 MB. With 50,000 ClinVar hits (far more than a genome has) it takes 0.6 s but 460 MB, where the Python took 3.5 s and 113 MB, because og builds the whole summary as a JSON tree before it writes it. What differs from the Python: + - A clinical or slivar VCF named `.vcf.gz` that is plain text is read. The Python marked the section unreadable. + - A PharmCAT report that is not UTF-8 is still unreadable, but the note names the fault in Rust's words. + - The three builders of the committed asset tables now run in the Rust `og` binary: `og dev build cn-intervals` (the interval rows of `bin/cn_loci.tsv`), `og dev build callability` (`assets/callability_genes.tsv`) and `og dev build sv-rule-regions` (`assets/sv_rule_regions.tsv`). They replace `scripts/ci/build-cn-intervals.py`, which is gone (the script is now `scripts/ci/build-cn-intervals.sh`), and the Python inside `scripts/ci/build-callability.sh` and `scripts/ci/build-sv-rule-regions.sh`. The scripts still download the pinned files, and the three tables are the same, byte for byte, as the committed ones. With the downloads in the cache, the callability table takes 2.5 s and 27 MB instead of 5.5 s and 39 MB, and the SV rule table 4.9 s and 30 MB instead of 10.5 s and 217 MB. What changes on purpose: - The three scripts share one download function, `scripts/ci/lib/fetch-pinned.sh`. A file the cache folder already holds with the right sha256 is read where it is, not copied. A 403 or 404 fails at once, including after an earlier retry. Other failures, stalled transfers included, get at most six retries, 30 s apart, within a five-minute retry window. `build-cn-intervals.py` tried three times, 30 s apart. - The marker line of `bin/cn_loci.tsv` names `scripts/ci/build-cn-intervals.sh`. No row changes. @@ -19,6 +23,7 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n - `versions.env` is read, not run in bash, and a line that is not `NAME="value"` stops the check with exit 2. - Two `*_IMAGE` variables with the same value are both named, joined by `|`. The Python named only the alphabetically last one. No two share a value today. - `scripts/ci/image-smoke.sh` probes a changed image's helpers only when `bin/og` is built, and says so when it is not. Guard's `helper-binaries` job still probes the helpers of every image on each pull request. + - The Rust `og` binary has a first part of the summary collector: `og summary collect --only a`. It reads seven sections of summary.json (coverage, callability, sex_check, sample_qc, mosaic_sex_loss, read_count and variants), the run status, the manifest and each section's state. Nothing calls it yet: the reports still read the summary that `bin/collect_summary.py` writes. The sections are the same as the Python's on 33 parity cases, compared as JSON with key order. It walks the sample folder once instead of running a `glob` for each candidate path, and it reads step 33's table once for its two sections. On a synthetic 1,000,000-record VCF the variant count takes 0.05 s and 7 MB, where the Python collector took 0.9 s and 38 MB. What differs from the Python: - A file named `.vcf.gz` that is plain text is read. The Python marked the variants section unreadable. - A sample name with `*`, `?` or `[` is matched as written. The Python read it as a file pattern. diff --git a/docs/rust.md b/docs/rust.md index c624d8f9..e7f77e24 100644 --- a/docs/rust.md +++ b/docs/rust.md @@ -22,7 +22,8 @@ rust/ crates/og-setup/ what setup.sh runs once: the STRchive catalog, zips, the genetic map crates/og-steps/ the step helpers: constraint-join, prs format, step 40's consumer export; og family (step 38); og pgx outside-calls|hla-concordance|pypgx-summary (steps 08, 32, 36); og paraphase (step 35b) crates/og-summary/ the step tables of steps 20, 25 and 33: og mito, og prs table, og qc sample; - og summary collect --only a, the first seven summary sections (no caller yet) + og summary collect --only a, the first seven summary sections (no caller yet); + --only b, the clinical sections: ClinVar, carrier, PGx, CPSR, ACMG SF, two hits crates/og-pedigree/ the pedigree rules: families files, members, kinship, PARs crates/og-chip/ step 49, clonal haematopoiesis: og chip prepare, report, check-table crates/og-cn/ copy number from depth (steps 32, 35, 41, 46): og cn bed, call, intervals, check-bam, marker-regions, marker-read diff --git a/rust/Cargo.lock b/rust/Cargo.lock index 44bf8f81..352e0313 100644 --- a/rust/Cargo.lock +++ b/rust/Cargo.lock @@ -840,6 +840,8 @@ dependencies = [ "flate2", "indexmap", "memchr", + "og-bam", + "og-clinvar", "og-io", "og-steps", "pycompat", diff --git a/rust/README.md b/rust/README.md index 394e83e1..c42b40b6 100644 --- a/rust/README.md +++ b/rust/README.md @@ -10,7 +10,7 @@ OwnGenome's own tools, ported from `bin/*.py` one at a time into one static bina | `crates/og-bam` | BAM reading as pysam does it (fetch, `reference_end`, `cigarstring`, `get_blocks`), on htslib linked statically, and the tools that read BAMs: the reads at a variant (`reads`: one CIGAR walk per read, SA pieces, filters per tool), `og cdna check`, `og cdna scan`, `og comphet phase`, `og bam contigs`, the CYP2D6 depth check and the RFC1 motif count | | `crates/og-pedigree` | the pedigree rules of step 38: families files, members, kinship, declared relationships, PARs | | `crates/og-qc` | `og qc sex-loss`: step 33's mosaic loss of Y or X, with CPython's medians, `fmean` and exact `variance` | -| `crates/og-summary` | the step tables: `og qc sample` (step 33), `og prs table` (step 25), `og mito annotate`, `og mito deletions`, `og mito copies` (step 20); the first sections of summary.json, `og summary collect --only a`, beside the Python collector, which still writes it | +| `crates/og-summary` | the step tables: `og qc sample` (step 33), `og prs table` (step 25), `og mito annotate`, `og mito deletions`, `og mito copies` (step 20); the sections of summary.json, `og summary collect --only a` (QC) and `--only b` (ClinVar, carrier list, PGx, CPSR, the clinical filter with ACMG SF and two hits, slivar), beside the Python collector, which still writes it | | `crates/og-chip` | step 49, clonal haematopoiesis: `og chip prepare`, `og chip report`, `og chip check-table` | | `crates/og-cn` | copy number from depth, steps 32, 35, 41 and 46: `og cn bed`, `og cn call`, `og cn intervals`, `og cn check-bam`, `og cn marker-regions`, `og cn marker-read` | | `crates/og-bloodtype` | step 41, blood groups: `og bloodtype report`, `og bloodtype check-markers`, `og bloodtype systems` | diff --git a/rust/crates/og-summary/Cargo.toml b/rust/crates/og-summary/Cargo.toml index 1a7fa879..d1758441 100644 --- a/rust/crates/og-summary/Cargo.toml +++ b/rust/crates/og-summary/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "og-summary" -description = "The summary cluster: the step tools that turn a third-party tool's output into one of OwnGenome's tables (step 33's sample QC verdict, step 25's PRS table, step 20's chrM tables)" +description = "The summary cluster: the step tools that turn a third-party tool's output into one of OwnGenome's tables (step 33's sample QC verdict, step 25's PRS table, step 20's chrM tables), and og summary collect, summary.json section by section" version.workspace = true edition.workspace = true license.workspace = true @@ -9,6 +9,8 @@ publish.workspace = true [dependencies] pycompat.workspace = true og-io.workspace = true +og-bam.workspace = true +og-clinvar.workspace = true og-steps.workspace = true indexmap.workspace = true memchr.workspace = true diff --git a/rust/crates/og-summary/src/collect.rs b/rust/crates/og-summary/src/collect.rs index 0ae77923..15a29828 100644 --- a/rust/crates/og-summary/src/collect.rs +++ b/rust/crates/og-summary/src/collect.rs @@ -7,7 +7,9 @@ //! switches to it. Group `a` is the model and the first seven sections: //! coverage, callability, sex_check, sample_qc, mosaic_sex_loss, //! read_count and variants, with the run status, the manifest and each -//! section's state. +//! section's state. Group `b` is the clinical sections ([`group_b`]): +//! clinvar, carrier, pharmcat, cpic, pypgx, cyrius, cpsr, clinical (with +//! the ACMG SF tier and the two-hit rule) and slivar, and the CYP2D6 block. //! //! How it reads a folder: one walk of the folder ([`scan::FolderIndex`]) //! instead of a `glob` per candidate path; each input read once (the @@ -16,10 +18,20 @@ //! sections are read. A reader's error stays in its section, which is //! then `unreadable` with a note that names the file. +pub mod acmg; +pub mod assets; +pub mod carrier; +pub mod clinical; +pub mod clinvar; +pub mod cpsr; +pub mod group_b; pub mod model; +pub mod pgx; pub mod scan; pub mod sections; +pub mod text; pub mod time; +pub mod two_hits; use crate::args::{self, Kind, Opt}; use model::{GROUP_A, Manifest, RunStatus, Section, Spec, Summary}; @@ -29,7 +41,7 @@ use sections::{Callability, Coverage, Data, Mosaic, ReadCount, SampleQc, SexChec use std::io::Write; use std::path::Path; -const USAGE: &str = "usage: og summary collect --sample SAMPLE --sample-dir SAMPLE_DIR --only a +const USAGE: &str = "usage: og summary collect --sample SAMPLE --sample-dir SAMPLE_DIR --only {a,b} --out OUT [--declared-sex DECLARED_SEX] [--assets ASSETS] "; @@ -39,7 +51,7 @@ const OPTS: &[Opt] = &[ Opt("--declared-sex", Kind::Value, false), Opt("--out", Kind::Value, true), Opt("--assets", Kind::Value, false), - Opt("--only", Kind::Choice(&["a"]), true), + Opt("--only", Kind::Choice(&["a", "b"]), true), ]; /// What one collect reads. @@ -75,12 +87,21 @@ pub fn cli(argv: &[String], out: &mut impl Write, err: &mut impl Write) -> u8 { // A part of the sections must never replace the sample's complete // summary.json, so there is no default path while `--only` exists. let target = a.text("--out").to_owned(); - let summary = match collect(&input) { + let read = match a.text("--only") { + "b" => group_b::collect, + _ => collect, + }; + let summary = match read(&input) { Ok(s) => s, Err(e) => return crate::finish(Err(e), err), }; + // The sections' own data is freed before the text is made, and the + // tree before the file is written: a large ClinVar list is in memory in + // one form at a time. + let value = summary.to_json(); + drop(summary); let body = match dumps( - &summary.to_json(), + &value, &Dump { indent: Some(1), ..Dump::default() @@ -89,6 +110,7 @@ pub fn cli(argv: &[String], out: &mut impl Write, err: &mut impl Write) -> u8 { Ok(b) => b + "\n", Err(e) => return crate::finish(Err(e.to_string()), err), }; + drop(value); let r = crate::write_file(&target, body.as_bytes()); if r.is_ok() { let _ = writeln!(out, "Summary written: {target}"); @@ -195,6 +217,7 @@ pub fn collect(input: &Input) -> Result { manifest, declared_sex, sections, + cyp2d6: None, }) } diff --git a/rust/crates/og-summary/src/collect/acmg.rs b/rust/crates/og-summary/src/collect/acmg.rs new file mode 100644 index 00000000..db307729 --- /dev/null +++ b/rust/crates/og-summary/src/collect/acmg.rs @@ -0,0 +1,923 @@ +//! The ACMG secondary-findings tier of the clinical section: the variants +//! that are secondary findings by the rule of their gene +//! (`assets/acmg_sf.tsv`, ACMG SF v3.3, Lee et al., Genet Med 2025). +//! +//! Three sources: step 06's ClinVar P/LP hits, the clinical filter's rare +//! HIGH-impact records (step 23, ranked MANE Select, then canonical, then +//! any transcript, as `og vcf csq-rank` ranks them) and step 05's class 4 +//! and 5 SVs in the genes whose rule allows them. A variant of an SF gene +//! that its rule does not report (one MUTYH variant, a heterozygous HFE +//! p.C282Y) is listed under `not_reported` with the reason. A called +//! deletion counts as the second hit of a gene whose rule needs two +//! ([`two_hits`](super::two_hits)), unless read-backed phase puts both on +//! one copy. + +use super::assets::{ACMG_SF_VERSION, Key, Rule, acmg_sf_rules, carrier_genes}; +use super::clinvar::{ + COMMON_AF, Floor, Hit, HitRecords, allele_copies, floor_values, obj, opt_int, opt_s, +}; +use super::text::{self, digits, format_map, gt_field, info_map}; +use super::two_hits::{DeletionInputs, Spans, Variant, gene_deletions, two_hits}; +use indexmap::IndexMap; +use pycompat::json::Value; +use pycompat::text::isdigit; +use std::collections::{HashMap, HashSet}; + +/// The genotype check's floors: step 06's GQ and DP, the alt read share a +/// heterozygote must lie within and a homozygous alt call must reach. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct ClinicalFloor { + pub call: Floor, + pub het_ab: (f64, f64), + pub hom_ab: f64, +} + +impl Default for ClinicalFloor { + fn default() -> Self { + ClinicalFloor { + call: Floor::default(), + het_ab: (0.2, 0.8), + hom_ab: 0.8, + } + } +} + +/// The length of the `[0-9.]+` run at the start of `s`. +fn run(s: &str) -> usize { + s.bytes() + .take_while(|c| c.is_ascii_digit() || *c == b'.') + .count() +} + +/// `re.search(KEY([0-9.]+)-([0-9.]+))` (with `pair`) or `KEY([0-9.]+)`: +/// the first place the whole pattern matches. +fn search<'a>(line: &'a str, key: &str, pair: bool) -> Option<(&'a str, Option<&'a str>)> { + line.match_indices(key).find_map(|(i, _)| { + let rest = &line[i + key.len()..]; + let n = run(rest); + if n == 0 { + return None; + } + if !pair { + return Some((&rest[..n], None)); + } + let hi = rest[n..].strip_prefix('-')?; + let m = run(hi); + (m > 0).then(|| (&rest[..n], Some(&hi[..m]))) + }) +} + +fn float(s: &str) -> Result { + crate::float(s).map_err(|e| format!("ValueError: {e}")) +} + +impl ClinicalFloor { + /// Reads the floors `og vcf csq-rank` wrote into the header. + pub fn read(lines: &[String]) -> Result { + let mut f = ClinicalFloor::default(); + for l in lines.iter().take_while(|l| l.starts_with("##")) { + if !l.starts_with("##PGP_CALL_FLOOR=") { + continue; + } + floor_values(l, true, &mut f.call); + if let Some((lo, Some(hi))) = search(l, "HET_AB=", true) { + f.het_ab = (float(lo)?, float(hi)?); + } + if let Some((v, _)) = search(l, "HOM_AB=", false) { + f.hom_ab = float(v)?; + } + } + Ok(f) + } +} + +/// One ranked CSQ entry. +#[derive(Debug, Clone, PartialEq)] +pub struct CsqEntry { + pub impact: String, + pub gene: String, + pub consequence: String, +} + +/// `csq_rank`: the transcript class (MANE, CANONICAL or ANY), its entries, +/// and the frequency word (rare, common, unknown, or none when the CSQ has +/// no frequency field) from MAX_AF, else gnomADe_AF and gnomADg_AF, over +/// every entry. +pub fn csq_rank(fmt: &[String], csq: &str) -> (&'static str, Vec, &'static str) { + let ix: HashMap<&str, usize> = fmt + .iter() + .enumerate() + .map(|(i, k)| (k.as_str(), i)) + .collect(); + let entries: Vec> = csq + .split(',') + .filter(|t| !t.is_empty()) + .map(|t| t.split('|').collect::>()) + .filter(|v| v.len() == fmt.len()) + .collect(); + let get = |v: &[&str], k: &str| ix.get(k).map_or("", |&i| v[i]).to_owned(); + let mane: Vec<&Vec<&str>> = entries + .iter() + .filter(|v| !get(v, "MANE_SELECT").is_empty()) + .collect(); + let canon: Vec<&Vec<&str>> = entries + .iter() + .filter(|v| get(v, "CANONICAL") == "YES") + .collect(); + let (cls, chosen) = if !mane.is_empty() { + ("MANE", mane) + } else if !canon.is_empty() { + ("CANONICAL", canon) + } else { + ("ANY", entries.iter().collect()) + }; + let fields: Vec<&str> = if ix.contains_key("MAX_AF") { + vec!["MAX_AF"] + } else { + ["gnomADe_AF", "gnomADg_AF"] + .into_iter() + .filter(|f| ix.contains_key(f)) + .collect() + }; + let afs: Vec = entries + .iter() + .flat_map(|v| { + fields + .iter() + .flat_map(|f| { + get(v, f) + .split('&') + .filter_map(|a| pycompat::num::parse_float(a).ok()) + .collect::>() + }) + .collect::>() + }) + .collect(); + let freq = if fields.is_empty() { + "none" + } else { + match text::py_max(afs) { + None => "unknown", + Some(m) if m < COMMON_AF => "rare", + Some(_) => "common", + } + }; + let out = chosen + .iter() + .map(|v| CsqEntry { + impact: get(v, "IMPACT"), + gene: get(v, "SYMBOL"), + consequence: get(v, "Consequence").replace('&', ","), + }) + .collect(); + (cls, out, freq) +} + +/// The genotype check of one call: weak, the reasons, GQ, DP and AD. +pub fn call_check( + fmt: &str, + sample: &str, + floor: &ClinicalFloor, +) -> ( + bool, + Vec, + Option, + Option, + Option, +) { + let f = format_map(fmt, sample); + let gq = digits(f.get("GQ").copied()); + let dp = digits(f.get("DP").copied()); + let ad = f.get("AD").copied().unwrap_or(""); + let mut why = Vec::new(); + if let Some(g) = gq + && g < floor.call.gq + { + why.push(format!("GQ {g} < {}", floor.call.gq)); + } + if let Some(d) = dp + && d < floor.call.dp + { + why.push(format!("DP {d} < {}", floor.call.dp)); + } + let gt = f.get("GT").copied().unwrap_or("").replace('|', "/"); + let al: Vec<&str> = gt.split('/').collect(); + let depths: Vec<&str> = ad.split(',').collect(); + if al.len() == 2 + && al.iter().all(|a| isdigit(a)) + && depths.len() > 1 + && depths.iter().all(|x| isdigit(x)) + { + let n: Vec = depths.iter().filter_map(|x| crate::int(x).ok()).collect(); + let tot: i128 = n.iter().sum(); + if tot != 0 { + let ab = n[1..].iter().sum::() as f64 / tot as f64; + let (lo, hi) = floor.het_ab; + if (al[0] == "0") != (al[1] == "0") && !(lo <= ab && ab <= hi) { + why.push(format!( + "alt reads {} of the depth on a heterozygote", + pycompat::fmt::fixed(ab, 2) + )); + } else if al[0] == al[1] && al[1] != "0" && ab < floor.hom_ab { + why.push(format!( + "alt reads {} of the depth on a homozygous call", + pycompat::fmt::fixed(ab, 2) + )); + } + } + } + ( + !why.is_empty(), + why, + gq, + dp, + (!ad.is_empty()).then(|| ad.to_owned()), + ) +} + +/// One record of the clinical VCF with a HIGH-impact ranked entry. +#[derive(Debug, Clone, PartialEq)] +pub struct HighRecord { + pub genotype: String, + pub ps: Option, + /// Gene to its first HIGH consequence. + pub genes: IndexMap, + pub transcript: &'static str, + pub frequency: &'static str, + pub weak: bool, + pub why: Vec, + pub gq: Option, + pub dp: Option, + pub ad: Option, +} + +/// The clinical VCF read once: its records count and its HIGH-impact +/// records. An error is the note of a file that cannot be read. +#[derive(Debug, Clone)] +pub struct ClinicalVcf { + pub records: usize, + /// Keyed by `(chrom, pos, ref, alt)`; an error is the collector's note + /// (a POS that is no integer on a HIGH-impact record, a bad floor). + pub high: Result, String>, +} + +impl ClinicalVcf { + /// One pass over the file, line by line: a whole-genome clinical VCF + /// is never held in memory. + pub fn read(path: &str) -> Result { + let mut records = 0; + let mut scan = HighScan::default(); + text::each_line(path, |l| { + if !l.starts_with('#') { + records += 1; + } + scan.line(l); + }) + .map_err(|e| e.unreadable(path))?; + Ok(ClinicalVcf { + records, + high: scan.finish(), + }) + } +} + +/// `csq_high_impact_genes`, fed one line at a time. The floors come from +/// the header (`clinical_call_floor`), which ends at the first line that +/// does not start with `##`. +#[derive(Default)] +struct HighScan { + header: Vec, + floor: Option>, + fmt: Option>, + out: IndexMap, + err: Option, +} + +impl HighScan { + fn line(&mut self, line: &str) { + if self.floor.is_none() { + if line.starts_with("##") { + self.header.push(line.to_owned()); + } else { + self.floor = Some(ClinicalFloor::read(&std::mem::take(&mut self.header))); + } + } + if self.err.is_some() { + return; + } + if line.starts_with("##INFO= Result, String> { + let floor = self + .floor + .take() + .unwrap_or_else(|| ClinicalFloor::read(&self.header)); + floor?; + match self.err { + Some(e) => Err(e), + None => Ok(self.out), + } + } +} + +/// One data line of `csq_high_impact_genes`. +fn high_record( + fm: &[String], + line: &str, + floor: &ClinicalFloor, + out: &mut IndexMap, +) -> Result<(), String> { + let r: Vec = line.split('\t').map(str::to_owned).collect(); + if r.len() < 8 { + return Ok(()); + } + let (cls, entries, freq) = csq_rank(fm, info_map(&r[7]).get("CSQ").copied().unwrap_or("")); + let mut genes: IndexMap = IndexMap::new(); + for e in entries.into_iter().filter(|e| e.impact == "HIGH") { + genes.entry(e.gene).or_insert(e.consequence); + } + if genes.is_empty() { + return Ok(()); + } + let gt = text::gt_of(&r).to_owned(); + let (weak, why, gq, dp, ad) = if r.len() > 9 { + call_check(&r[8], &r[9], floor) + } else { + (false, Vec::new(), None, None, None) + }; + let pos = text::pos(&r)?; + out.insert( + (r[0].clone(), pos, r[3].clone(), r[4].clone()), + HighRecord { + genotype: gt, + ps: gt_field(&r, "PS"), + genes, + transcript: cls, + frequency: freq, + weak, + why, + gq, + dp, + ad, + }, + ); + Ok(()) +} + +/// A source of a candidate variant. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +enum Kind { + ClinVar, + High, + Sv, +} + +impl Kind { + fn source(self) -> &'static str { + match self { + Kind::ClinVar => "ClinVar", + Kind::High => "HIGH impact", + Kind::Sv => "AnnotSV", + } + } +} + +/// A variant key the verdicts are counted by: a small variant, an SV by +/// its span, or a gene's called deletion. +#[derive(Debug, Clone, PartialEq, Eq, Hash)] +enum AKey { + Small(Key), + Sv(String, String, String, String), + Deletion(String), +} + +/// One candidate: the fields it is printed with (in their order), its key, +/// copies and ploidy, and its call's GT and phase set. +#[derive(Debug, Clone)] +struct Cand { + kind: Kind, + gene: String, + fields: IndexMap, + key: AKey, + copies: usize, + ploidy: usize, + gt: String, + ps: Option, +} + +/// `sf_verdicts`: for each distinct variant of a gene, `None` (reported) or +/// the reason it is not, and the partial-deficiency note of a reported gene +/// with a hypomorph allele. +fn sf_verdicts( + gene: &str, + rule: &Rule, + alleles: &IndexMap, +) -> (HashMap>, Option) { + let sev = |k: &AKey| match k { + AKey::Small(k) => rule.severity.get(k), + _ => None, + }; + let mut out: IndexMap> = IndexMap::new(); + let mut kept: IndexMap = IndexMap::new(); + for (k, cp) in alleles { + match sev(k) { + Some(a) if a.kind.as_deref() == Some("pseudodeficiency") => { + out.insert( + k.clone(), + Some(format!( + "{}: a pseudodeficiency allele, not counted", + a.name + )), + ); + } + _ => { + kept.insert(k.clone(), *cp); + } + } + } + let n: usize = kept.values().map(|(c, _)| c).sum(); + let ri = rule.report_if.as_str(); + let hemi = kept.values().any(|(_, p)| *p == 1); + for (k, (c, p)) in &kept { + let good = match ri { + "any" | "het_female_too" => true, + "biallelic" => n >= 2, + "hemizygous" => hemi || n >= 2, + _ => matches!(k, AKey::Small(s) if rule.variants.contains_key(s)) && *c >= 2 && *p == 2, + }; + let why = if good { + None + } else if ri == "hom_variant" { + Some(format!("{gene} is reported only for: {}", rule.table_text)) + } else if ri == "hemizygous" { + Some(format!( + "one heterozygous variant; {gene} is reported for {}", + rule.table_text + )) + } else { + Some(format!( + "one heterozygous variant; {gene} is reported with {}", + rule.table_text + )) + }; + out.insert(k.clone(), why); + } + let hypo: Vec<&str> = out + .iter() + .filter(|(_, why)| why.is_none()) + .filter_map(|(k, _)| sev(k).filter(|a| a.kind.as_deref() == Some("hypomorph"))) + .map(|a| a.name.as_str()) + .collect(); + let partial = (!hypo.is_empty()).then(|| { + format!( + "partial, not profound: {} is a hypomorph allele, so the deficiency it causes is partial", + hypo.join("; ") + ) + }); + (out.into_iter().collect(), partial) +} + +/// `sf_structural`: step 05's class 4 and 5 SVs in genes whose rule allows +/// them, and the class 4 or 5 deletions over a gene's sv_regions. +fn sf_structural( + inp: &DeletionInputs, + rules: &IndexMap, +) -> (Vec, bool, Option) { + let rows = match &inp.annotsv { + Err(e) => { + return ( + Vec::new(), + false, + Some(format!("the AnnotSV table could not be read: {e}")), + ); + } + Ok(None) => return (Vec::new(), false, None), + Ok(Some(rows)) => rows, + }; + let mut out = Vec::new(); + for x in rows { + if x.common || !matches!(x.class(), "4" | "5") { + continue; + } + let (chrom, span) = x.sv.split_once(':').unwrap_or((&x.sv, "")); + let (a, b) = span.split_once('-').unwrap_or((span, "")); + let chrom = if chrom.starts_with("chr") { + chrom.to_owned() + } else { + format!("chr{chrom}") + }; + let mut hits: IndexMap> = IndexMap::new(); + for g in &x.genes { + if rules.get(g).is_some_and(|r| r.sv_cn) { + hits.insert(g.clone(), None); + } + } + let sv_type = x.sv_type.as_deref().unwrap_or("None"); + let (an, bn) = ( + isdigit(a).then(|| crate::int(a).ok()).flatten(), + isdigit(b).then(|| crate::int(b).ok()).flatten(), + ); + for (g, rule) in rules { + for (c, rs, re, name) in &rule.sv_regions { + if let (Some(an), Some(bn)) = (an, bn) + && sv_type == "DEL" + && *c == chrom + && an <= *re + && bn >= *rs + && !hits.contains_key(g) + { + hits.insert(g.clone(), Some(name.clone())); + } + } + } + hits.sort_keys(); + for (g, via) in hits { + let mut f = IndexMap::new(); + f.insert("gene".into(), Value::Str(g.clone())); + f.insert( + "variant".into(), + Value::Str(format!("{sv_type} {chrom}:{a}-{b}")), + ); + f.insert("acmg_class".into(), Value::Str(x.class().to_owned())); + f.insert("caller".into(), opt_s(x.caller.as_deref())); + f.insert("support".into(), opt_int(x.support)); + f.insert("via".into(), opt_s(via.as_deref())); + out.push(Cand { + kind: Kind::Sv, + gene: g, + fields: f, + key: AKey::Sv( + chrom.clone(), + a.to_owned(), + b.to_owned(), + sv_type.to_owned(), + ), + copies: 1, + ploidy: 2, + gt: String::new(), + ps: None, + }); + } + } + (out, true, None) +} + +fn zyg_map(gt: &str) -> String { + match gt { + "0/1" | "1/0" | "0|1" | "1|0" => "het".into(), + "1/1" | "1|1" => "hom".into(), + _ => gt.into(), + } +} + +/// What the ACMG tier reads. +pub struct Inputs<'a> { + pub assets: &'a str, + /// Step 06's floor and hits, or the note of a hits VCF that cannot be + /// read; `None` without one. + pub hits: Option<&'a Result>, + /// The clinical VCF, or its note; `None` without one. + pub clinical: Option<&'a Result>, + pub deletions: &'a DeletionInputs, +} + +/// `acmg_sf_or_note`: the tier, or `Ok` of the unavailable form when a +/// list or an input cannot be read. An `Err` is the collector's note of the +/// whole clinical section. +pub fn acmg_sf(i: &Inputs) -> Result { + match tier(i)? { + Ok(v) => Ok(v), + Err(why) => Ok(obj([ + ("version", Value::Str(ACMG_SF_VERSION.into())), + ("unavailable", Value::Str(why)), + ("clinvar_hits", Value::Array(Vec::new())), + ("high_impact", Value::Array(Vec::new())), + ("structural", Value::Array(Vec::new())), + ("not_reported", Value::Array(Vec::new())), + ])), + } +} + +/// `acmg_sf_tier`: the outer error is the section's note, the inner one +/// the reason the tier is unavailable. +fn tier(i: &Inputs) -> Result, String> { + let rules = match acmg_sf_rules(i.assets) { + Ok(r) => r, + Err(e) => return Ok(Err(e)), + }; + let mut clinvar: Vec = Vec::new(); + if let Some(h) = i.hits { + let (floor, recs) = match h { + Ok(x) => x, + Err(e) => return Ok(Err(e.clone())), + }; + for r in recs.iter().filter(|r| r.len() >= 8) { + let row = Hit::parse(r, *floor)?; + let gt = text::gt_of(r).to_owned(); + let (copies, ploidy) = allele_copies(>); + for g in row.gene.split(',') { + if !rules.contains_key(g) { + continue; + } + let mut f = IndexMap::new(); + f.insert("gene".into(), Value::Str(g.to_owned())); + f.insert("variant".into(), Value::Str(row.variant())); + f.insert("genotype".into(), Value::Str(row.genotype.clone())); + f.insert("significance".into(), Value::Str(row.significance.clone())); + f.insert("stars".into(), Value::Int(row.stars.into())); + f.insert("low_confidence".into(), Value::Bool(row.low_confidence)); + clinvar.push(Cand { + kind: Kind::ClinVar, + gene: g.to_owned(), + fields: f, + key: AKey::Small(row.key()), + copies, + ploidy, + gt: gt.clone(), + ps: gt_field(r, "PS"), + }); + } + } + } + let mut high: Vec = Vec::new(); + if let Some(c) = i.clinical { + let c = match c { + Ok(c) => c, + Err(e) => return Ok(Err(e.clone())), + }; + let mut recs: Vec<(&Key, &HighRecord)> = + c.high.as_ref().map_err(Clone::clone)?.iter().collect(); + recs.sort_by(|a, b| a.0.cmp(b.0)); + for ((ch, pos, r, a), x) in recs { + if x.frequency == "common" { + continue; + } + let mut genes: Vec<(&String, &String)> = x.genes.iter().collect(); + genes.sort(); + for (g, cons) in genes { + if !rules.contains_key(g) { + continue; + } + let (copies, ploidy) = allele_copies(&x.genotype); + let mut f = IndexMap::new(); + f.insert("gene".into(), Value::Str(g.clone())); + f.insert("variant".into(), Value::Str(format!("{ch}:{pos} {r}>{a}"))); + f.insert("genotype".into(), Value::Str(zyg_map(&x.genotype))); + f.insert("consequence".into(), Value::Str(cons.clone())); + f.insert("transcript".into(), Value::Str(x.transcript.into())); + f.insert("frequency".into(), Value::Str(x.frequency.into())); + f.insert( + "frequency_unknown".into(), + Value::Bool(x.frequency == "unknown"), + ); + f.insert("weak_call".into(), Value::Bool(x.weak)); + f.insert( + "call_note".into(), + if x.why.is_empty() { + Value::Null + } else { + Value::Str(x.why.join("; ")) + }, + ); + f.insert("gq".into(), opt_int(x.gq)); + f.insert("dp".into(), opt_int(x.dp)); + f.insert("ad".into(), opt_s(x.ad.as_deref())); + high.push(Cand { + kind: Kind::High, + gene: g.clone(), + fields: f, + key: AKey::Small((ch.clone(), *pos, r.clone(), a.clone())), + copies, + ploidy, + gt: x.genotype.clone(), + ps: x.ps.clone(), + }); + } + } + } + let (svs, sv_checked, sv_note) = sf_structural(i.deletions, &rules); + + // One verdict per gene, over its distinct variants. + let mut alleles: IndexMap> = IndexMap::new(); + for x in clinvar.iter().chain(&high) { + let p = if x.ploidy == 0 { 2 } else { x.ploidy }; + alleles + .entry(x.gene.clone()) + .or_default() + .entry(x.key.clone()) + .or_insert((x.copies.max(1), p)); + } + for x in &svs { + alleles + .entry(x.gene.clone()) + .or_default() + .entry(x.key.clone()) + .or_insert((1, 2)); + } + let (spans, two_hit_note): (Spans, Option) = match carrier_genes(i.assets) { + Ok(g) => (g.into_iter().map(|(k, v)| (k, v.span)).collect(), None), + Err(e) => ( + IndexMap::new(), + Some(format!( + "A deletion as the second hit: only AnnotSV's gene names were read; step 35's \ + deletions need the gene spans of {e}" + )), + ), + }; + let (dels, _, _) = gene_deletions(i.deletions, &spans); + let mut seen = HashSet::new(); + let mut variants = Vec::new(); + for x in clinvar.iter().chain(&high) { + let AKey::Small(k) = &x.key else { continue }; + let pseudo = rules[&x.gene].severity_kind(k) == Some("pseudodeficiency"); + if !pseudo && seen.insert((x.gene.clone(), k.clone())) { + variants.push(Variant { + gene: x.gene.clone(), + key: k.clone(), + chrom: k.0.clone(), + pos: k.1, + text: match &x.fields["variant"] { + Value::Str(s) => s.clone(), + _ => String::new(), + }, + genotype: match &x.fields["genotype"] { + Value::Str(s) => s.clone(), + _ => String::new(), + }, + copies: x.copies, + ploidy: if x.ploidy == 0 { 2 } else { x.ploidy }, + gt: x.gt.clone(), + ps: x.ps.clone(), + allele: None, + }); + } + } + let found = two_hits(&variants, &dels); + let mut two: HashMap<(String, Key), String> = HashMap::new(); + for g in found.entries.keys() { + if matches!(rules[g].report_if.as_str(), "biallelic" | "hemizygous") { + if let Some(a) = alleles.get_mut(g) { + a.insert(AKey::Deletion(g.clone()), (1, 2)); + } + two.extend( + found + .pairs + .iter() + .filter(|(k, _)| &k.0 == g) + .map(|(k, t)| (k.clone(), t.clone())), + ); + } + } + let verdict: HashMap<&String, _> = alleles + .iter() + .map(|(g, a)| (g, sf_verdicts(g, &rules[g], a))) + .collect(); + let mut reported: [Vec; 3] = Default::default(); + let mut not_reported = Vec::new(); + for x in clinvar.into_iter().chain(high).chain(svs) { + let (whys, partial) = &verdict[&x.gene]; + let why = whys.get(&x.key).cloned().flatten(); + let rule = &rules[&x.gene]; + let mut f = x.fields; + if let AKey::Small(k) = &x.key + && found.hemi.contains_key(&(x.gene.clone(), k.clone())) + { + f.insert("genotype".into(), Value::Str("hemizygous".into())); + } + let slot = match x.kind { + Kind::ClinVar => 0, + Kind::High => 1, + Kind::Sv => 2, + }; + match why { + None => { + f.insert("rule".into(), Value::Str(rule.table_text.clone())); + f.insert("rule_note".into(), opt_s(rule.note.as_deref())); + f.insert("severity".into(), opt_s(partial.as_deref())); + if x.kind != Kind::Sv + && let AKey::Small(k) = &x.key + && let Some(t) = two.get(&(x.gene.clone(), k.clone())) + { + f.insert("two_hit".into(), Value::Str(t.clone())); + } + reported[slot].push(Value::Object(f)); + } + Some(reason) => not_reported.push(obj([ + ("gene", Value::Str(x.gene.clone())), + ("variant", f["variant"].clone()), + ( + "genotype", + f.get("genotype").cloned().unwrap_or(Value::Null), + ), + ("source", Value::Str(x.kind.source().into())), + ("reason", Value::Str(reason)), + ])), + } + } + let [clinvar_rows, high_rows, sv_rows] = reported; + let flag = + |x: &Value, k: &str| matches!(x, Value::Object(m) if m.get(k) == Some(&Value::Bool(true))); + let sure = high_rows + .iter() + .filter(|x| !(flag(x, "weak_call") || flag(x, "frequency_unknown"))) + .count(); + let weak = high_rows.iter().filter(|x| flag(x, "weak_call")).count(); + let unknown = high_rows + .iter() + .filter(|x| flag(x, "frequency_unknown")) + .count(); + Ok(Ok(obj([ + ("version", Value::Str(ACMG_SF_VERSION.into())), + ("genes_on_list", Value::Int(rules.len() as i128)), + ("rules", Value::Str("assets/acmg_sf.tsv".into())), + ("clinvar_hits", Value::Array(clinvar_rows)), + ("clinvar_checked", Value::Bool(i.hits.is_some())), + ("high_impact", Value::Array(high_rows)), + ("high_impact_checked", Value::Bool(i.clinical.is_some())), + ("high_impact_rare", Value::Int(sure as i128)), + ("high_impact_weak", Value::Int(weak as i128)), + ("high_impact_frequency_unknown", Value::Int(unknown as i128)), + ("structural", Value::Array(sv_rows)), + ("structural_checked", Value::Bool(sv_checked)), + ("structural_note", opt_s(sv_note.as_deref())), + ("not_reported", Value::Array(not_reported)), + ("two_hit_note", opt_s(two_hit_note.as_deref())), + ]))) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn fmt(s: &str) -> Vec { + s.split('|').map(str::to_owned).collect() + } + + #[test] + fn mane_then_canonical_then_any() { + let f = fmt("Allele|IMPACT|SYMBOL|Consequence|CANONICAL|MANE_SELECT|gnomADe_AF"); + let (c, e, q) = csq_rank( + &f, + "T|HIGH|A|stop_gained|YES||0.5,T|LOW|A|synonymous_variant||NM_1|0.001", + ); + assert_eq!((c, q), ("MANE", "common")); + assert_eq!(e.len(), 1); + assert_eq!(e[0].impact, "LOW"); + let (c, _, q) = csq_rank(&f, "T|HIGH|A|stop_gained&splice|YES||"); + assert_eq!((c, q), ("CANONICAL", "unknown")); + let f2 = fmt("Allele|IMPACT|SYMBOL|Consequence"); + let (c, e, q) = csq_rank(&f2, "T|HIGH|A|a&b"); + assert_eq!((c, q), ("ANY", "none")); + assert_eq!(e[0].consequence, "a,b"); + } + + #[test] + fn the_genotype_check() { + let f = ClinicalFloor::default(); + let (weak, why, gq, dp, ad) = call_check("GT:GQ:DP:AD", "0/1:20:30:25,5", &f); + assert!(weak); + assert_eq!( + why, + [ + "GQ 20 < 30", + "alt reads 0.17 of the depth on a heterozygote" + ] + ); + assert_eq!((gq, dp, ad.as_deref()), (Some(20), Some(30), Some("25,5"))); + let (weak, why, ..) = call_check("GT:AD", "1|1:5,15", &f); + assert!(weak); + assert_eq!(why, ["alt reads 0.75 of the depth on a homozygous call"]); + let (weak, ..) = call_check("GT:GQ:DP:AD", "0/1:40:30:15,15", &f); + assert!(!weak); + } + + #[test] + fn the_floors_of_the_header() { + let l = |s: &str| vec![s.to_owned(), "#CHROM".to_owned()]; + let f = ClinicalFloor::read(&l( + "##PGP_CALL_FLOOR=GQ=20,DP=5,HET_AB=0.25-0.75,HOM_AB=0.9", + )) + .unwrap(); + assert_eq!(f.call, Floor { gq: 20, dp: 5 }); + assert_eq!((f.het_ab, f.hom_ab), ((0.25, 0.75), 0.9)); + assert!(ClinicalFloor::read(&l("##PGP_CALL_FLOOR=HOM_AB=0..9")).is_err()); + } +} diff --git a/rust/crates/og-summary/src/collect/assets.rs b/rust/crates/og-summary/src/collect/assets.rs new file mode 100644 index 00000000..a71eeb61 --- /dev/null +++ b/rust/crates/og-summary/src/collect/assets.rs @@ -0,0 +1,297 @@ +//! The two gene lists the clinical sections read from `assets/`: +//! `acmg_sf.tsv` (the ACMG SF v3.3 genes and the reporting rule of each, +//! Lee et al., Genet Med 2025, PMID 40568962) and `carrier_genes.tsv` (the +//! genes of the carrier list and their spans). +//! +//! A missing or broken list is an error with the words the reports print; +//! the section that needs it says the list is unavailable and keeps its +//! other values. + +use crate::{basename, join}; +use indexmap::IndexMap; +use pycompat::csv::{DictReader, DictRow, Quoting}; + +/// The list version both reports print. +pub const ACMG_SF_VERSION: &str = "ACMG SF v3.3"; +/// The rules a gene of the ACMG list can have. +pub const SF_REPORT_IF: [&str; 5] = [ + "any", + "biallelic", + "hom_variant", + "hemizygous", + "het_female_too", +]; +/// The kinds of a severity entry. +pub const SEVERITY_KINDS: [&str; 2] = ["hypomorph", "pseudodeficiency"]; + +/// A variant as `(chrom, pos, ref, alt)`. +pub type Key = (String, i128, String, String); + +/// What a severity entry or a variants entry names. +#[derive(Debug, Clone, PartialEq)] +pub struct Allele { + /// `hypomorph` or `pseudodeficiency`; `None` in a variants cell. + pub kind: Option, + pub name: String, +} + +/// The reporting rule of one ACMG gene. +#[derive(Debug, Clone)] +pub struct Rule { + pub report_if: String, + pub variants: IndexMap, + /// AnnotSV's class 4 and 5 SVs count in this gene. + pub sv_cn: bool, + /// Regions where a deletion counts for the gene: chrom, start, end, name. + pub sv_regions: Vec<(String, i128, i128, String)>, + pub severity: IndexMap, + pub table_text: String, + pub note: Option, +} + +impl Rule { + pub fn severity_kind(&self, key: &Key) -> Option<&str> { + self.severity.get(key).and_then(|a| a.kind.as_deref()) + } +} + +/// One gene of the carrier list. +#[derive(Debug, Clone)] +pub struct CarrierGene { + pub inheritance: String, + pub condition: String, + pub severity: IndexMap, + /// GRCh38, 1-based, inclusive. + pub span: (String, i128, i128), +} + +/// The rows of `assets/` without its `#` lines, cells as written +/// (no quoting). +pub fn asset_rows(assets: &str, name: &str) -> Result, String> { + let p = join(assets, name); + let missing = || { + format!( + "{name} is not in {assets}: the report needs the repository's assets/ folder \ + (render_report.py --assets)" + ) + }; + if !crate::is_file(&p) { + return Err(missing()); + } + let bytes = std::fs::read(&p).map_err(|e| crate::os_error(&e, &p))?; + let text = String::from_utf8(bytes).map_err(|e| { + format!( + "UnicodeDecodeError: 'utf-8' codec can't decode {}: {e}", + basename(&p) + ) + })?; + let lines = text + .split_inclusive('\n') + .filter(|l| !l.starts_with('#')) + .map(|l| Ok::<_, std::io::Error>(l.to_owned())); + DictReader::new(lines, '\t', Quoting::None) + .collect::, _>>() + .map_err(|e| format!("{}: {e}", e.exception())) +} + +/// `row[key]`, its error as `str(KeyError)`. +fn col<'a>(r: &'a DictRow, key: &str) -> Result, String> { + r.get(key).ok_or_else(|| crate::repr(key)) +} + +fn int(s: Option<&str>) -> Result { + match s { + None => Err( + "int() argument must be a string, a bytes-like object or a real number, not 'NoneType'" + .into(), + ), + Some(s) => crate::int(s), + } +} + +/// `'chr6:26092913:G:A'` as a key. +pub fn allele_key(locus: &str) -> Result { + let f: Vec<&str> = locus.split(':').collect(); + match f.as_slice() { + [c, pos, r, a] => Ok(( + (*c).to_owned(), + crate::int(pos)?, + (*r).to_owned(), + (*a).to_owned(), + )), + _ if f.len() < 4 => Err(format!( + "not enough values to unpack (expected 4, got {})", + f.len() + )), + _ => Err(format!( + "too many values to unpack (expected 4, got {})", + f.len() + )), + } +} + +/// The variants or severity cell: `.` or entries separated by `;`, each +/// `[kind ]chrom:pos:ref:alt name`. +pub fn parse_alleles(text: Option<&str>, kinds: bool) -> Result, String> { + let mut out = IndexMap::new(); + for e in text.filter(|t| !t.is_empty()).unwrap_or(".").split(';') { + let e = pycompat::text::strip(e, None); + if e.is_empty() || e == "." { + continue; + } + let mut parts: Vec<&str> = e.splitn(if kinds { 3 } else { 2 }, ' ').collect(); + if kinds && parts.len() < 2 { + return Err(format!( + "severity entry needs a kind and a chrom:pos:ref:alt: {e}" + )); + } + let kind = if kinds { + let k = parts.remove(0); + if !SEVERITY_KINDS.contains(&k) { + return Err(format!( + "severity '{k}' is not one of {}: {e}", + SEVERITY_KINDS.join(", ") + )); + } + Some(k.to_owned()) + } else { + None + }; + let name = parts.get(1).copied().unwrap_or(parts[0]).to_owned(); + out.insert(allele_key(parts[0])?, Allele { kind, name }); + } + Ok(out) +} + +/// `{gene: rule}` of `acmg_sf.tsv`, in the file's order. +pub fn acmg_sf_rules(assets: &str) -> Result, String> { + let mut rules: IndexMap = IndexMap::new(); + for (n, r) in asset_rows(assets, "acmg_sf.tsv")?.iter().enumerate() { + let one = || -> Result<(String, Rule), String> { + let g = col(r, "gene")?.unwrap_or("None").to_owned(); + let report_if = col(r, "report_if")?; + let sv_cn = col(r, "sv_cn")?; + let ri_ok = report_if.is_some_and(|x| SF_REPORT_IF.contains(&x)); + if !ri_ok || !matches!(sv_cn, Some("yes" | "no")) { + return Err(format!( + "report_if must be one of {}, sv_cn yes or no", + SF_REPORT_IF.join(", ") + )); + } + let report_if = report_if.unwrap_or_default(); + if sv_cn == Some("yes") && !matches!(report_if, "any" | "het_female_too") { + return Err( + "sv_cn yes needs a rule that one variant meets (any, het_female_too)".into(), + ); + } + let mut regions = Vec::new(); + for e in col(r, "sv_regions")? + .filter(|x| !x.is_empty()) + .unwrap_or(".") + .split(';') + { + let e = pycompat::text::strip(e, None); + if e.is_empty() || e == "." { + continue; + } + let (loc, name) = e.split_once(' ').unwrap_or((e, "")); + let (c, span) = loc.split_once(':').unwrap_or((loc, "")); + let (a, b) = span.split_once('-').unwrap_or((span, "")); + regions.push(( + c.to_owned(), + crate::int(a)?, + crate::int(b)?, + if name.is_empty() { loc } else { name }.to_owned(), + )); + } + let variants = parse_alleles(col(r, "variants")?, false)?; + if report_if == "hom_variant" && variants.is_empty() { + return Err("report_if hom_variant needs its variants".into()); + } + if rules.contains_key(&g) { + return Err("the gene is on two rows".into()); + } + let inheritance = col(r, "inheritance")?; + let _ = inheritance; + let severity = parse_alleles(col(r, "severity")?, true)?; + let table_text = col(r, "table_text")?.unwrap_or("None").to_owned(); + let note = col(r, "note")? + .filter(|x| !matches!(*x, "" | ".")) + .map(str::to_owned); + Ok(( + g, + Rule { + report_if: report_if.to_owned(), + variants, + sv_cn: sv_cn == Some("yes"), + sv_regions: regions, + severity, + table_text, + note, + }, + )) + }; + let (g, rule) = one().map_err(|e| format!("acmg_sf.tsv row {}: {e}", n + 1))?; + rules.insert(g, rule); + } + if rules.is_empty() { + return Err("acmg_sf.tsv has no gene".into()); + } + Ok(rules) +} + +/// `{gene: row}` of `carrier_genes.tsv`; a gene on two rows keeps its +/// first place and its last row. +pub fn carrier_genes(assets: &str) -> Result, String> { + let mut out = IndexMap::new(); + for (n, r) in asset_rows(assets, "carrier_genes.tsv")?.iter().enumerate() { + let one = || -> Result<(String, CarrierGene), String> { + let g = col(r, "gene")?.unwrap_or("None").to_owned(); + let inheritance = col(r, "inheritance")?.unwrap_or("None").to_owned(); + col(r, "lists")?; + let condition = col(r, "condition")?.unwrap_or("None").to_owned(); + let severity = parse_alleles(r.get("severity").flatten(), true)?; + let chrom = col(r, "chrom")?.unwrap_or("None").to_owned(); + let span = (chrom, int(col(r, "start")?)?, int(col(r, "end")?)?); + Ok(( + g, + CarrierGene { + inheritance, + condition, + severity, + span, + }, + )) + }; + let (g, x) = one().map_err(|e| format!("carrier_genes.tsv row {}: {e}", n + 1))?; + out.insert(g, x); + } + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn alleles_read_with_and_without_kinds() { + let v = parse_alleles(Some("chr6:10:G:A C282Y; chr6:20:C:G"), false).unwrap(); + assert_eq!(v.len(), 2); + assert_eq!( + v[&("chr6".into(), 10, "G".into(), "A".into())].name, + "C282Y" + ); + assert_eq!( + v[&("chr6".into(), 20, "C".into(), "G".into())].name, + "chr6:20:C:G" + ); + let s = parse_alleles(Some("hypomorph chr1:5:A:T p.X"), true).unwrap(); + assert_eq!( + s.values().next().unwrap().kind.as_deref(), + Some("hypomorph") + ); + assert!(parse_alleles(Some("odd chr1:5:A:T"), true).is_err()); + assert!(parse_alleles(Some("chr1:5:A"), false).is_err()); + assert!(parse_alleles(None, true).unwrap().is_empty()); + } +} diff --git a/rust/crates/og-summary/src/collect/carrier.rs b/rust/crates/og-summary/src/collect/carrier.rs new file mode 100644 index 00000000..852e35bf --- /dev/null +++ b/rust/crates/og-summary/src/collect/carrier.rs @@ -0,0 +1,146 @@ +//! The sample's own carrier list: step 06's ClinVar hits in the genes of +//! `assets/carrier_genes.tsv`, one person, no partner needed. An allele the +//! list marks hypomorph or pseudodeficiency is set aside, never a carrier +//! finding on its own. The genes of autosomal recessive conditions that a +//! called deletion and a hit both reach are listed as two hits, and a hit +//! called homozygous inside a deletion is shown as hemizygous. + +use super::assets::carrier_genes; +use super::clinvar::{Floor, Hit, HitRecords, allele_copies, carrier_genotype, obj}; +use super::text::{self, gt_field}; +use super::two_hits::{DeletionInputs, Spans, Variant, gene_deletions, two_hits}; +use indexmap::IndexMap; +use pycompat::json::Value; + +/// The section's data: the list, or why it could not be made. +#[derive(Debug, Clone, PartialEq)] +pub struct Carrier(pub Value); + +impl Carrier { + /// `carrier_list_or_note`. An `Err` is the collector's note (a hit + /// whose POS is no integer). + pub fn read( + assets: &str, + hits: &Result, + deletions: &DeletionInputs, + ) -> Result { + let unavailable = |e: String| Carrier(obj([("unavailable", Value::Str(e))])); + let genes = match carrier_genes(assets) { + Ok(g) => g, + Err(e) => return Ok(unavailable(e)), + }; + let recs = match hits { + Ok((_, r)) => r, + Err(e) => return Ok(unavailable(e.clone())), + }; + let mut items: Vec> = Vec::new(); + // Which list each item is on: carriers (false) or set aside (true). + let mut aside: Vec = Vec::new(); + let mut variants: Vec<(Variant, usize)> = Vec::new(); + for r in recs.iter().filter(|r| r.len() >= 8) { + let row = Hit::parse(r, Floor::default())?; + let gt = text::gt_of(r).to_owned(); + let (copies, ploidy) = allele_copies(>); + let key = row.key(); + for g in row.gene.split(',') { + let Some(x) = genes.get(g) else { continue }; + let mut item = IndexMap::new(); + item.insert("gene".into(), Value::Str(g.to_owned())); + item.insert("variant".into(), Value::Str(row.variant())); + item.insert("genotype".into(), Value::Str(carrier_genotype(>, &row))); + item.insert("significance".into(), Value::Str(row.significance.clone())); + item.insert("stars".into(), Value::Int(row.stars.into())); + item.insert("inheritance".into(), Value::Str(x.inheritance.clone())); + item.insert("condition".into(), Value::Str(x.condition.clone())); + let sev = x.severity.get(&key); + if let Some(s) = sev { + item.insert( + "kind".into(), + s.kind.clone().map_or(Value::Null, Value::Str), + ); + item.insert("name".into(), Value::Str(s.name.clone())); + } + let genotype = match &item["genotype"] { + Value::Str(s) => s.clone(), + _ => String::new(), + }; + items.push(item); + aside.push(sev.is_some()); + if sev.is_none_or(|s| s.kind.as_deref() == Some("hypomorph")) { + variants.push(( + Variant { + gene: g.to_owned(), + key: key.clone(), + chrom: row.chrom.clone(), + pos: row.pos, + text: row.variant(), + genotype, + copies, + ploidy, + gt: gt.clone(), + ps: gt_field(r, "PS"), + allele: sev.map(|s| s.name.clone()), + }, + items.len() - 1, + )); + } + } + } + let spans: Spans = genes + .iter() + .map(|(g, x)| (g.clone(), x.span.clone())) + .collect(); + let (dels, read, unread) = gene_deletions(deletions, &spans); + let vs: Vec = variants.iter().map(|(v, _)| v.clone()).collect(); + let found = two_hits(&vs, &dels); + for (v, at) in &variants { + if let Some(t) = found.hemi.get(&(v.gene.clone(), v.key.clone())) { + let item = &mut items[*at]; + item.insert("genotype".into(), Value::Str("hemizygous".into())); + item.insert( + "hemizygous".into(), + Value::Str(format!("called homozygous inside {t}")), + ); + } + } + let mut rows: Vec<(&String, _)> = found.entries.iter().collect(); + rows.sort_by(|a, b| a.0.cmp(b.0)); + let two: Vec = rows + .into_iter() + .filter(|(g, _)| genes[*g].inheritance == "AR") + .map(|(g, e)| { + let mut m = e.to_json(); + m.insert("condition".into(), Value::Str(genes[g].condition.clone())); + m.insert("inheritance".into(), Value::Str("AR".into())); + Value::Object(m) + }) + .collect(); + let (mut carriers, mut set_aside) = (Vec::new(), Vec::new()); + for (item, a) in items.into_iter().zip(aside) { + if a { + set_aside.push(Value::Object(item)); + } else { + carriers.push(Value::Object(item)); + } + } + let strs = |v: Vec| Value::Array(v.into_iter().map(Value::Str).collect()); + Ok(Carrier(obj([ + ("genes_on_list", Value::Int(genes.len() as i128)), + ("list", Value::Str("assets/carrier_genes.tsv".into())), + ("carriers", Value::Array(carriers)), + ("set_aside", Value::Array(set_aside)), + ( + "two_hits", + obj([ + ("genes", Value::Array(two)), + ("deletions_read", strs(read)), + ("deletions_not_read", strs(unread)), + ]), + ), + ]))) + } + + pub fn to_json(&self) -> Value { + self.0.clone() + } +} diff --git a/rust/crates/og-summary/src/collect/clinical.rs b/rust/crates/og-summary/src/collect/clinical.rs new file mode 100644 index 00000000..5936475e --- /dev/null +++ b/rust/crates/og-summary/src/collect/clinical.rs @@ -0,0 +1,210 @@ +//! The clinical variant filter (step 23) with its ACMG secondary-findings +//! tier, and the variant prioritization (step 31, slivar) with its +//! compound-het pairs phased from the reads. +//! +//! The bash pipeline writes summary tables; the Nextflow report gets the +//! VCFs instead, so each section reads the table when there is one and +//! counts the VCF's records otherwise. + +use super::acmg::{self, ClinicalVcf}; +use super::clinvar::obj; +use super::text::{self, cell, digits}; +use indexmap::IndexMap; +use pycompat::csv::DictRow; +use pycompat::json::Value; +use std::collections::BTreeSet; + +/// The trans pairs the report lists by name; the rest are counted. +pub const COMPHET_TRANS_SHOWN: usize = 20; + +fn rows_or_note(path: &str) -> Result, String> { + text::rows(path).map_err(|e| e.note()) +} + +fn int(n: usize) -> Value { + Value::Int(n as i128) +} + +/// The data of the clinical section. +#[derive(Debug, Clone, PartialEq)] +pub struct Clinical(pub Value); + +impl Clinical { + /// `sec_clinical`: `table` is step 23's summary table, `vcf` the + /// clinical VCF read once (both full paths). + pub fn read( + table: Option<&str>, + vcf: Option<&Result>, + tier: &acmg::Inputs, + ) -> Result { + let Some(table) = table else { + // The Nextflow report: the record count only (one record per + // variant, as in the table). + let v = match vcf { + Some(Ok(v)) => v, + Some(Err(e)) => return Err(e.clone()), + None => return Ok(Clinical(Value::Null)), + }; + return Ok(Clinical(obj([ + ("variants", int(v.records)), + ("acmg_sf", acmg::acmg_sf(tier)?), + ]))); + }; + let rows = rows_or_note(table)?; + let mut by: IndexMap = IndexMap::new(); + for r in &rows { + let k = cell(r, "IMPACT").filter(|v| !v.is_empty()).unwrap_or("."); + *by.entry(k.to_owned()).or_default() += 1; + } + let genes: BTreeSet<&str> = rows + .iter() + .filter_map(|r| cell(r, "GENE")) + .filter(|g| !matches!(*g, "" | ".")) + .collect(); + let mut m: IndexMap = IndexMap::new(); + m.insert("variants".into(), int(rows.len())); + m.insert( + "by_impact".into(), + Value::Object(by.into_iter().map(|(k, v)| (k, int(v))).collect()), + ); + m.insert("genes".into(), int(genes.len())); + m.insert("acmg_sf".into(), acmg::acmg_sf(tier)?); + // Step 23's tier and call check columns; an older table has neither. + if let Some(r) = rows.first() { + if r.get("FREQ_TIER").is_some() { + let n = rows + .iter() + .filter(|r| cell(r, "FREQ_TIER") == Some("unknown")) + .count(); + m.insert("frequency_unknown".into(), int(n)); + } + if r.get("CALL_CHECK").is_some() { + let n = rows + .iter() + .filter(|r| cell(r, "CALL_CHECK") == Some("weak")) + .count(); + m.insert("weak_calls".into(), int(n)); + } + } + Ok(Clinical(Value::Object(m))) + } + + pub fn to_json(&self) -> Value { + self.0.clone() + } +} + +/// The data of the slivar section. +#[derive(Debug, Clone, PartialEq)] +pub struct Slivar(pub Value); + +impl Slivar { + /// `sec_slivar` with step 31's summary table (`table`) or, for the + /// Nextflow report, its prioritized VCF (`vcf`); full paths. + pub fn read(table: Option<&str>, vcf: Option<&str>, s: &str) -> Result { + let mut m: IndexMap = IndexMap::new(); + let Some(p) = table else { + let Some(v) = vcf else { + return Ok(Slivar(Value::Null)); + }; + let mut n = 0; + text::each_line(v, |l| n += usize::from(!l.starts_with('#'))) + .map_err(|e| e.unreadable(v))?; + m.insert("prioritized".into(), int(n)); + comphet_phase( + &crate::join(text::dirname(v), &format!("{s}_compound_het_pairs.tsv")), + &mut m, + )?; + return Ok(Slivar(Value::Object(m))); + }; + let rows = rows_or_note(p)?; + m.insert("prioritized".into(), int(rows.len())); + // Step 31's frequency tier: prioritized for the compound-het search, not rare. + if rows.first().is_some_and(|r| r.get("FREQ_TIER").is_some()) { + let n = rows + .iter() + .filter(|r| cell(r, "FREQ_TIER") == Some("unknown")) + .count(); + m.insert("frequency_unknown".into(), int(n)); + } + let dir = text::dirname(p); + let ch = crate::join(dir, &format!("{s}_compound_hets.tsv")); + if crate::is_file(&ch) { + let c = rows_or_note(&ch)?; + let genes: std::collections::HashSet<&str> = c + .iter() + .filter_map(|r| cell(r, "GENE")) + .filter(|g| !g.is_empty()) + .collect(); + m.insert("compound_het_variants".into(), int(c.len())); + m.insert("compound_het_genes".into(), int(genes.len())); + } + comphet_phase( + &crate::join(dir, &format!("{s}_compound_het_pairs.tsv")), + &mut m, + )?; + Ok(Slivar(Value::Object(m))) + } + + pub fn to_json(&self) -> Value { + self.0.clone() + } +} + +/// Step 31's pairs table: the pairs by phase from the reads, and the trans +/// ones (both copies of the gene hit) by name. A run from before the table +/// adds none of these keys. +fn comphet_phase(path: &str, m: &mut IndexMap) -> Result<(), String> { + if !crate::is_file(path) { + return Ok(()); + } + let rows = rows_or_note(path)?; + let phase = obj(["trans", "cis", "unknown"].map(|k| { + ( + k, + int(rows.iter().filter(|r| cell(r, "PHASE") == Some(k)).count()), + ) + })); + let py = |r: &DictRow, k: &str| cell(r, k).unwrap_or("None").to_owned(); + let n = |r: &DictRow, k: &str| digits(cell(r, k)).unwrap_or(0); + let trans: Vec = rows + .iter() + .filter(|r| cell(r, "PHASE") == Some("trans")) + .take(COMPHET_TRANS_SHOWN) + .map(|r| { + let gene = match r.get("GENE") { + None => Value::Str(String::new()), + Some(v) => v.map_or(Value::Null, |g| Value::Str(g.to_owned())), + }; + obj([ + ("gene", gene), + ( + "variants", + Value::Str(format!( + "{}:{} {}>{} + {}:{} {}>{}", + py(r, "CHROM"), + py(r, "POS_1"), + py(r, "REF_1"), + py(r, "ALT_1"), + py(r, "CHROM"), + py(r, "POS_2"), + py(r, "REF_2"), + py(r, "ALT_2") + )), + ), + ( + "reads", + Value::Str(format!( + "{} of {}", + py(r, "TRANS_READS"), + n(r, "TRANS_READS") + n(r, "CIS_READS") + )), + ), + ]) + }) + .collect(); + m.insert("compound_het_pairs".into(), int(rows.len())); + m.insert("compound_het_phase".into(), phase); + m.insert("compound_het_trans".into(), Value::Array(trans)); + Ok(()) +} diff --git a/rust/crates/og-summary/src/collect/clinvar.rs b/rust/crates/og-summary/src/collect/clinvar.rs new file mode 100644 index 00000000..f930c25e --- /dev/null +++ b/rust/crates/og-summary/src/collect/clinvar.rs @@ -0,0 +1,800 @@ +//! The ClinVar screen (step 06): each P/LP hit with its review stars, its +//! call quality against step 06's floor, gnomAD AF and rsIDs from step 13's +//! VEP output, and the cDNA read check, the cDNA scan and the always-report +//! deep intronic sites when step 06 wrote them. +//! +//! The hits VCF is small and three sections read it (this one, the ACMG +//! tier of the clinical section and the carrier list), so it is read once +//! ([`HitsVcf`]) and each section parses the rows it needs. + +use super::text::{self, ReadError, digits, format_map, gt_of, info_map, py_max}; +use indexmap::IndexMap; +use og_clinvar::{ReviewStatus, combined_note, half_call, zygosity}; +use pycompat::json::Value; +use std::collections::HashMap; + +/// gnomAD AF above which a P/LP hit is marked common. +pub const COMMON_AF: f64 = 0.01; + +/// Step 06's call floor: a hit whose GQ or DP is below it is low confidence. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Floor { + pub gq: i128, + pub dp: i128, +} + +impl Default for Floor { + fn default() -> Self { + Floor { gq: 30, dp: 10 } + } +} + +impl Floor { + fn to_json(self) -> Value { + obj([("GQ", Value::Int(self.gq)), ("DP", Value::Int(self.dp))]) + } +} + +/// `(GQ|DP)=(\d+)` over a header line: every match, the last one wins. +pub fn floor_values(line: &str, word_start: bool, floor: &mut Floor) { + let b = line.as_bytes(); + for key in ["GQ=", "DP="] { + let mut from = 0; + while let Some(i) = line[from..].find(key).map(|i| i + from) { + from = i + 1; + if word_start && i > 0 && (b[i - 1].is_ascii_alphanumeric() || b[i - 1] == b'_') { + continue; + } + let rest = &line[i + 3..]; + let n = rest.bytes().take_while(u8::is_ascii_digit).count(); + if n > 0 + && let Ok(v) = crate::int(&rest[..n]) + { + // Python takes the matches in text order; with one key per + // pass, the last of each key is the one that stays. + if key == "GQ=" { + floor.gq = v; + } else { + floor.dp = v; + } + } + } + } +} + +/// A variant as `(chrom, pos, ref, alt)`. +pub type VarKey = (String, i128, String, String); +/// The hits VCF's call floor and data lines. +pub type HitRecords = (Floor, Vec>); + +/// The hits VCF, read once: its call floor and its data lines. +#[derive(Debug, Clone)] +pub struct HitsVcf { + /// The path (full) and the path relative to the sample folder. + pub path: String, + pub rel: String, + /// The floor and the records, or the note of a file that cannot be read. + pub read: Result, +} + +impl HitsVcf { + pub fn read(path: String, rel: String) -> HitsVcf { + let mut floor = Floor::default(); + let mut header = true; + let mut recs: Vec> = Vec::new(); + let read = text::each_line(&path, |l| { + header &= l.starts_with("##"); + if header && l.starts_with("##PGP_CALL_FLOOR=") { + floor_values(l, false, &mut floor); + } + if !l.starts_with('#') { + recs.push(l.split('\t').map(str::to_owned).collect()); + } + }) + .map(|()| (floor, recs)) + .map_err(|e: ReadError| e.unreadable(&path)); + HitsVcf { path, rel, read } + } +} + +/// One hit as `clinvar_row` reads it. +#[derive(Debug, Clone, PartialEq)] +pub struct Hit { + pub stars: u8, + pub chrom: String, + pub pos: i128, + pub reference: String, + pub alt: String, + pub genotype: String, + pub gene: String, + pub significance: String, + pub review_status: String, + pub combined_allele: String, + pub dp: Option, + pub gq: Option, + pub low_confidence: bool, + pub frequency: Option, + pub cdna: Option, +} + +/// gnomAD and rsIDs of one hit from VEP. +#[derive(Debug, Clone, PartialEq)] +pub struct Frequency { + pub exome: Option, + pub genome: Option, + pub rsids: Vec, +} + +impl Hit { + /// `clinvar_row`: an error is the collector's note (a POS that is no + /// integer). + pub fn parse(r: &[String], floor: Floor) -> Result { + let info = if r.len() > 7 { + info_map(&r[7]) + } else { + HashMap::new() + }; + let get = |k: &str| info.get(k).copied().unwrap_or(""); + let rev = get("CLNREVSTAT"); + let genes: Vec<&str> = get("GENEINFO") + .split('|') + .filter(|g| !g.is_empty()) + .map(|g| g.split(':').next().unwrap_or("")) + .collect(); + let gt = gt_of(r); + let fmt = if r.len() > 9 { + format_map(&r[8], &r[9]) + } else { + HashMap::new() + }; + let dp = digits(fmt.get("DP").copied()); + let gq = digits(fmt.get("GQ").copied()); + let pos = text::pos(r)?; + let s = |v: &str| v.replace('_', " "); + let or_dot = |v: String| if v.is_empty() { ".".to_owned() } else { v }; + Ok(Hit { + stars: ReviewStatus::parse(rev.as_bytes()).stars(), + chrom: r[0].clone(), + pos, + reference: r[3].clone(), + alt: r[4].clone(), + genotype: String::from_utf8_lossy(&zygosity(gt.as_bytes())).into_owned(), + gene: or_dot(genes.join(",")), + significance: or_dot(s(get("CLNSIG"))), + review_status: or_dot(s(rev)), + combined_allele: String::from_utf8_lossy(&combined_note( + get("PGP_JOINED_FROM").as_bytes(), + get("PGP_PART_OF").as_bytes(), + get("PGP_ONE_ALLELE").as_bytes(), + gt.as_bytes(), + )) + .into_owned(), + dp, + gq, + low_confidence: gq.is_some_and(|v| v < floor.gq) || dp.is_some_and(|v| v < floor.dp), + frequency: None, + cdna: None, + }) + } + + pub fn key(&self) -> (String, i128, String, String) { + ( + self.chrom.clone(), + self.pos, + self.reference.clone(), + self.alt.clone(), + ) + } + + /// `chrom:pos ref>alt`. + pub fn variant(&self) -> String { + format!( + "{}:{} {}>{}", + self.chrom, self.pos, self.reference, self.alt + ) + } + + fn to_json(&self) -> Value { + let mut m = IndexMap::new(); + let put = |m: &mut IndexMap, k: &str, v: Value| { + m.insert(k.to_owned(), v); + }; + put(&mut m, "stars", Value::Int(self.stars.into())); + put(&mut m, "chrom", Value::Str(self.chrom.clone())); + put(&mut m, "pos", Value::Int(self.pos)); + put(&mut m, "ref", Value::Str(self.reference.clone())); + put(&mut m, "alt", Value::Str(self.alt.clone())); + put(&mut m, "genotype", Value::Str(self.genotype.clone())); + put(&mut m, "gene", Value::Str(self.gene.clone())); + put( + &mut m, + "significance", + Value::Str(self.significance.clone()), + ); + put( + &mut m, + "review_status", + Value::Str(self.review_status.clone()), + ); + put( + &mut m, + "combined_allele", + Value::Str(self.combined_allele.clone()), + ); + put(&mut m, "dp", opt_int(self.dp)); + put(&mut m, "gq", opt_int(self.gq)); + put(&mut m, "low_confidence", Value::Bool(self.low_confidence)); + if let Some(f) = &self.frequency { + let af = py_max(f.exome.into_iter().chain(f.genome)); + put(&mut m, "gnomad_exome_af", opt_f(f.exome)); + put(&mut m, "gnomad_genome_af", opt_f(f.genome)); + put( + &mut m, + "rsids", + Value::Array(f.rsids.iter().cloned().map(Value::Str).collect()), + ); + put(&mut m, "af", opt_f(af)); + put( + &mut m, + "common", + Value::Bool(af.is_some_and(|a| a > COMMON_AF)), + ); + } + if let Some(c) = &self.cdna { + put(&mut m, "cdna", c.clone()); + } + Value::Object(m) + } +} + +/// A hit's genotype as the carrier list words it. +pub fn carrier_genotype(gt: &str, hit: &Hit) -> String { + if half_call(gt.as_bytes()) { + hit.genotype.clone() + } else { + let (copies, ploidy) = allele_copies(gt); + zygosity_word(copies, ploidy).to_owned() + } +} + +/// `(alt copies, ploidy)` of a GT: a missing allele is a copy not called, +/// not a missing copy, so the half of a split call is never hemizygous. +pub fn allele_copies(gt: &str) -> (usize, usize) { + let a: Vec<&str> = gt.split(['/', '|']).filter(|x| !x.is_empty()).collect(); + ( + a.iter().filter(|x| !matches!(**x, "0" | ".")).count(), + a.len(), + ) +} + +pub fn zygosity_word(copies: usize, ploidy: usize) -> &'static str { + if copies == 0 { + "unknown" + } else if ploidy == 1 { + "hemizygous" + } else if copies >= 2 { + "hom" + } else { + "het" + } +} + +pub fn obj<'a>(pairs: impl IntoIterator) -> Value { + Value::Object(pairs.into_iter().map(|(k, v)| (k.to_owned(), v)).collect()) +} + +pub fn opt_int(x: Option) -> Value { + x.map_or(Value::Null, Value::Int) +} + +pub fn opt_f(x: Option) -> Value { + x.map_or(Value::Null, Value::Float) +} + +pub fn opt_s(x: Option<&str>) -> Value { + x.map_or(Value::Null, |s| Value::Str(s.to_owned())) +} + +// --- the cDNA read check, the scan, the deep intronic sites ---------------------- + +/// A row key of a cDNA check table: chrom and the alleles may be missing. +pub type CdnaKey = (Option, i128, Option, Option); + +/// `clinvar_cdna_checks`: the read check of each row, keyed by variant. +pub fn cdna_checks(path: &str) -> Result, String> { + let rows = text::rows(path).map_err(|e| match e { + ReadError::Os(_) => e.unreadable(path), + ReadError::Other(n) => n, + })?; + let mut out = HashMap::new(); + for r in &rows { + let Some(pos) = digits(text::cell(r, "pos")) else { + continue; + }; + let c = |k: &str| text::cell(r, k).map(str::to_owned); + let or = |k: &str, d: &str| { + Value::Str( + text::cell(r, k) + .filter(|v| !v.is_empty()) + .unwrap_or(d) + .to_owned(), + ) + }; + let n = |k: &str| opt_int(digits(text::cell(r, k))); + out.insert( + (c("chrom"), pos, c("ref"), c("alt")), + obj([ + ("check", or("cdna_check", "not_checked")), + ("note", or("cdna_note", "")), + ("alt_reads", n("alt_reads")), + ("alt_clipped", n("alt_clipped")), + ("ref_reads", n("ref_reads")), + ("ref_clipped", n("ref_clipped")), + ( + "construct", + opt_s(text::cell(r, "construct_site").filter(|v| !v.is_empty())), + ), + ]), + ); + } + Ok(out) +} + +/// The check of a variant the table does not hold. +pub fn not_in_table() -> Value { + obj([ + ("check", Value::Str("not_checked".into())), + ("note", Value::Str("not in the read-check table".into())), + ("alt_reads", Value::Null), + ("alt_clipped", Value::Null), + ("ref_reads", Value::Null), + ("ref_clipped", Value::Null), + ]) +} + +pub fn cdna_key(chrom: &str, pos: i128, r: &str, a: &str) -> CdnaKey { + ( + Some(chrom.to_owned()), + pos, + Some(r.to_owned()), + Some(a.to_owned()), + ) +} + +pub fn check_word(c: &Value) -> Option<&str> { + match c { + Value::Object(m) => match m.get("check") { + Some(Value::Str(s)) => Some(s), + _ => None, + }, + _ => None, + } +} + +fn has_construct(c: &Value) -> bool { + matches!(c, Value::Object(m) if m.get("construct").is_some_and(pycompat::value::truthy)) +} + +/// `clinvar_cdna_scan`: the scan's note and loci. +fn cdna_scan(path: &str) -> Result { + let lines = text::lines(path).map_err(|e| e.unreadable(path))?; + let note = lines + .first() + .and_then(|l| l.strip_prefix("# cdna_scan: ")) + .unwrap_or("") + .to_owned(); + let body = lines + .iter() + .filter(|l| !l.starts_with('#')) + .map(|l| Ok::<_, std::io::Error>(l.clone())); + let mut loci = Vec::new(); + for r in pycompat::csv::DictReader::new(body, '\t', pycompat::csv::Quoting::Minimal) { + let r = r.map_err(|e| format!("{}: {e}", e.exception()))?; + let Some(start) = digits(text::cell(&r, "start")) else { + continue; + }; + let int0 = |k: &str| -> Result { + match text::cell(&r, k).filter(|v| !v.is_empty()) { + None => Ok(Value::Int(0)), + Some(v) => crate::int(v) + .map(Value::Int) + .map_err(|e| format!("ValueError: {e}")), + } + }; + loci.push(obj([ + ("chrom", opt_s(text::cell(&r, "chrom"))), + ("start", Value::Int(start)), + ("end", int0("end")?), + ("junctions", int0("junctions")?), + ("reads", int0("reads")?), + ( + "junction_list", + Value::Str(text::cell_or(&r, "junction_list").to_owned()), + ), + ])); + } + Ok(obj([ + ("note", Value::Str(note)), + ("loci", Value::Array(loci)), + ])) +} + +/// `clinvar_deep_intronic`: one row per always-report site. +fn deep_intronic(path: &str) -> Result { + let rows = text::rows(path).map_err(|e| match e { + ReadError::Os(_) => e.unreadable(path), + ReadError::Other(n) => n, + })?; + let mut out = Vec::new(); + for r in &rows { + let Some(pos) = digits(text::cell(r, "pos")) else { + continue; + }; + let or = |k: &str, d: &str| { + Value::Str( + text::cell(r, k) + .filter(|v| !v.is_empty()) + .unwrap_or(d) + .to_owned(), + ) + }; + out.push(obj([ + ("chrom", opt_s(text::cell(r, "chrom"))), + ("pos", Value::Int(pos)), + ("ref", opt_s(text::cell(r, "ref"))), + ("alt", opt_s(text::cell(r, "alt"))), + ("gene", or("gene", ".")), + ("hgvs", or("hgvs", "")), + ("clinvar_id", or("clinvar_id", "")), + ("call", or("call", "no call")), + ("detail", or("detail", "")), + ])); + } + Ok(Value::Array(out)) +} + +// --- gnomAD AF and rsIDs from step 13's VEP output ----------------------------------- + +/// Where the frequencies came from, or why there are none. +#[derive(Debug, Clone, PartialEq)] +pub struct FrequencySource { + pub source: Option, + pub note: Option, +} + +/// The VEP output's frequencies of every hit, in one pass per contig over +/// the hits sorted by position: one indexed query per hit position, served +/// by the reader's block cache when hits share a block. +pub fn vep_frequencies(vep: &str, keys: &[VarKey]) -> Result, String> { + let mut fmt: Option> = None; + { + let mut lines = og_io::Lines::new( + og_io::open(vep).map_err(|e| crate::os_error(&e, vep))?, + og_io::Utf8::Replace, + ); + let mut line = String::new(); + while lines.read_line(&mut line).map_err(|e| e.to_string())? { + if !line.starts_with('#') { + break; + } + if line.starts_with("##INFO= = keys.iter().collect(); + sorted.sort(); + sorted.dedup(); + let mut out = HashMap::new(); + if sorted.is_empty() { + return Ok(out); + } + let mut reader = + og_io::tabix::Reader::with_index(vep, format!("{vep}.tbi")).map_err(|e| e.to_string())?; + for key in sorted { + let (chrom, pos, reference, alt) = key; + let Ok(p) = u64::try_from(*pos) else { continue }; + if p == 0 { + continue; + } + let lines: Vec> = reader + .query(chrom, p - 1, p) + .collect::>() + .map_err(|e| e.to_string())?; + for raw in lines { + let line = String::from_utf8_lossy(&raw); + let r: Vec<&str> = line.split('\t').collect(); + if r.len() < 3 || r[1].parse::().ok() != Some(*pos) { + continue; + } + let alts: Vec<&str> = if r.len() > 7 { + r[4].split(',').collect() + } else { + Vec::new() + }; + if r.get(3) != Some(&reference.as_str()) || !alts.contains(&alt.as_str()) { + continue; + } + // VEP's Allele drops the base an indel's alleles share. + let first = |s: &str| s.chars().next(); + let snv = reference.chars().count() == 1 && alts.iter().all(|a| a.chars().count() == 1); + let allele = if snv || alts.iter().any(|a| first(a) != first(reference)) { + alt.clone() + } else { + let rest: String = alt.chars().skip(1).collect(); + if rest.is_empty() { + "-".to_owned() + } else { + rest + } + }; + let mut got = [None, None]; + let mut ids: Vec = Vec::new(); + for tr in info_map(r[7]).get("CSQ").copied().unwrap_or("").split(',') { + let v: Vec<&str> = tr.split('|').collect(); + if v.len() != fmt.len() || (alts.len() > 1 && v[allele_ix] != allele) { + continue; + } + for (i, c) in cols.iter().enumerate() { + if let Some(c) = c + && got[i].is_none() + { + got[i] = py_max( + v[*c] + .split('&') + .filter_map(|x| pycompat::num::parse_float(x).ok()), + ); + } + } + if let Some(e) = ev { + for x in v[e].split('&') { + if is_rsid(x) { + let x = x.to_lowercase(); + if !ids.contains(&x) { + ids.push(x); + } + } + } + } + } + ids.sort_by_key(|x| x[2..].parse::().unwrap_or(u128::MAX)); + out.insert( + key.clone(), + Frequency { + exome: got[0], + genome: got[1], + rsids: ids, + }, + ); + } + } + Ok(out) +} + +/// `(?i)rs\d+`, the whole field. +fn is_rsid(x: &str) -> bool { + let b = x.as_bytes(); + b.len() > 2 && b[..2].eq_ignore_ascii_case(b"rs") && b[2..].iter().all(u8::is_ascii_digit) +} + +// --- the section -------------------------------------------------------------------- + +/// The data of the ClinVar section. +#[derive(Debug, Clone, PartialEq)] +pub struct ClinVar { + pub hits: Vec, + pub floor: Floor, + pub frequency: FrequencySource, + pub cdna: bool, + pub cdna_scan: Option, + pub deep_intronic: Option, +} + +impl ClinVar { + /// `sec_clinvar` over the hits read once. `side(rel)` is the file of + /// another of step 06's or 13's outputs, found in the folder. + pub fn read( + hits: &HitsVcf, + dir: &str, + side: &dyn Fn(&str) -> Option<(String, String)>, + s: &str, + ) -> Result { + let (floor, recs) = match &hits.read { + Ok((f, r)) => (*f, r), + Err(e) => return Err(e.clone()), + }; + let mut rows = Vec::new(); + for r in recs.iter().filter(|r| r.len() >= 8) { + rows.push(Hit::parse(r, floor)?); + } + rows.sort_by(|a, b| { + b.stars + .cmp(&a.stars) + .then_with(|| a.chrom.cmp(&b.chrom)) + .then_with(|| a.pos.cmp(&b.pos)) + }); + let frequency = frequencies(dir, side, s, &mut rows); + let mut cdna = false; + if let Some((p, _)) = side(&format!("clinvar/{s}_clinvar_cdna_check.tsv")) { + let table = cdna_checks(&p)?; + for x in &mut rows { + x.cdna = Some( + table + .get(&cdna_key(&x.chrom, x.pos, &x.reference, &x.alt)) + .cloned() + .unwrap_or_else(not_in_table), + ); + } + cdna = true; + } + let cdna_scan = match side(&format!("clinvar/{s}_clinvar_cdna_scan.tsv")) { + Some((p, _)) => Some(cdna_scan(&p)?), + None => None, + }; + let deep_intronic = match side(&format!("clinvar/{s}_clinvar_deep_intronic.tsv")) { + Some((p, _)) => Some(deep_intronic(&p)?), + None => None, + }; + Ok(ClinVar { + hits: rows, + floor, + frequency, + cdna, + cdna_scan, + deep_intronic, + }) + } + + pub fn to_json(&self) -> Value { + let sure: Vec<&Hit> = self.hits.iter().filter(|x| !x.low_confidence).collect(); + let by = (0..=4u8) + .rev() + .map(|k| { + ( + k.to_string(), + Value::Int(sure.iter().filter(|x| x.stars == k).count() as i128), + ) + }) + .collect(); + let mut m: IndexMap = IndexMap::new(); + m.insert("count".into(), Value::Int(sure.len() as i128)); + m.insert("by_stars".into(), Value::Object(by)); + m.insert( + "hits".into(), + Value::Array(self.hits.iter().map(Hit::to_json).collect()), + ); + m.insert( + "low_confidence".into(), + Value::Int((self.hits.len() - sure.len()) as i128), + ); + m.insert("call_floor".into(), self.floor.to_json()); + m.insert( + "frequency".into(), + obj([ + ("source", opt_s(self.frequency.source.as_deref())), + ("note", opt_s(self.frequency.note.as_deref())), + ("common_af", Value::Float(COMMON_AF)), + ]), + ); + if self.cdna { + let checks: Vec<&Value> = self.hits.iter().filter_map(|x| x.cdna.as_ref()).collect(); + m.insert( + "cdna_flagged".into(), + Value::Int( + checks + .iter() + .filter(|c| check_word(c) == Some("flagged")) + .count() as i128, + ), + ); + m.insert( + "construct_sites".into(), + Value::Int(checks.iter().filter(|c| has_construct(c)).count() as i128), + ); + } + if let Some(v) = &self.cdna_scan { + m.insert("cdna_scan".into(), v.clone()); + } + if let Some(v) = &self.deep_intronic { + m.insert("deep_intronic".into(), v.clone()); + } + Value::Object(m) + } +} + +/// `clinvar_frequencies`: puts the frequencies on each hit and says where +/// they came from. +fn frequencies( + dir: &str, + side: &dyn Fn(&str) -> Option<(String, String)>, + s: &str, + rows: &mut [Hit], +) -> FrequencySource { + let vep = side(&format!("vep/{s}_vep.vcf.gz")); + let _ = dir; + let mut note = None; + let mut found = HashMap::new(); + match &vep { + None => note = Some("no VEP output (step 13) in this folder".to_owned()), + Some((p, rel)) if !std::path::Path::new(&format!("{p}.tbi")).exists() => { + note = Some(format!("{rel} has no .tbi index")); + } + Some((p, rel)) => { + let keys: Vec<_> = rows.iter().map(Hit::key).collect(); + match vep_frequencies(p, &keys) { + Ok(f) => found = f, + Err(e) => note = Some(format!("{rel} could not be read: {e}")), + } + } + } + for x in rows.iter_mut() { + x.frequency = Some(found.get(&x.key()).cloned().unwrap_or(Frequency { + exome: None, + genome: None, + rsids: Vec::new(), + })); + } + FrequencySource { + source: if note.is_none() { + vep.map(|(_, rel)| rel) + } else { + None + }, + note, + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn floor_takes_the_last_value_of_each_key() { + let mut f = Floor::default(); + floor_values("##PGP_CALL_FLOOR=GQ=20,DP=8,GQ=25", false, &mut f); + assert_eq!(f, Floor { gq: 25, dp: 8 }); + let mut f = Floor::default(); + floor_values("##PGP_CALL_FLOOR=XGQ=5,DP=7", true, &mut f); + assert_eq!(f, Floor { gq: 30, dp: 7 }); + } + + #[test] + fn copies_and_words() { + assert_eq!(allele_copies("0/1"), (1, 2)); + assert_eq!(allele_copies("1|1"), (2, 2)); + assert_eq!(allele_copies("1"), (1, 1)); + assert_eq!(allele_copies("1/."), (1, 2)); + assert_eq!(allele_copies(""), (0, 0)); + assert_eq!(zygosity_word(1, 1), "hemizygous"); + assert_eq!(zygosity_word(0, 2), "unknown"); + } + + #[test] + fn rsid_fields() { + assert!(is_rsid("rs123")); + assert!(is_rsid("RS9")); + assert!(!is_rsid("rs")); + assert!(!is_rsid("COSV1")); + } +} diff --git a/rust/crates/og-summary/src/collect/cpsr.rs b/rust/crates/og-summary/src/collect/cpsr.rs new file mode 100644 index 00000000..ffd178d9 --- /dev/null +++ b/rust/crates/og-summary/src/collect/cpsr.rs @@ -0,0 +1,264 @@ +//! Cancer predisposition (step 17, CPSR): the class counts, and CPSR's own +//! class beside the final class and ClinVar's, so a final P/LP class that +//! ClinVar set over CPSR's own benign call is not hidden. + +use super::clinvar::{CdnaKey, Floor, cdna_checks, check_word, not_in_table, opt_int, opt_s}; +use super::text::{self, cell, digits, stripped}; +use indexmap::IndexMap; +use og_bam::cdna::{cpsr_variant, norm_class}; +use pycompat::csv::DictRow; +use pycompat::json::Value; +use std::collections::HashMap; + +const PLP: [&str; 2] = ["pathogenic", "likely pathogenic"]; +const CLASS_COLUMNS: [&str; 3] = [ + "CLASSIFICATION", + "FINAL_CLASSIFICATION", + "CPSR_CLASSIFICATION", +]; + +#[derive(Debug, Clone, PartialEq)] +pub struct Cpsr { + pub html_report: bool, + /// The class column and the counts, when the table was read. + pub table: Option, +} + +#[derive(Debug, Clone, PartialEq)] +pub struct Table { + /// `None`: the table has no class column. + pub column: Option, + pub classification: IndexMap, + pub compare: Option, +} + +#[derive(Debug, Clone, PartialEq)] +pub struct Compare { + pub plp: Vec, + pub disagreements: Vec, + pub own_benign: usize, + pub own: Option, + pub floor_dp: i128, + pub cdna: bool, +} + +/// Where the inputs beside the CPSR table are. +pub struct Inputs { + pub html: Option, + /// The classification table (full path). + pub tsv: Option, + /// Step 06's floor, the default without a hits VCF; an error is the + /// note of a hits VCF that cannot be read. + pub floor: Result, + /// Step 17's cDNA read check (full path). + pub cdna: Option, +} + +impl Cpsr { + pub fn read(i: &Inputs) -> Result { + let html_report = i.html.is_some(); + let Some(tsv) = &i.tsv else { + return Ok(Cpsr { + html_report, + table: None, + }); + }; + let rows = text::rows(tsv).map_err(|e| match e { + text::ReadError::Os(_) => e.unreadable(tsv), + text::ReadError::Other(n) => n, + })?; + // `rows[0].keys()`: the header, and the None key of a long first row. + let cols: Vec<&str> = rows + .first() + .map(|r| r.iter().map(|(k, _)| k).collect()) + .unwrap_or_default(); + let col = CLASS_COLUMNS.iter().find(|c| cols.contains(c)).copied(); + let Some(col) = col else { + return Ok(Cpsr { + html_report, + table: Some(Table { + column: None, + classification: IndexMap::new(), + compare: None, + }), + }); + }; + let mut counts: IndexMap = IndexMap::new(); + for r in &rows { + let k = cell(r, col).filter(|v| !v.is_empty()).unwrap_or("."); + *counts.entry(k.to_owned()).or_default() += 1; + } + counts.sort_by(|_, a, _, b| b.cmp(a)); + let compare = compare(i, &rows, &cols, col)?; + Ok(Cpsr { + html_report, + table: Some(Table { + column: Some(col.to_owned()), + classification: counts, + compare: Some(compare), + }), + }) + } + + pub fn to_json(&self) -> Value { + let mut m: IndexMap = IndexMap::new(); + m.insert("html_report".into(), Value::Bool(self.html_report)); + let Some(t) = &self.table else { + return Value::Object(m); + }; + m.insert("classification_column".into(), opt_s(t.column.as_deref())); + if t.column.is_some() { + m.insert( + "classification".into(), + Value::Object( + t.classification + .iter() + .map(|(k, v)| (k.clone(), Value::Int(*v as i128))) + .collect(), + ), + ); + } + if let Some(c) = &t.compare { + let field = |x: &Value, k: &str| match x { + Value::Object(m) => m.get(k).cloned(), + _ => None, + }; + let low = c + .plp + .iter() + .filter(|x| field(x, "low_confidence") == Some(Value::Bool(true))) + .count(); + m.insert("plp".into(), Value::Array(c.plp.clone())); + m.insert( + "disagreements".into(), + Value::Array(c.disagreements.clone()), + ); + m.insert( + "final_plp_own_benign".into(), + Value::Int(c.own_benign as i128), + ); + m.insert("own_class_column".into(), opt_s(c.own.as_deref())); + m.insert("low_confidence".into(), Value::Int(low as i128)); + m.insert("call_floor_dp".into(), Value::Int(c.floor_dp)); + if c.cdna { + let flagged = c + .plp + .iter() + .filter(|x| field(x, "cdna").as_ref().and_then(check_word) == Some("flagged")) + .count(); + m.insert("cdna_flagged".into(), Value::Int(flagged as i128)); + } + } + Value::Object(m) + } +} + +/// `cpsr_compare`. +fn compare(i: &Inputs, rows: &[DictRow], cols: &[&str], col: &str) -> Result { + let own = (cols.contains(&"CPSR_CLASSIFICATION") && col != "CPSR_CLASSIFICATION") + .then_some("CPSR_CLASSIFICATION"); + let floor = i.floor.clone()?; + let cdna: Option> = match &i.cdna { + Some(p) => Some(cdna_checks(p)?), + None => None, + }; + let has_clinvar_col = cols.contains(&"CLINVAR_CLASSIFICATION"); + let (mut plp, mut disagree, mut own_benign) = (Vec::new(), Vec::new(), 0); + for r in rows { + let final_ = norm_class(cell(r, col)); + let mine = own.map(|o| norm_class(cell(r, o))).unwrap_or_default(); + let authority = stripped(r, "ASSERTION_AUTHORITY"); + let clinvar_raw = if has_clinvar_col { + stripped(r, "CLINVAR_CLASSIFICATION") + } else if authority.to_lowercase() == "clinvar" { + stripped(r, col) + } else { + "" + }; + let clinvar = norm_class(Some(clinvar_raw)); + let is = |c: &str| PLP.contains(&c); + if is(&final_) && matches!(mine.as_str(), "benign" | "likely benign") { + own_benign += 1; + } + let is_plp = is(&final_) || is(&mine); + let differ = + own.is_some() && !clinvar.is_empty() && !mine.is_empty() && is(&clinvar) != is(&mine); + if !(is_plp || differ) { + continue; + } + let v = cpsr_variant(cell(r, "VAR_ID"), cell(r, "GENOMIC_CHANGE")) + .map_err(|e| format!("ValueError: {e}"))?; + let dp = digits(Some(stripped(r, "DP_CONTROL"))); + let or_dot = + |x: Option<&str>| Value::Str(x.filter(|v| !v.is_empty()).unwrap_or(".").to_owned()); + let mut x: IndexMap = IndexMap::new(); + x.insert( + "variant".into(), + match &v { + Some(h) => Value::Str(format!("{}:{} {}>{}", h.chrom, h.pos, h.reference, h.alt)), + None => or_dot(cell(r, "GENOMIC_CHANGE")), + }, + ); + x.insert("gene".into(), or_dot(cell(r, "SYMBOL"))); + x.insert("final".into(), or_dot(cell(r, col))); + x.insert( + "cpsr".into(), + own.map_or(Value::Null, |o| or_dot(cell(r, o))), + ); + x.insert( + "clinvar".into(), + if clinvar.is_empty() { + Value::Null + } else { + Value::Str(clinvar_raw.to_owned()) + }, + ); + x.insert( + "authority".into(), + opt_s(Some(authority).filter(|a| !a.is_empty())), + ); + x.insert( + "acmg_codes".into(), + opt_s(Some(stripped(r, "ACMG_CODE")).filter(|a| !a.is_empty())), + ); + x.insert("dp".into(), opt_int(dp)); + x.insert( + "low_confidence".into(), + Value::Bool(dp.is_some_and(|d| d < floor.dp)), + ); + if let Some(t) = &cdna + && is_plp + { + let found = v.as_ref().and_then(|h| { + t.get(&( + Some(h.chrom.clone()), + i128::from(h.pos), + Some(h.reference.clone()), + Some(h.alt.clone()), + )) + }); + x.insert( + "cdna".into(), + found + .filter(|c| pycompat::value::truthy(c)) + .cloned() + .unwrap_or_else(not_in_table), + ); + } + let x = Value::Object(x); + if differ { + disagree.push(x.clone()); + } + if is_plp { + plp.push(x); + } + } + Ok(Compare { + plp, + disagreements: disagree, + own_benign, + own: own.map(str::to_owned), + floor_dp: floor.dp, + cdna: cdna.is_some(), + }) +} diff --git a/rust/crates/og-summary/src/collect/group_b.rs b/rust/crates/og-summary/src/collect/group_b.rs new file mode 100644 index 00000000..cb0ba8a1 --- /dev/null +++ b/rust/crates/og-summary/src/collect/group_b.rs @@ -0,0 +1,284 @@ +//! Group b of `og summary collect`: the clinical sections. ClinVar, the +//! carrier list, PharmCAT, CPIC, pypgx, Cyrius, CPSR, the clinical filter +//! with its ACMG secondary-findings tier, slivar, and the `cyp2d6` block. +//! +//! Each input is read once: the ClinVar hits VCF feeds the ClinVar +//! section, the ACMG tier, the carrier list and CPSR's call floor; the +//! clinical VCF feeds the record count and the ACMG tier; AnnotSV's table +//! and step 35's copy-number tables feed both two-hit readings. The +//! ClinVar section (with its VEP lookups) and the clinical section (with +//! the clinical VCF) are read on their own threads. + +use super::acmg::{self, ClinicalVcf}; +use super::carrier::Carrier; +use super::clinical::{Clinical, Slivar}; +use super::clinvar::{ClinVar, Floor, HitsVcf}; +use super::cpsr::{self, Cpsr}; +use super::model::{Manifest, RunStatus, Section, Spec, Summary}; +use super::pgx::{Cpic, Cyp2d6, Cyrius, PharmCat, Pypgx}; +use super::scan::FolderIndex; +use super::sections::{Data, SampleQc}; +use super::two_hits::{self, DeletionInputs}; +use super::{Input, Read, read_at, read_optional, with_mtime}; +use std::path::Path; + +/// The sections `og summary collect --only b` reads, in summary.json's order. +pub const GROUP_B: [Spec; 9] = [ + Spec { + key: "clinvar", + title: "ClinVar screen", + step: "06", + }, + Spec { + key: "carrier", + title: "Carrier list", + step: "06", + }, + Spec { + key: "pharmcat", + title: "PharmCAT", + step: "07", + }, + Spec { + key: "cpic", + title: "CPIC drug-gene recommendations", + step: "27", + }, + Spec { + key: "pypgx", + title: "pypgx", + step: "32", + }, + Spec { + key: "cyrius", + title: "Cyrius (CYP2D6)", + step: "21", + }, + Spec { + key: "cpsr", + title: "Cancer predisposition (CPSR)", + step: "17", + }, + Spec { + key: "clinical", + title: "Clinical variant filter", + step: "23", + }, + Spec { + key: "slivar", + title: "Variant prioritization (slivar)", + step: "31", + }, +]; + +/// Reads the folder into the summary of group b. An error is a run status +/// or a manifest that is there and cannot be read. +pub fn collect(input: &Input) -> Result { + let dir = Path::new(&input.sample_dir); + let d = input.sample_dir.as_str(); + let index = FolderIndex::walk(dir); + let run = + read_optional(&dir.join("logs").join("run_status.tsv"))?.map(|t| RunStatus::parse(&t)); + let manifest = read_optional(&dir.join("run_manifest.tsv"))?.map(|t| Manifest::parse(&t)); + let s = pycompat::mtime::escape(&input.sample); + let raw = input.sample.as_str(); + let file = |rel: &str| -> (String, String) { (crate::join(d, rel), rel.to_owned()) }; + let find = |pats: &[String]| index.first(pats).map(file); + let one = |p: String| find(&[p]); + + // The shared inputs, each read once. + let hits = one(format!("clinvar/{s}_clinvar_hits.vcf")).map(|(p, rel)| HitsVcf::read(p, rel)); + let clinical_table = one(format!("clinical/{s}_clinical_summary.tsv")); + let clinical_vcf_at = one(format!("clinical/{s}_clinical.vcf.gz")); + let deletions = DeletionInputs { + annotsv: two_hits::read_annotsv( + one(format!("annotsv/{s}_sv_class4_5_checked.tsv")), + one(format!("sv_merged/{s}_sv_consensus.vcf.gz")).map(|(p, _)| p), + ), + step35: two_hits::read_step35(d, raw, &input.assets), + }; + let side = |rel: &str| one(rel.to_owned()); + + let (clinvar, clinical, rest) = std::thread::scope(|scope| { + let clinvar = scope.spawn(|| -> Read { + let Some(h) = &hits else { return Ok(None) }; + let data = ClinVar::read(h, d, &side, &s)?; + with_mtime(&h.path, h.rel.clone(), Data::ClinVar(data)) + }); + let clinical = scope.spawn(|| -> Read { + let at = clinical_table.clone().or_else(|| clinical_vcf_at.clone()); + let Some((p, rel)) = at else { return Ok(None) }; + let vcf = clinical_vcf_at.as_ref().map(|(v, _)| ClinicalVcf::read(v)); + let tier = acmg::Inputs { + assets: &input.assets, + hits: hits.as_ref().map(|h| &h.read), + clinical: vcf.as_ref(), + deletions: &deletions, + }; + let data = Clinical::read( + clinical_table.as_ref().map(|(p, _)| p.as_str()), + vcf.as_ref(), + &tier, + )?; + with_mtime(&p, rel, Data::Clinical(data)) + }); + let rest = read_rest(input, &s, &index, &find, &hits, &deletions); + let join = |h: std::thread::ScopedJoinHandle<'_, Read>| { + h.join() + .unwrap_or_else(|_| Err("internal error: the section reader stopped".into())) + }; + (join(clinvar), join(clinical), rest) + }); + let [carrier, pharmcat, cpic, pypgx, cyrius, cpsr, slivar] = rest; + let reads = [ + clinvar, carrier, pharmcat, cpic, pypgx, cyrius, cpsr, clinical, slivar, + ]; + let sections: Vec
= GROUP_B + .iter() + .zip(reads) + .map(|(spec, r)| Section::new(*spec, r, run.as_ref())) + .collect(); + let sec = |k: &str| sections.iter().find(|x| x.spec.key == k); + let cyp2d6 = match (sec("cpic"), sec("pypgx"), sec("cyrius")) { + (Some(c), Some(p), Some(y)) => Some(Cyp2d6::new( + c, + p, + y, + one(format!("pgx_consensus/{s}_pgx_consensus.tsv")), + d, + )), + _ => None, + }; + + // The run block is group a's: the declared sex from the command line, + // the run status, the manifest or step 33's table, the first that has one. + let qc_declared = one(format!("qc/{s}_sample_qc.tsv")).and_then(|(p, _)| { + let rows = super::sections::read_tsv(&p).ok()?; + SampleQc::read(&p, &rows) + .ok()? + .declared_sex() + .map(str::to_owned) + }); + let declared_sex = [ + input.declared_sex.clone(), + run.as_ref().and_then(|r| r.declared_sex.clone()), + manifest + .as_ref() + .and_then(|m| m.value("run", "declared_sex").map(str::to_owned)), + qc_declared, + ] + .into_iter() + .flatten() + .find(|v| !v.is_empty()); + let now = std::time::SystemTime::now() + .duration_since(std::time::UNIX_EPOCH) + .map_or(0.0, |d| d.as_secs_f64()); + Ok(Summary { + sample: input.sample.clone(), + generated_utc: super::time::utc_minute(now).unwrap_or_default(), + run, + manifest, + declared_sex, + sections, + cyp2d6, + }) +} + +/// Finds a section's file: the full path and the path in the folder. +type Find<'a> = dyn Fn(&[String]) -> Option<(String, String)> + 'a; + +/// The seven sections read on the calling thread, in order. +fn read_rest( + input: &Input, + s: &str, + index: &FolderIndex, + find: &Find, + hits: &Option, + deletions: &DeletionInputs, +) -> [Read; 7] { + let one = |p: String| find(&[p]); + let carrier = match hits { + None => Ok(None), + Some(h) => Carrier::read(&input.assets, &h.read, deletions) + .and_then(|c| with_mtime(&h.path, h.rel.clone(), Data::Carrier(c))), + }; + let pharmcat = pharmcat_at(index, input, s) + .and_then(|at| read_at(at, |p| PharmCat::read(p).map(Data::PharmCat))); + let cpic = read_at(one(format!("cpic/{s}_phenotypes.tsv")), |p| { + Cpic::read(p, &input.sample).map(Data::Cpic) + }); + let pypgx = read_at(one(format!("pypgx/{s}_pypgx_summary.tsv")), |p| { + Pypgx::read(p, &input.sample).map(Data::Pypgx) + }); + let cyrius = read_at( + find(&[ + format!("cyrius/{s}_cyp2d6.tsv"), + "cyrius/*_cyp2d6.tsv".to_owned(), + ]), + |p| Cyrius::read(p).map(Data::Cyrius), + ); + let html = find(&[ + format!("cpsr/{s}.cpsr.grch38.html"), + "cpsr/*.cpsr.*.html".to_owned(), + ]); + let tsv = find(&[ + format!("cpsr/{s}.cpsr.grch38.classification.tsv.gz"), + "cpsr/*.cpsr.*.classification.tsv.gz".to_owned(), + ]); + let cpsr = read_at(tsv.clone().or_else(|| html.clone()), |_| { + let floor = match hits { + None => Ok(Floor::default()), + Some(h) => h.read.as_ref().map(|(f, _)| *f).map_err(Clone::clone), + }; + Cpsr::read(&cpsr::Inputs { + html: html.as_ref().map(|(p, _)| p.clone()), + tsv: tsv.as_ref().map(|(p, _)| p.clone()), + floor, + cdna: one(format!("cpsr/{s}_cpsr_cdna_check.tsv")).map(|(p, _)| p), + }) + .map(Data::Cpsr) + }); + let slivar_table = one(format!("slivar/{s}_slivar_summary.tsv")); + let slivar_vcf = if slivar_table.is_none() { + one(format!("slivar/{s}_prioritized.vcf.gz")) + } else { + None + }; + let slivar = read_at(slivar_table.clone().or(slivar_vcf.clone()), |_| { + Slivar::read( + slivar_table.as_ref().map(|(p, _)| p.as_str()), + slivar_vcf.as_ref().map(|(p, _)| p.as_str()), + &input.sample, + ) + .map(Data::Slivar) + }); + [carrier, pharmcat, cpic, pypgx, cyrius, cpsr, slivar] +} + +/// PharmCAT's report: the newest of the files its four names match. +fn pharmcat_at( + index: &FolderIndex, + input: &Input, + s: &str, +) -> Result, String> { + let d = input.sample_dir.as_str(); + let found: Vec<&str> = [ + format!("pharmcat/{s}.report.json"), + "pharmcat/*.report.json".to_owned(), + format!("vcf/{s}.report.json"), + "vcf/*.report.json".to_owned(), + ] + .iter() + .flat_map(|p| index.all(p)) + .collect(); + let paths: Vec<(String, &str)> = found + .iter() + .map(|rel| (crate::join(d, rel), *rel)) + .collect(); + let newest = pycompat::mtime::newest(paths.iter().map(|(p, _)| p.as_str())) + .map_err(|e| e.to_string())?; + Ok(newest.map(|p| { + let rel = paths.iter().find(|(q, _)| q == p).map_or("", |(_, r)| r); + (p.to_owned(), rel.to_owned()) + })) +} diff --git a/rust/crates/og-summary/src/collect/model.rs b/rust/crates/og-summary/src/collect/model.rs index 689d8c08..c6d6e22a 100644 --- a/rust/crates/og-summary/src/collect/model.rs +++ b/rust/crates/og-summary/src/collect/model.rs @@ -295,6 +295,8 @@ pub struct Summary { pub manifest: Option, pub declared_sex: Option, pub sections: Vec
, + /// The CYP2D6 callers side by side (group b). + pub cyp2d6: Option, } impl Summary { @@ -334,6 +336,9 @@ impl Summary { .collect(), ), ); + if let Some(c) = &self.cyp2d6 { + m.insert("cyp2d6".into(), c.to_json()); + } Value::Object(m) } } diff --git a/rust/crates/og-summary/src/collect/pgx.rs b/rust/crates/og-summary/src/collect/pgx.rs new file mode 100644 index 00000000..d59db94c --- /dev/null +++ b/rust/crates/og-summary/src/collect/pgx.rs @@ -0,0 +1,664 @@ +//! The pharmacogenomics sections: PharmCAT (step 07), the CPIC table (step +//! 27), pypgx (step 32) and Cyrius (step 21), and the summary's `cyp2d6` +//! block, which sets the three CYP2D6 callers side by side by step 36's +//! rule (`og_steps::pgx_rules`). + +use super::clinvar::{obj, opt_int, opt_s}; +use super::model::{Section, State}; +use super::sections::Data; +use super::text::{self, cell, stripped}; +use indexmap::IndexMap; +use og_steps::pgx_rules::{NO_CALL, is_call, same_diplotype}; +use og_steps::pypgx_summary::{LOW_GQ, NOT_CALLABLE}; +use pycompat::csv::DictRow; +use pycompat::json::Value; +use pycompat::text::{lower, strip, upper}; + +/// The label of pypgx's CYP2D6 call when it is the only usable one. +pub const PROBABLE: &str = "probable (one caller: pypgx)"; +/// The rows of the CYP2D6 depth check the reports show. +pub const DEPTH_CHECK_KEYS: [&str; 5] = [ + "status", + "gene_mapq1_share", + "gene_mapq20_share", + "body_mapq1_share", + "body_mapq20_share", +]; + +fn rows_or_note(path: &str) -> Result, String> { + text::rows(path).map_err(|e| e.note()) +} + +/// `row.get(key, "")`: `""` without the column, `None` for a short row. +fn get_default(r: &DictRow, key: &str) -> Value { + match r.get(key) { + None => Value::Str(String::new()), + Some(v) => opt_s(v), + } +} + +/// `row[key]` as JSON: a `KeyError` note without the column. +fn need(r: &DictRow, key: &str) -> Result { + text::need(r, key).map(opt_s) +} + +// --- PharmCAT ------------------------------------------------------------------------ + +#[derive(Debug, Clone, PartialEq)] +pub struct PharmCat { + pub version: String, + pub html_report: bool, +} + +/// `json.load` of a file opened as text: universal newlines, and Python's +/// `JSONDecodeError` text (`MSG: line L column C (char N)`, in code points). +fn json_load(raw: &str) -> Result { + let text = raw.replace("\r\n", "\n").replace('\r', "\n"); + if text.starts_with('\u{feff}') { + return Err( + "Unexpected UTF-8 BOM (decode using utf-8-sig): line 1 column 1 (char 0)".into(), + ); + } + pycompat::json::loads(&text).map_err(|e| match e { + pycompat::json::JsonError::Decode(what, at) => { + let before = &text[..at.min(text.len())]; + let pos = before.chars().count(); + let line = before.matches('\n').count() + 1; + let col = match before.rfind('\n') { + Some(i) => before[i + 1..].chars().count() + 1, + None => pos + 1, + }; + format!("{what}: line {line} column {col} (char {pos})") + } + other => other.to_string(), + }) +} + +impl PharmCat { + pub fn read(path: &str) -> Result { + let raw = std::fs::read(path) + .map_err(|e| format!("{}: {}", crate::basename(path), crate::os_error(&e, path)))?; + let body = String::from_utf8(raw) + .map_err(|e| format!("{}: 'utf-8' codec can't decode: {e}", crate::basename(path)))?; + let data = json_load(&body).map_err(|e| format!("{}: {e}", crate::basename(path)))?; + let version = match &data { + Value::Object(m) => match m.get("pharmcatVersion") { + None => String::new(), + Some(v) => pycompat::value::py_str(v).map_err(|e| e.to_string())?, + }, + _ => String::new(), + }; + let html = path.strip_suffix(".json").map(|p| format!("{p}.html")); + Ok(PharmCat { + version, + html_report: html.is_some_and(crate::is_file), + }) + } + + pub fn to_json(&self) -> Value { + obj([ + ("version", Value::Str(self.version.clone())), + ("html_report", Value::Bool(self.html_report)), + ]) + } +} + +// --- CPIC ---------------------------------------------------------------------------- + +#[derive(Debug, Clone, PartialEq)] +pub struct CpicGene { + pub gene: Value, + pub diplotype: Value, + pub phenotype: Value, + pub status: Value, +} + +#[derive(Debug, Clone, PartialEq)] +pub struct Cpic { + pub genes: Vec, + pub parse_failed: bool, + pub warnings: Vec, +} + +impl Cpic { + /// `sec_cpic`: the phenotype table and, beside it, the recommendation + /// text's parse failure and warnings. + pub fn read(path: &str, s: &str) -> Result { + let genes: Vec = rows_or_note(path)? + .iter() + .map(|r| CpicGene { + gene: get_default(r, "Gene"), + diplotype: get_default(r, "Diplotype"), + phenotype: get_default(r, "Phenotype"), + status: get_default(r, "Status"), + }) + .collect(); + let rec = crate::join( + text::dirname(path), + &format!("{s}_cpic_recommendations.txt"), + ); + let (mut warnings, mut failed) = (Vec::new(), false); + if crate::is_file(&rec) { + for line in text::lines(&rec).map_err(|e| e.note())? { + if line.contains("PARSING FAILED") { + failed = true; + } + let t = strip(&line, None); + if let Some(w) = t.strip_prefix("WARNING:") { + warnings.push(strip(w, None).to_owned()); + } + } + } + let parse_failed = failed || genes.is_empty(); + Ok(Cpic { + genes, + parse_failed, + warnings, + }) + } + + fn count(&self, status: &str) -> Value { + Value::Int( + self.genes + .iter() + .filter(|g| matches!(&g.status, Value::Str(s) if s == status)) + .count() as i128, + ) + } + + /// The CYP2D6 row, if the table has one. + pub fn cyp2d6(&self) -> Option<&CpicGene> { + self.genes + .iter() + .find(|g| matches!(&g.gene, Value::Str(s) if s == "CYP2D6")) + } + + pub fn to_json(&self) -> Value { + obj([ + ( + "genes", + Value::Array( + self.genes + .iter() + .map(|g| { + obj([ + ("gene", g.gene.clone()), + ("diplotype", g.diplotype.clone()), + ("phenotype", g.phenotype.clone()), + ("status", g.status.clone()), + ]) + }) + .collect(), + ), + ), + ("non_normal", self.count("non-normal")), + ("ambiguous", self.count("ambiguous")), + ("not_called", self.count("not called")), + ("unclassified", self.count("unclassified")), + ("parse_failed", Value::Bool(self.parse_failed)), + ( + "warnings", + Value::Array(self.warnings.iter().cloned().map(Value::Str).collect()), + ), + ]) + } +} + +// --- pypgx --------------------------------------------------------------------------- + +#[derive(Debug, Clone, PartialEq)] +pub struct Pypgx { + pub genes_total: usize, + pub genes_called: usize, + pub cyp2d6: String, + /// Gene (as JSON's key text) to the copy-number call or why none. + pub cnv_call: IndexMap, + pub not_callable: Vec, + pub low_gq: Vec, + /// phenotype, defining_filter, quality of the CYP2D6 row. + pub cyp2d6_evidence: Option<[String; 3]>, + pub depth_check: Option>, + /// conflicts, one_tool_only, concordant. + pub comparison: Option<[usize; 3]>, +} + +/// A JSON object key from a cell that may be `None` (`json.dumps` writes `null`). +fn key_text(v: Option<&str>) -> String { + v.unwrap_or("null").to_owned() +} + +impl Pypgx { + pub fn read(path: &str, s: &str) -> Result { + let rows = rows_or_note(path)?; + let called = rows + .iter() + .filter(|r| !NO_CALL.contains(&lower(stripped(r, "Diplotype")).as_str())) + .count(); + let cyp_row = rows.iter().find(|r| cell(r, "Gene") == Some("CYP2D6")); + let cyp = cyp_row.and_then(|r| cell(r, "Diplotype")).unwrap_or(""); + let mut cnv_call = IndexMap::new(); + for r in &rows { + let status = stripped(r, "Status"); + if upper(stripped(r, "Source")) != "BAM" || status == NOT_CALLABLE { + continue; + } + let v = stripped(r, "CNV_call"); + let gene = key_text(text::need(r, "Gene")?); + let call = if !matches!(v, "" | ".") && !v.starts_with("N/A") { + v.to_owned() + } else if lower(status) == "failed" || upper(stripped(r, "Diplotype")) == "FAILED" { + "not called: pypgx failed".to_owned() + } else if status.contains("depth check") || v.contains("depth check") { + "not called: the CYP2D6 depth check failed".to_owned() + } else { + "not called: pypgx gave no copy-number call".to_owned() + }; + cnv_call.insert(gene, call); + } + let mut not_callable = Vec::new(); + for r in rows + .iter() + .filter(|r| stripped(r, "Status") == NOT_CALLABLE) + { + not_callable.push(need(r, "Gene")?); + } + let mut low_gq = Vec::new(); + for r in rows.iter().filter(|r| stripped(r, "Quality") == LOW_GQ) { + let gene = need(r, "Gene")?; + let min = stripped(r, "Min_GQ"); + let min_gq = if pycompat::text::isdigit(min) { + opt_int(crate::int(min).ok()) + } else { + Value::Null + }; + low_gq.push(obj([ + ("gene", gene), + ("diplotype", Value::Str(stripped(r, "Diplotype").to_owned())), + ("min_gq", min_gq), + ( + "weakest", + Value::Str(stripped(r, "Weakest_variant").to_owned()), + ), + ])); + } + let cyp2d6_evidence = cyp_row.map(|r| { + [ + stripped(r, "Phenotype").to_owned(), + stripped(r, "Defining_filter").to_owned(), + stripped(r, "Quality").to_owned(), + ] + }); + let dir = text::dirname(path); + let dc = crate::join(dir, &format!("{s}_cyp2d6_depth_check.tsv")); + let depth_check = if crate::is_file(&dc) { + let mut kv: IndexMap = IndexMap::new(); + for line in text::lines(&dc).map_err(|e| e.note())? { + if let Some((k, v)) = line.split_once('\t') { + kv.insert(k.to_owned(), v.to_owned()); + } + } + Some( + DEPTH_CHECK_KEYS + .iter() + .filter_map(|k| kv.get(*k).map(|v| ((*k).to_owned(), v.clone()))) + .collect(), + ) + } else { + None + }; + let comp = crate::join(dir, &format!("{s}_pharmcat_comparison.tsv")); + let comparison = if crate::is_file(&comp) { + let c = rows_or_note(&comp)?; + let n = + |f: &dyn Fn(Option<&str>) -> bool| c.iter().filter(|r| f(cell(r, "Match"))).count(); + Some([ + n(&|m| m == Some("No")), + n(&|m| matches!(m, Some("pypgx only" | "PharmCAT only"))), + n(&|m| m == Some("Yes")), + ]) + } else { + None + }; + Ok(Pypgx { + genes_total: rows.len(), + genes_called: called, + cyp2d6: if cyp.is_empty() { "not called" } else { cyp }.to_owned(), + cnv_call, + not_callable, + low_gq, + cyp2d6_evidence, + depth_check, + comparison, + }) + } + + pub fn to_json(&self) -> Value { + let mut m: IndexMap = IndexMap::new(); + m.insert("genes_total".into(), Value::Int(self.genes_total as i128)); + m.insert("genes_called".into(), Value::Int(self.genes_called as i128)); + m.insert("cyp2d6".into(), Value::Str(self.cyp2d6.clone())); + m.insert( + "cnv_call".into(), + Value::Object( + self.cnv_call + .iter() + .map(|(k, v)| (k.clone(), Value::Str(v.clone()))) + .collect(), + ), + ); + m.insert( + "not_callable".into(), + Value::Array(self.not_callable.clone()), + ); + m.insert("low_gq".into(), Value::Array(self.low_gq.clone())); + if let Some([p, f, q]) = &self.cyp2d6_evidence { + m.insert( + "cyp2d6_evidence".into(), + obj([ + ("phenotype", Value::Str(p.clone())), + ("defining_filter", Value::Str(f.clone())), + ("quality", Value::Str(q.clone())), + ]), + ); + } + if let Some(d) = &self.depth_check { + m.insert( + "depth_check".into(), + Value::Object( + d.iter() + .map(|(k, v)| (k.clone(), Value::Str(v.clone()))) + .collect(), + ), + ); + } + if let Some([c, o, y]) = self.comparison { + m.insert( + "comparison".into(), + obj([ + ("conflicts", Value::Int(c as i128)), + ("one_tool_only", Value::Int(o as i128)), + ("concordant", Value::Int(y as i128)), + ]), + ); + } + Value::Object(m) + } +} + +// --- Cyrius -------------------------------------------------------------------------- + +#[derive(Debug, Clone, PartialEq)] +pub struct Cyrius { + pub genotype: String, + /// `None` for a row cut short of the column. + pub filter: Option, +} + +impl Cyrius { + pub fn read(path: &str) -> Result { + let rows = rows_or_note(path)?; + let r = rows.first(); + let genotype = r.and_then(|r| cell(r, "Genotype")).unwrap_or(""); + let filter = match r.map(|r| r.get("Filter")) { + None | Some(None) => Some(String::new()), + Some(Some(v)) => v.map(str::to_owned), + }; + Ok(Cyrius { + genotype: if genotype.is_empty() { + "none" + } else { + genotype + } + .to_owned(), + filter, + }) + } + + pub fn to_json(&self) -> Value { + obj([ + ("genotype", Value::Str(self.genotype.clone())), + ("filter", opt_s(self.filter.as_deref())), + ]) + } +} + +// --- the cyp2d6 block ---------------------------------------------------------------- + +/// Step 36's CYP2D6 row, as the summary keeps it. +#[derive(Debug, Clone, PartialEq)] +pub struct Consensus { + pub result: String, + pub passed_to_pharmcat: bool, + pub reason: String, + pub source: Option, +} + +/// pypgx's call alone, labelled. +#[derive(Debug, Clone, PartialEq)] +pub struct Probable { + pub diplotype: String, + pub phenotype: String, +} + +/// The `cyp2d6` block of summary.json. +#[derive(Debug, Clone, PartialEq)] +pub struct Cyp2d6 { + pub pharmcat: Value, + pub pypgx: Value, + pub cyrius: Value, + pub agree: Option, + pub consensus: Option, + pub probable: Option, +} + +/// `probable_cyp2d6`: pypgx's call when it is the only usable one. +fn probable( + diplotype: &str, + evidence: Option<&[String; 3]>, + cyrius: Option<&Cyrius>, + depth_status: Option<&str>, +) -> Option { + if !is_call(diplotype) { + return None; + } + if let Some(c) = cyrius + && is_call(&c.genotype) + && c.filter.as_deref() == Some("PASS") + { + return None; + } + if depth_status.is_some_and(|d| d != "ok") { + return None; + } + let [phenotype, filter, quality] = evidence.cloned().unwrap_or_default(); + // A row without the evidence passes None, which neither list holds. + let have = evidence.is_some(); + if !have + || !matches!(filter.as_str(), "PASS" | "none") + || !matches!(quality.as_str(), "ok" | "N/A") + { + return None; + } + Some(Probable { + diplotype: strip(diplotype, None).to_owned(), + phenotype: strip(&phenotype, None).to_owned(), + }) +} + +/// `read_consensus`: step 36's CYP2D6 row, or `None` when the table is +/// missing, older than a caller's result it was made from, or has no +/// CYP2D6 row. +fn read_consensus( + found: Option<(String, String)>, + dir: &str, + sources: [Option<&str>; 2], +) -> Result, String> { + let Some((p, rel)) = found else { + return Ok(None); + }; + let mtime = |x: &str| pycompat::mtime::mtime(x).map_err(|e| crate::os_error(&e, x)); + let own = mtime(&p)?; + for src in sources.into_iter().flatten() { + if mtime(&crate::join(dir, src))? > own { + return Ok(None); + } + } + let rows = text::rows(&p).map_err(|e| e.text())?; + let Some(row) = rows.iter().find(|r| stripped(r, "Gene") == "CYP2D6") else { + return Ok(None); + }; + Ok(Some(Consensus { + result: text::cell_or(row, "Result").to_owned(), + passed_to_pharmcat: cell(row, "Outside_call") == Some("yes"), + reason: text::cell_or(row, "Reason").to_owned(), + source: Some(rel), + })) +} + +impl Cyp2d6 { + /// `cyp2d6_block` over the three caller sections. + pub fn new( + cpic: &Section, + pypgx: &Section, + cyrius: &Section, + consensus_at: Option<(String, String)>, + dir: &str, + ) -> Cyp2d6 { + let pc = match &cpic.data { + Some(Data::Cpic(c)) => c.cyp2d6(), + _ => None, + }; + let pd = match (&pypgx.data, pypgx.state) { + (Some(Data::Pypgx(p)), State::Ok) => Some(p), + _ => None, + }; + let cy = match (&cyrius.data, cyrius.state) { + (Some(Data::Cyrius(c)), State::Ok) => Some(c), + _ => None, + }; + let failed = |s: &Section| { + (s.state == State::Failed).then(|| Value::Str(format!("failed (step {})", s.spec.step))) + }; + let pg = pd.map(|p| p.cyp2d6.as_str()); + let pg_ok = pg.is_some_and(is_call); + let cy_ok = cy.is_some_and(|c| is_call(&c.genotype) && c.filter.as_deref() == Some("PASS")); + let cy_txt = cy.map(|c| { + let mut t = c.genotype.clone(); + if is_call(&t) && !cy_ok { + let f = c + .filter + .as_deref() + .filter(|f| !f.is_empty()) + .unwrap_or("missing"); + t.push_str(&format!(" (not usable: Filter {f})")); + } + t + }); + let pharmcat = match pc { + Some(g) => { + let st = &g.status; + if matches!(st, Value::Str(s) if s == "not called" || s == "ambiguous") { + st.clone() + } else { + g.diplotype.clone() + } + } + None => failed(cpic).unwrap_or(Value::Null), + }; + let agree = match (pg_ok && cy_ok, pg, cy) { + (true, Some(a), Some(c)) => Some(same_diplotype(a, &c.genotype)), + _ => None, + }; + let consensus = match read_consensus( + consensus_at, + dir, + [pypgx.source.as_deref(), cyrius.source.as_deref()], + ) { + Ok(c) => c, + Err(e) => Some(Consensus { + result: "unreadable".into(), + passed_to_pharmcat: false, + reason: e, + source: None, + }), + }; + let probable = match (pg_ok, pd) { + (true, Some(p)) => probable( + &p.cyp2d6, + p.cyp2d6_evidence.as_ref(), + cy, + p.depth_check + .as_ref() + .and_then(|d| d.get("status")) + .map(String::as_str), + ), + _ => None, + }; + Cyp2d6 { + pharmcat, + pypgx: failed(pypgx).unwrap_or_else(|| opt_s(pg)), + cyrius: failed(cyrius).unwrap_or_else(|| opt_s(cy_txt.as_deref())), + agree, + consensus, + probable, + } + } + + pub fn to_json(&self) -> Value { + obj([ + ( + "calls", + obj([ + ("PharmCAT", self.pharmcat.clone()), + ("pypgx", self.pypgx.clone()), + ("Cyrius", self.cyrius.clone()), + ]), + ), + ("agree", self.agree.map_or(Value::Null, Value::Bool)), + ( + "consensus", + self.consensus.as_ref().map_or(Value::Null, |c| { + obj([ + ("result", Value::Str(c.result.clone())), + ("passed_to_pharmcat", Value::Bool(c.passed_to_pharmcat)), + ("reason", Value::Str(c.reason.clone())), + ("source", opt_s(c.source.as_deref())), + ]) + }), + ), + ( + "probable", + self.probable.as_ref().map_or(Value::Null, |p| { + obj([ + ("diplotype", Value::Str(p.diplotype.clone())), + ("phenotype", Value::Str(p.phenotype.clone())), + ("label", Value::Str(PROBABLE.into())), + ]) + }), + ), + ]) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn json_errors_read_as_python_words_them() { + assert_eq!( + json_load("{\"pharmcatVersion\": ").unwrap_err(), + "Expecting value: line 1 column 21 (char 20)" + ); + assert_eq!( + json_load("{\r\n\"é\": x").unwrap_err(), + "Expecting value: line 2 column 6 (char 7)" + ); + assert!( + json_load("\u{feff}{}") + .unwrap_err() + .starts_with("Unexpected UTF-8 BOM") + ); + assert!(json_load("[1]").is_ok()); + } +} diff --git a/rust/crates/og-summary/src/collect/scan.rs b/rust/crates/og-summary/src/collect/scan.rs index 26d7c008..2792c67d 100644 --- a/rust/crates/og-summary/src/collect/scan.rs +++ b/rust/crates/og-summary/src/collect/scan.rs @@ -87,6 +87,25 @@ impl FolderIndex { } None } + + /// Every file one pattern names, in sorted order (`glob` with the + /// order made stable), for the readers that take the newest of them. + pub fn all(&self, pattern: &str) -> Vec<&str> { + if !has_magic(pattern) { + return self + .files + .get(pattern) + .map(String::as_str) + .into_iter() + .collect(); + } + let pat: Vec<&str> = pattern.split('/').collect(); + self.sorted + .iter() + .filter(|f| matches(&f.split('/').collect::>(), &pat)) + .map(String::as_str) + .collect() + } } fn has_magic(s: &str) -> bool { diff --git a/rust/crates/og-summary/src/collect/sections.rs b/rust/crates/og-summary/src/collect/sections.rs index 8aa04d53..a719a593 100644 --- a/rust/crates/og-summary/src/collect/sections.rs +++ b/rust/crates/og-summary/src/collect/sections.rs @@ -20,6 +20,15 @@ pub enum Data { MosaicSexLoss(Mosaic), ReadCount(ReadCount), Variants(Variants), + ClinVar(super::clinvar::ClinVar), + Carrier(super::carrier::Carrier), + PharmCat(super::pgx::PharmCat), + Cpic(super::pgx::Cpic), + Pypgx(super::pgx::Pypgx), + Cyrius(super::pgx::Cyrius), + Cpsr(super::cpsr::Cpsr), + Clinical(super::clinical::Clinical), + Slivar(super::clinical::Slivar), } impl Data { @@ -32,6 +41,15 @@ impl Data { Data::MosaicSexLoss(d) => d.to_json(), Data::ReadCount(d) => d.to_json(), Data::Variants(d) => d.to_json(), + Data::ClinVar(d) => d.to_json(), + Data::Carrier(d) => d.to_json(), + Data::PharmCat(d) => d.to_json(), + Data::Cpic(d) => d.to_json(), + Data::Pypgx(d) => d.to_json(), + Data::Cyrius(d) => d.to_json(), + Data::Cpsr(d) => d.to_json(), + Data::Clinical(d) => d.to_json(), + Data::Slivar(d) => d.to_json(), } } } diff --git a/rust/crates/og-summary/src/collect/text.rs b/rust/crates/og-summary/src/collect/text.rs new file mode 100644 index 00000000..275e236d --- /dev/null +++ b/rust/crates/og-summary/src/collect/text.rs @@ -0,0 +1,270 @@ +//! The small readers the clinical sections (group b) share: a VCF's data +//! lines as cells, INFO and FORMAT as maps, a table whose read errors keep +//! their kind, and Python's way of reading a digit cell. + +use crate::{basename, os_error, text_lines}; +use pycompat::csv::{DictReader, DictRow, Quoting, Reader}; +use pycompat::text::{isdigit, strip}; +use std::collections::HashMap; + +/// Why a file could not be read: an `OSError` (its `str()`), or a fault of +/// another kind, already written as the section's note. +#[derive(Debug, Clone)] +pub enum ReadError { + Os(String), + Other(String), +} + +impl ReadError { + /// The note when the Python wrapped the error in `Unreadable(f"{name}: {e}")`. + pub fn unreadable(&self, path: &str) -> String { + match self { + ReadError::Os(e) => format!("{}: {e}", basename(path)), + ReadError::Other(n) => n.clone(), + } + } + + /// The note of the collector's catch-all: `type(e).__name__: e`. + pub fn note(&self) -> String { + match self { + ReadError::Os(e) => format!("{}: {e}", os_kind(e)), + ReadError::Other(n) => n.clone(), + } + } + + /// `str(e)`, as a message that embeds the error prints it. + pub fn text(&self) -> String { + match self { + ReadError::Os(e) => e.clone(), + ReadError::Other(n) => n.split_once(": ").map_or(n.clone(), |(_, m)| m.to_owned()), + } + } +} + +fn os_kind(e: &str) -> &'static str { + if e == GZ_EOF { + "EOFError" + } else if e.contains("No such file") { + "FileNotFoundError" + } else if e.contains("Permission denied") { + "PermissionError" + } else if e.contains("Is a directory") { + "IsADirectoryError" + } else { + "OSError" + } +} + +/// The lines of a text file (gzip by its magic, universal newlines, bad +/// UTF-8 as U+FFFD), each without its line ending. +pub fn lines(path: &str) -> Result, ReadError> { + let it = og_io::Lines::new( + og_io::open(path).map_err(|e| ReadError::Os(os_error(&e, path)))?, + og_io::Utf8::Replace, + ); + it.map(|l| l.map_err(|e| ReadError::Os(gz_text(&e)))) + .collect() +} + +/// Python's gzip on a stream cut short (an `EOFError`). +const GZ_EOF: &str = "Compressed file ended before the end-of-stream marker was reached"; + +/// gzip's own words for a stream cut short; the reader's otherwise. +fn gz_text(e: &std::io::Error) -> String { + if e.kind() == std::io::ErrorKind::UnexpectedEof { + GZ_EOF.to_owned() + } else { + e.to_string() + } +} + +/// Calls `f` with each line of a text file in turn (as [`lines`] reads +/// them), never holding the whole file. +pub fn each_line(path: &str, mut f: impl FnMut(&str)) -> Result<(), ReadError> { + let it = og_io::Lines::new( + og_io::open(path).map_err(|e| ReadError::Os(os_error(&e, path)))?, + og_io::Utf8::Replace, + ); + for l in it { + f(&l.map_err(|e| ReadError::Os(gz_text(&e)))?); + } + Ok(()) +} + +/// Table lines keep their endings so quoted fields may span lines. +fn table_lines(path: &str) -> Result>, ReadError> { + let mut lines = text_lines(path).map_err(ReadError::Os)?; + Ok(std::iter::from_fn(move || { + let mut s = String::new(); + match lines.read_line(&mut s) { + Ok(true) => Some(Ok(s)), + Ok(false) => None, + Err(e) if e.kind() == std::io::ErrorKind::UnexpectedEof => { + Some(Err(std::io::Error::new(e.kind(), gz_text(&e)))) + } + Err(e) => Some(Err(e)), + } + })) +} + +fn csv_error(e: pycompat::csv::Error) -> ReadError { + match e.exception() { + "OSError" => ReadError::Os(e.to_string()), + kind => ReadError::Other(format!("{kind}: {e}")), + } +} + +/// The rows of a tab-separated table with a header (`csv.DictReader`). +pub fn rows(path: &str) -> Result, ReadError> { + DictReader::new(table_lines(path)?, '\t', Quoting::Minimal) + .collect::, _>>() + .map_err(csv_error) +} + +/// The header of a table as `csv.reader`'s first row reads it. +pub fn columns(path: &str) -> Result, ReadError> { + Reader::new(table_lines(path)?, '\t', Quoting::Minimal) + .next() + .transpose() + .map(Option::unwrap_or_default) + .map_err(csv_error) +} + +/// `row.get(key)` where a missing key and a short row both read as `None`. +pub fn cell<'a>(row: &'a DictRow, key: &str) -> Option<&'a str> { + row.get(key).flatten() +} + +/// `row.get(key) or ""`. +pub fn cell_or<'a>(row: &'a DictRow, key: &str) -> &'a str { + cell(row, key).unwrap_or("") +} + +/// `(row.get(key) or "").strip()`. +pub fn stripped<'a>(row: &'a DictRow, key: &str) -> &'a str { + strip(cell_or(row, key), None) +} + +/// `row[key]`: a `KeyError` note when the table has no such column. A +/// short row gives `None`. +pub fn need<'a>(row: &'a DictRow, key: &str) -> Result, String> { + row.get(key) + .ok_or_else(|| format!("KeyError: {}", crate::repr(key))) +} + +/// `int(x) if x.isdigit() else None`. +pub fn digits(x: Option<&str>) -> Option { + x.filter(|v| isdigit(v)).and_then(|v| crate::int(v).ok()) +} + +/// `info_map`: a VCF INFO field as a map, a later key over an earlier one. +pub fn info_map(field: &str) -> HashMap<&str, &str> { + field + .split(';') + .map(|kv| kv.split_once('=').unwrap_or((kv, ""))) + .collect() +} + +/// `dict(zip(FORMAT.split(":"), SAMPLE.split(":")))`. +pub fn format_map<'a>(fmt: &'a str, sample: &'a str) -> HashMap<&'a str, &'a str> { + fmt.split(':').zip(sample.split(':')).collect() +} + +/// One FORMAT field of the first sample, `None` when absent, empty or `.`. +pub fn gt_field(r: &[String], key: &str) -> Option { + if r.len() < 10 { + return None; + } + format_map(&r[8], &r[9]) + .get(key) + .filter(|v| !v.is_empty() && **v != ".") + .map(|v| (*v).to_owned()) +} + +/// The GT of the first sample, `""` without one. +pub fn gt_of(r: &[String]) -> &str { + r.get(9).map_or("", |s| s.split(':').next().unwrap_or("")) +} + +/// `int(r[1])`, as the note of the collector's catch-all when it is not one. +pub fn pos(r: &[String]) -> Result { + crate::int(&r[1]).map_err(|e| format!("ValueError: {e}")) +} + +/// Python's `max()` over floats: the first of the largest, NaN kept only +/// when it comes first. +pub fn py_max(xs: impl IntoIterator) -> Option { + xs.into_iter().fold(None, |best, x| match best { + Some(b) if x.partial_cmp(&b) != Some(std::cmp::Ordering::Greater) => Some(b), + _ => Some(x), + }) +} + +/// `os.path.relpath(path, base)` for a path found under `base`. +pub fn relpath(base: &str, path: &str) -> String { + let b = base.trim_end_matches('/'); + path.strip_prefix(b) + .and_then(|r| r.strip_prefix('/')) + .unwrap_or(path) + .to_owned() +} + +/// `os.path.dirname(p)`. +pub fn dirname(p: &str) -> &str { + p.rfind('/').map_or("", |i| &p[..i]) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::Dir; + + #[test] + fn quoted_columns_match_the_row_keys() { + let d = Dir::new("text-columns"); + for header in [ + "\"chrom\"\tstart\tend", + "\"a\tb\"\t\"a\"\"b\"\t\"a\nb\"", + "\"\"\tstart\tend", + ] { + let p = d.put("quoted.tsv", format!("{header}\n1\t2\t3\n").as_bytes()); + let r = rows(&p).unwrap(); + let keys: Vec<_> = r[0].iter().map(|(k, _)| k.to_owned()).collect(); + assert_eq!(columns(&p).unwrap(), keys, "{header}"); + } + } + + #[test] + fn columns_read_only_the_first_record() { + let d = Dir::new("text-first-record"); + // A row is larger than the CSV field limit. Reading the header + // must not parse or retain that unrelated row. + let p = d.put( + "large.tsv", + format!("\"chrom\"\tstart\n{}\t1\n", "a".repeat(140_000)).as_bytes(), + ); + assert_eq!(columns(&p).unwrap(), ["chrom", "start"]); + for (header, expected) in [ + ("", vec![]), + ("\n", vec![]), + ("\"\"\n", vec![""]), + ("\"chrom\"\tstart", vec!["chrom", "start"]), + ] { + let p = d.put("empty.tsv", header.as_bytes()); + assert_eq!(columns(&p).unwrap(), expected); + } + } + + #[test] + fn a_table_cut_short_is_python_s_eof_error() { + let d = Dir::new("text-cut"); + let body = "A\tB\n".to_owned() + &"1\t2\n".repeat(2000); + let full = std::fs::read(d.put_gz("t.tsv.gz", &body)).unwrap(); + for n in [30, full.len() / 2, full.len() - 4] { + let p = d.put("cut.tsv.gz", &full[..n]); + let e = rows(&p).unwrap_err(); + assert_eq!(e.unreadable(&p), format!("cut.tsv.gz: {GZ_EOF}"), "{n}"); + assert_eq!(e.note(), format!("EOFError: {GZ_EOF}"), "{n}"); + } + } +} diff --git a/rust/crates/og-summary/src/collect/two_hits.rs b/rust/crates/og-summary/src/collect/two_hits.rs new file mode 100644 index 00000000..b022c841 --- /dev/null +++ b/rust/crates/og-summary/src/collect/two_hits.rs @@ -0,0 +1,1001 @@ +//! Two hits in one gene across callers: a called deletion plus a small +//! variant. +//! +//! Small variants alone miss a gene hit twice when one hit is a deletion: +//! an exon deletion of one copy of MUTYH plus a pathogenic variant on the +//! other, an HBA deletion plus an HBA point variant, one SMN1 copy plus an +//! SMN1 variant, a 1q21.1 deletion plus the RBM8A low-expression allele +//! (TAR syndrome: Albers et al., Nat Genet 2012, PMID 22366785). The +//! deletions come from the engines that call them: AnnotSV's class 4 and 5 +//! deletions (step 05), step 35's interval table, its HBA count and +//! Parascopy's SMN1 count. +//! +//! The two are in trans (biallelic) only when the reads show it: the +//! variant sits inside the deletion and is called homozygous (one copy is +//! left, and it carries the variant), or both calls carry read-backed phase +//! (a `|` GT and the same PS) on different copies. On the same copy (cis) +//! they are no second hit; otherwise they are possibly in trans. +//! +//! The readers here take from step 05's and step 35's tables only what the +//! rule needs. The AnnotSV and Parascopy sections themselves are read by +//! their own modules. + +use super::clinvar::obj; +use super::text::{self, ReadError, cell, columns, digits, rows}; +use indexmap::IndexMap; +use pycompat::csv::DictRow; +use pycompat::json::Value; +use pycompat::text::isdigit; +use std::collections::{BTreeSet, HashMap}; + +/// Gene spans: chrom, start, end (1-based, inclusive). +pub type Spans = IndexMap; + +/// A called deletion that reaches a gene. +#[derive(Debug, Clone, PartialEq)] +pub struct Deletion { + pub chrom: String, + pub start: i128, + pub end: i128, + pub text: String, + /// A variant inside the span is on the one copy left. False for the + /// HBA count, which cannot say which gene of the pair is gone. + pub inside: bool, + /// The call's GT and phase set (AnnotSV only). + pub gt: Option, + pub ps: Option, +} + +/// A small variant of a gene, as the rule reads it. +#[derive(Debug, Clone, PartialEq)] +pub struct Variant { + pub gene: String, + pub key: VKey, + pub chrom: String, + pub pos: i128, + pub text: String, + pub genotype: String, + pub copies: usize, + pub ploidy: usize, + pub gt: String, + pub ps: Option, + /// The severity name of a hypomorph allele. + pub allele: Option, +} + +/// A variant key: a small variant `(chrom, pos, ref, alt)`. +pub type VKey = (String, i128, String, String); + +/// How a deletion and a small variant in one gene pair up. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)] +pub enum Verdict { + /// Both copies are hit. + Biallelic, + /// Phase unknown, or the calls disagree: possibly both copies. + PossiblyBiallelic, + /// Both on the same copy: no second hit. + Cis, +} + +impl Verdict { + pub fn word(self) -> &'static str { + match self { + Verdict::Biallelic => "biallelic", + Verdict::PossiblyBiallelic => "possibly biallelic", + Verdict::Cis => "cis", + } + } +} + +/// The best pair of a gene. +#[derive(Debug, Clone, PartialEq)] +pub struct Entry { + pub gene: String, + pub verdict: Verdict, + pub why: String, + pub deletion: String, + pub variant: String, + pub genotype: String, + pub allele: Option, +} + +impl Entry { + pub fn to_json(&self) -> IndexMap { + let Value::Object(m) = obj([ + ("gene", Value::Str(self.gene.clone())), + ("verdict", Value::Str(self.verdict.word().into())), + ("why", Value::Str(self.why.clone())), + ("deletion", Value::Str(self.deletion.clone())), + ("variant", Value::Str(self.variant.clone())), + ("genotype", Value::Str(self.genotype.clone())), + ( + "allele", + self.allele.clone().map_or(Value::Null, Value::Str), + ), + ]) else { + unreachable!() + }; + m + } +} + +/// The copies a phased GT puts the alt allele on, or `None` when the GT +/// is not phased or a copy is not called (the missing copy could carry the +/// other allele, so no cis or trans can be read from it). +pub fn phase_haps(gt: Option<&str>) -> Option> { + let gt = gt.filter(|g| !g.is_empty() && g.contains('|'))?; + let alleles: Vec<&str> = gt.split('|').collect(); + if alleles.contains(&".") { + return None; + } + Some( + alleles + .iter() + .enumerate() + .filter(|(_, a)| **a != "0") + .map(|(i, _)| i) + .collect(), + ) +} + +/// The verdict of one deletion and one small variant in one gene, and why. +pub fn pair_phase(dl: &Deletion, v: &Variant) -> (Verdict, String) { + let inside = dl.inside && dl.chrom == v.chrom && dl.start <= v.pos && v.pos <= dl.end; + // Read-backed phase first: on the same copy (cis) is no second hit, + // wherever the variant sits. + let hd = phase_haps(dl.gt.as_deref()).filter(|h| !h.is_empty()); + let hv = phase_haps(Some(&v.gt)).filter(|h| !h.is_empty()); + let mut phased: Option<(Verdict, String)> = None; + if let (Some(hd), Some(hv)) = (&hd, &hv) + && let Some(ps) = dl.ps.as_deref().filter(|p| !p.is_empty()) + && Some(ps) == v.ps.as_deref() + { + if hd == hv { + return ( + Verdict::Cis, + format!("in cis: read-backed phase puts both on the same copy (phase set {ps})"), + ); + } + phased = Some(if hd.is_disjoint(hv) { + ( + Verdict::Biallelic, + format!( + "in trans: read-backed phase puts them on different copies (phase set {ps})" + ), + ) + } else { + ( + Verdict::PossiblyBiallelic, + format!( + "read-backed phase puts the variant on both copies (phase set {ps}); check the reads" + ), + ) + }); + } + if inside && v.copies >= 2 { + return ( + Verdict::Biallelic, + "the variant is called homozygous inside the deletion: it is hemizygous, on the one \ + copy the deletion leaves" + .into(), + ); + } + if inside && v.ploidy == 2 { + return ( + Verdict::PossiblyBiallelic, + "a heterozygous call inside a called deletion: one of the two calls may be wrong; \ + check the reads" + .into(), + ); + } + if let Some(p) = phased { + return p; + } + ( + Verdict::PossiblyBiallelic, + "phase unknown, possibly in trans (no read-backed phase links the two calls)".into(), + ) +} + +/// What [`two_hits`] finds. +#[derive(Debug, Default)] +pub struct TwoHits { + /// The best pair of each gene hit twice, in the order genes were found. + pub entries: IndexMap, + /// Variants called homozygous inside a deletion (hemizygous): the + /// deletion's text. + pub hemi: HashMap<(String, VKey), String>, + /// The best pair of each variant that has one: ` with + /// : `. + pub pairs: HashMap<(String, VKey), String>, +} + +/// The genes hit twice. +pub fn two_hits(variants: &[Variant], dels: &HashMap>) -> TwoHits { + let mut out = TwoHits::default(); + let mut best: HashMap<(String, VKey), Verdict> = HashMap::new(); + for v in variants { + let k = (v.gene.clone(), v.key.clone()); + for dl in dels.get(&v.gene).map_or(&[][..], Vec::as_slice) { + let (verdict, why) = pair_phase(dl, v); + if verdict == Verdict::Biallelic && why.starts_with("the variant is called homozygous") + { + out.hemi.insert(k.clone(), dl.text.clone()); + } + if verdict == Verdict::Cis { + continue; + } + let e = Entry { + gene: v.gene.clone(), + verdict, + why: why.clone(), + deletion: dl.text.clone(), + variant: v.text.clone(), + genotype: if out.hemi.contains_key(&k) { + "hemizygous".into() + } else { + v.genotype.clone() + }, + allele: v.allele.clone(), + }; + match out.entries.get(&v.gene) { + Some(old) if old.verdict <= verdict => {} + _ => { + out.entries.insert(v.gene.clone(), e); + } + } + if best.get(&k).is_none_or(|b| verdict < *b) { + best.insert(k.clone(), verdict); + out.pairs.insert( + k.clone(), + format!("{} with {}: {why}", verdict.word(), dl.text), + ); + } + } + } + out +} + +// --- the deletions -------------------------------------------------------------------- + +/// One AnnotSV row as the two-hit and ACMG rules read it. +#[derive(Debug, Clone, PartialEq)] +pub struct SvRow { + /// `SV_chrom:SV_start-SV_end` as AnnotSV writes it. + pub sv: String, + pub sv_type: Option, + pub gt: Option, + pub ps: Option, + pub acmg_class: Option, + pub genes: Vec, + pub caller: Option, + /// Callers supporting it in step 22's consensus; `None` without one. + pub support: Option, + pub common: bool, +} + +impl SvRow { + /// `str(acmg_class)`. + pub fn class(&self) -> &str { + self.acmg_class.as_deref().unwrap_or("None") + } +} + +/// Step 05's checked class 4 and 5 table: its class 4 and 5 rows, `None` +/// when there is no table, or the error text the rules print. +pub fn read_annotsv( + at: Option<(String, String)>, + consensus: Option, +) -> Result>, String> { + let Some((p, _)) = at else { + return Ok(None); + }; + let size = std::fs::metadata(&p) + .map_err(|e| crate::os_error(&e, &p))? + .len(); + let rows = if size > 0 { + rows(&p).map_err(|e| e.text())? + } else { + Vec::new() + }; + const NEED: [&str; 9] = [ + "ACMG_class", + "Caller", + "Common", + "Exon_edge", + "Population_AF", + "SV_chrom", + "SV_end", + "SV_start", + "SV_type", + ]; + if let Some(r) = rows.first() + && NEED.iter().any(|k| r.get(k).is_none()) + { + let list: Vec = NEED.iter().map(|k| crate::repr(k)).collect(); + return Err(format!( + "{} lacks a column of [{}]", + crate::basename(&p), + list.join(", ") + )); + } + let support = consensus_support(&rows, consensus.as_deref()); + let py = |v: Option<&str>| v.unwrap_or("None").to_owned(); + let out = rows + .iter() + .zip(support) + .map(|(r, support)| { + let fmt = sample_format(r); + let some = |k: &str| { + fmt.get(k) + .filter(|v| !v.is_empty()) + .map(|v| (*v).to_owned()) + }; + SvRow { + sv: format!( + "{}:{}-{}", + py(cell(r, "SV_chrom")), + py(cell(r, "SV_start")), + py(cell(r, "SV_end")) + ), + sv_type: cell(r, "SV_type").map(str::to_owned), + gt: some("GT"), + ps: some("PS"), + acmg_class: cell(r, "ACMG_class").map(str::to_owned), + genes: text::cell_or(r, "Gene_name") + .split(';') + .filter(|g| !g.is_empty()) + .map(str::to_owned) + .collect(), + caller: cell(r, "Caller").map(str::to_owned), + support, + common: cell(r, "Common") == Some("yes"), + } + }) + .filter(|x| matches!(x.class(), "4" | "5")) + .collect(); + Ok(Some(out)) +} + +/// The first sample's FORMAT values of an AnnotSV row: the column after +/// FORMAT in the header, none when that is `Annotation_mode` or absent. +fn sample_format(r: &DictRow) -> HashMap<&str, &str> { + let keys: Vec<&str> = r.iter().map(|(k, _)| k).collect(); + let Some(i) = keys.iter().position(|k| *k == "FORMAT") else { + return HashMap::new(); + }; + match keys.get(i + 1) { + Some(&k) if k != "Annotation_mode" => { + text::format_map(text::cell_or(r, "FORMAT"), text::cell_or(r, k)) + } + _ => HashMap::new(), + } +} + +/// How far a consensus record may sit from an AnnotSV call, each end. +const CONSENSUS_BP: i128 = 1000; + +fn chrom_key(c: &str) -> &str { + if c.as_bytes() + .get(..3) + .is_some_and(|prefix| prefix.eq_ignore_ascii_case(b"chr")) + { + // An ASCII prefix ends at a UTF-8 boundary. + &c[3..] + } else { + c + } +} + +fn info_value<'a>(info: &'a str, key: &str) -> &'a str { + info.split(';') + .find_map(|kv| { + let (k, v) = kv.split_once('=').unwrap_or((kv, "")); + (k == key).then_some(v) + }) + .unwrap_or("") +} + +fn dup_or(t: &str) -> String { + let t = t.to_uppercase(); + if t.starts_with("DUP") { + "DUP".into() + } else { + t + } +} + +/// The support of each row from step 22's consensus: the nearest record of +/// the same contig and type within 1 kb at both ends; 1 without one; `None` +/// for every row when there is no consensus or it cannot be read. +fn consensus_support(rows: &[DictRow], path: Option<&str>) -> Vec> { + let none = || vec![None; rows.len()]; + let Some(path) = path.filter(|p| crate::is_file(p)) else { + return none(); + }; + let Ok(lines) = text::lines(path) else { + return none(); + }; + // (contig, type) to the records' (pos, end, supp). + type Records = HashMap<(String, String), Vec<(i128, i128, i128)>>; + let mut by: Records = HashMap::new(); + for l in lines.iter().filter(|l| !l.starts_with('#')) { + let c: Vec<&str> = l.split('\t').collect(); + if c.len() < 8 { + continue; + } + let num = |s: &str| crate::int(s).ok(); + let Some(pos) = num(c[1]) else { continue }; + let end = match info_value(c[7], "END") { + "" => Some(pos), + e => num(e), + }; + let supp = match info_value(c[7], "SUPP") { + "" => Some(0), + e => num(e), + }; + let (Some(end), Some(supp)) = (end, supp) else { + continue; + }; + by.entry(( + chrom_key(c[0]).to_owned(), + dup_or(info_value(c[7], "SVTYPE")), + )) + .or_default() + .push((pos, end, supp)); + } + rows.iter() + .map(|r| { + let (Some(Ok(start)), Some(Ok(end))) = ( + cell(r, "SV_start").map(crate::int), + cell(r, "SV_end").map(crate::int), + ) else { + return Some(1); + }; + let key = ( + chrom_key(cell(r, "SV_chrom").unwrap_or("")).to_owned(), + dup_or(cell(r, "SV_type").unwrap_or("")), + ); + let mut best: Option<(i128, i128)> = None; + for &(p, e, supp) in by.get(&key).map_or(&[][..], Vec::as_slice) { + let (dp, de) = ((p - start).abs(), (e - end).abs()); + if dp <= CONSENSUS_BP && de <= CONSENSUS_BP && best.is_none_or(|(d, _)| dp + de < d) + { + best = Some((dp + de, supp)); + } + } + Some(best.map_or(1, |(_, s)| s)) + }) + .collect() +} + +/// What step 35 wrote that the rule reads. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Step35 { + /// The interval table's deletions and duplications; `None` without it. + pub intervals: Option>>, + /// The HBA table (metric to value); `None` without it. + pub hba: Option>, + /// Parascopy's resolved SMN1 copy number at exon 7. + pub smn1: Option, +} + +/// SMN1 exon 7, c.840, GRCh38 1-based. +pub const SMN1_EXON7: (&str, i128) = ("chr5", 70951946); +/// The tables in paralogs/ that are not Parascopy catalog tables. +const NOT_CATALOG_TABLES: [&str; 2] = ["hba", "paraphase"]; +const PARASCOPY_COLUMNS: [&str; 9] = [ + "chrom", + "start", + "end", + "agCN", + "agCN_filter", + "agCN_qual", + "psCN", + "psCN_qual", + "psCN_filter", +]; +const INTERVAL_COLUMNS: [&str; 9] = [ + "locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change", +]; +const DOSAGE_COLUMNS: [&str; 6] = ["locus", "status", "copies", "estimate", "reason", "reading"]; + +fn has_all(path: &str, cols: &[&str]) -> Result { + let have = columns(path).map_err(|e| e.text())?; + Ok(cols.iter().all(|c| have.iter().any(|h| h == c))) +} + +fn short_row(path: &str, rows: &[DictRow], cols: &[&str]) -> Result<(), String> { + if let Some(i) = rows + .iter() + .position(|r| cols.iter().any(|k| cell(r, k).is_none())) + { + return Err(format!( + "{} line {} lacks a field", + crate::basename(path), + i + 2 + )); + } + Ok(()) +} + +/// Step 35's HBA or CYP21A2 table as metric to value; a table without a +/// status row reads as not assessed. +fn read_hba(path: &str) -> Result, String> { + let mut out = IndexMap::new(); + for l in text::lines(path).map_err(|e| e.text())? { + if let Some((k, v)) = l.split_once('\t') + && k != "metric" + { + out.insert(k.to_owned(), v.to_owned()); + } + } + if !matches!( + out.get("status").map(String::as_str), + Some("called" | "not assessed") + ) { + out = IndexMap::from([ + ("status".to_owned(), "not assessed".to_owned()), + ( + "message".to_owned(), + format!("not assessed: {} has no status row", crate::basename(path)), + ), + ]); + } + Ok(out) +} + +/// The psCN of SMN1 at exon 7 when Parascopy resolved it (filter PASS, +/// quality 20 or more, a whole number), as `parascopy_cn_at` reads it. +fn smn1_cn(rows: &[DictRow]) -> Option { + let (chrom, pos) = SMN1_EXON7; + for r in rows { + let (Some(c), Some(Ok(s)), Some(Ok(e))) = ( + cell(r, "chrom"), + cell(r, "start").map(crate::int), + cell(r, "end").map(crate::int), + ) else { + continue; + }; + let mut regions = vec![(c.to_owned(), s + 1, e)]; + for h in text::cell(r, "homologous_regions") + .filter(|v| !v.is_empty()) + .unwrap_or("*") + .split(',') + { + let h = pycompat::text::strip(h, None); + if let Some((hc, rest)) = h.split_once(':') { + let a: String = rest.chars().take_while(char::is_ascii_digit).collect(); + let rest = &rest[a.len()..]; + if let Some(rest) = rest.strip_prefix('-') { + let b: String = rest.chars().take_while(char::is_ascii_digit).collect(); + if !hc.is_empty() && !a.is_empty() && !b.is_empty() { + regions.push((hc.to_owned(), crate::int(&a).ok()?, crate::int(&b).ok()?)); + } + } + } + } + let Some(idx) = regions + .iter() + .position(|(c, b, e)| c == chrom && *b <= pos && pos <= *e) + else { + continue; + }; + let cns: Vec<&str> = text::cell_or(r, "psCN").split(',').collect(); + let quals: Vec<&str> = text::cell_or(r, "psCN_qual").split(',').collect(); + let cn = cns.get(idx).copied().unwrap_or("?"); + let q = quals + .get(idx) + .and_then(|q| pycompat::num::parse_float(q).ok()); + let filt = cell(r, "psCN_filter") + .filter(|v| !v.is_empty()) + .unwrap_or("."); + let ok = filt == "PASS" && q.is_some_and(|q| q >= 20.0) && isdigit(cn); + return if ok { crate::int(cn).ok() } else { None }; + } + None +} + +/// What the rule takes from step 35's tables, in the order the Parascopy +/// section reads them, so the same fault stops it. `Ok(None)`: step 35 +/// wrote no table. +pub fn read_step35(dir: &str, s: &str, assets: &str) -> Result, String> { + let paralogs = crate::join(dir, "paralogs"); + let mut tables: Vec<(String, String)> = Vec::new(); + if let Ok(entries) = std::fs::read_dir(¶logs) { + for e in entries.flatten() { + let Ok(name) = e.file_name().into_string() else { + continue; + }; + let Some(t) = name + .strip_prefix(&format!("{s}_")) + .and_then(|n| n.strip_suffix("_copy_number.tsv")) + else { + continue; + }; + let p = crate::join(¶logs, &name); + if !NOT_CATALOG_TABLES.contains(&t) && crate::is_file(&p) { + tables.push((t.to_owned(), p)); + } + } + } + tables.sort_by(|a, b| pycompat::mtime::py_cmp(&a.1, &b.1)); + let side = |n: &str| crate::join(¶logs, &format!("{s}_{n}")); + let (hp, ip, cp, dp) = ( + side("hba_copy_number.tsv"), + side("cn_intervals.tsv"), + side("cyp21a2_psv.tsv"), + side("cn_dosage.tsv"), + ); + let smn = tables + .iter() + .find(|(t, _)| t == "smn") + .map(|(_, p)| p.clone()); + let any = + smn.is_some() || !tables.is_empty() || [&hp, &ip, &cp, &dp].iter().any(crate::is_file); + if !any { + return Ok(None); + } + let mut smn1 = None; + if let Some(p) = &smn { + let r = rows(p).map_err(|e| e.text())?; + if !has_all(p, &PARASCOPY_COLUMNS)? { + return Err(format!( + "{} lacks the chrom, start, end, agCN, agCN_filter, agCN_qual, psCN, psCN_qual or \ + psCN_filter column", + crate::basename(p) + )); + } + smn1 = smn1_cn(&r); + } + let others: Vec<&(String, String)> = tables.iter().filter(|(t, _)| t != "smn").collect(); + if !others.is_empty() { + let cat = super::assets::asset_rows(assets, "paralog_loci.tsv")?; + for r in &cat { + if r.get("table").is_none() { + return Err("'table'".into()); + } + } + } + let hba = if crate::is_file(&hp) { + Some(read_hba(&hp)?) + } else { + None + }; + let intervals = if crate::is_file(&ip) { + let r = rows(&ip).map_err(|e| e.text())?; + if !has_all(&ip, &INTERVAL_COLUMNS)? { + return Err(format!( + "{} lacks the locus, name, status, copies, expected or change column", + crate::basename(&ip) + )); + } + short_row(&ip, &r, &INTERVAL_COLUMNS)?; + Some( + r.iter() + .filter(|x| matches!(cell(x, "change"), Some("deletion" | "duplication"))) + .map(|x| { + [ + "locus", "name", "chrom", "start", "end", "copies", "expected", "change", + ] + .iter() + .map(|k| ((*k).to_owned(), text::cell_or(x, k).to_owned())) + .collect() + }) + .collect(), + ) + } else { + None + }; + if crate::is_file(&cp) { + read_hba(&cp)?; + } + if crate::is_file(&dp) { + if !has_all(&dp, &DOSAGE_COLUMNS)? { + return Err(format!( + "{} lacks the {} or reading column", + crate::basename(&dp), + DOSAGE_COLUMNS[..5].join(", ") + )); + } + let r = rows(&dp).map_err(|e| e.text())?; + short_row(&dp, &r, &DOSAGE_COLUMNS)?; + } + for (_, p) in others { + rows(p).map_err(|e: ReadError| e.text())?; + columns(p).map_err(|e| e.text())?; + } + Ok(Some(Step35 { + intervals, + hba, + smn1, + })) +} + +/// The deletion inputs of the sample, read once for every rule that needs +/// them. +#[derive(Debug, Clone)] +pub struct DeletionInputs { + pub annotsv: Result>, String>, + pub step35: Result, String>, +} + +/// The genes a called deletion reaches: the deletions of each gene, the +/// engines read, and the engines not read with why. +pub fn gene_deletions( + inp: &DeletionInputs, + spans: &Spans, +) -> (HashMap>, Vec, Vec) { + let mut out: HashMap> = HashMap::new(); + let (mut read, mut unread) = (Vec::new(), Vec::new()); + let over = |chrom: &str, a: i128, b: i128| -> Vec { + spans + .iter() + .filter(|(_, (c, gs, ge))| c == chrom && *gs <= b && *ge >= a) + .map(|(g, _)| g.clone()) + .collect() + }; + let mut add = |g: String, x: Deletion| out.entry(g).or_default().push(x); + match &inp.annotsv { + Err(e) => unread.push(format!("AnnotSV (step 05): {e}")), + Ok(None) => unread.push("AnnotSV (step 05): no table".into()), + Ok(Some(rows)) => { + read.push("AnnotSV class 4 and 5 deletions (step 05)".into()); + for x in rows { + let (chrom, span) = x.sv.split_once(':').unwrap_or((&x.sv, "")); + let (a, b) = span.split_once('-').unwrap_or((span, "")); + if x.sv_type.as_deref() != Some("DEL") || x.common || !(isdigit(a) && isdigit(b)) { + continue; + } + // A call genotyped 0/0 is no deletion in this sample. + if let Some(gt) = x.gt.as_deref() + && gt.replace('|', "/").split('/').all(|a| a == "0") + { + continue; + } + let chrom = if chrom.starts_with("chr") { + chrom.to_owned() + } else { + format!("chr{chrom}") + }; + let (Ok(a_n), Ok(b_n)) = (crate::int(a), crate::int(b)) else { + continue; + }; + // SV_start is the base before the deleted bases. + let first = a_n + 1; + let mut genes: BTreeSet = over(&chrom, first, b_n).into_iter().collect(); + genes.extend(x.genes.iter().filter(|g| !spans.contains_key(*g)).cloned()); + for g in genes { + add( + g, + Deletion { + chrom: chrom.clone(), + start: first, + end: b_n, + inside: true, + text: format!( + "DEL {chrom}:{a}-{b} (AnnotSV class {}, {})", + x.class(), + x.caller.as_deref().unwrap_or("None") + ), + gt: x.gt.clone(), + ps: x.ps.clone(), + }, + ); + } + } + } + } + let st = match &inp.step35 { + Ok(Some(st)) => st, + Ok(None) => { + unread.push("step 35's copy-number tables: none".into()); + return (out, read, unread); + } + Err(e) => { + unread.push(format!("step 35's copy-number tables: {e}")); + return (out, read, unread); + } + }; + match &st.intervals { + Some(rows) => { + read.push("step 35's interval table (read depth)".into()); + for r in rows { + let f = |k: &str| r.get(k).map_or("", String::as_str); + if f("change") != "deletion" || !(isdigit(f("start")) && isdigit(f("end"))) { + continue; + } + let (Ok(a), Ok(b)) = (crate::int(f("start")), crate::int(f("end"))) else { + continue; + }; + for g in over(f("chrom"), a + 1, b) { + add( + g, + Deletion { + chrom: f("chrom").to_owned(), + start: a + 1, + end: b, + inside: true, + gt: None, + ps: None, + text: format!( + "{} deletion, {} of {} copies (read depth, step 35)", + f("locus"), + f("copies"), + f("expected") + ), + }, + ); + } + } + } + None => unread.push("step 35's interval table: none".into()), + } + let hba = st.hba.clone().unwrap_or_default(); + let n = hba.get("alpha_globin_genes").map_or("", String::as_str); + if hba.get("status").map(String::as_str) == Some("called") && isdigit(n) { + read.push("step 35's HBA1/HBA2 count".into()); + // 0 genes: no alpha-globin gene is left to carry a variant. + if digits(Some(n)).is_some_and(|k| (1..4).contains(&k)) { + for g in ["HBA1", "HBA2"] { + if let Some((c, a, b)) = spans.get(g) { + add( + g.to_owned(), + Deletion { + chrom: c.clone(), + start: *a, + end: *b, + inside: false, + gt: None, + ps: None, + text: format!("{n} of 4 alpha-globin genes (read depth, step 35)"), + }, + ); + } + } + } + } else { + unread.push(format!( + "step 35's HBA1/HBA2 count: {}", + hba.get("message") + .filter(|m| !m.is_empty()) + .map_or("not called", String::as_str) + )); + } + match st.smn1 { + Some(cn) => { + read.push("Parascopy's SMN1 count (step 35)".into()); + if cn == 1 + && let Some((c, a, b)) = spans.get("SMN1") + { + add( + "SMN1".into(), + Deletion { + chrom: c.clone(), + start: *a, + end: *b, + inside: true, + gt: None, + ps: None, + text: "one SMN1 copy (Parascopy, step 35)".into(), + }, + ); + } + } + None => unread.push("Parascopy's SMN1 count: not resolved".into()), + } + (out, read, unread) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn chromosome_prefix_is_safe_for_non_ascii_contigs() { + for (input, expected) in [ + ("éé", "éé"), + ("é", "é"), + ("☃1", "☃1"), + ("chr☃", "☃"), + ("CHR1", "1"), + ("cHrX", "X"), + ("chr", ""), + ("ch", "ch"), + ("", ""), + ] { + assert_eq!(chrom_key(input), expected); + } + } + + fn del(gt: Option<&str>, ps: Option<&str>) -> Deletion { + Deletion { + chrom: "chr1".into(), + start: 100, + end: 200, + text: "DEL".into(), + inside: true, + gt: gt.map(str::to_owned), + ps: ps.map(str::to_owned), + } + } + + fn var(pos: i128, gt: &str, ps: Option<&str>) -> Variant { + let (copies, ploidy) = super::super::clinvar::allele_copies(gt); + Variant { + gene: "G".into(), + key: ("chr1".into(), pos, "A".into(), "T".into()), + chrom: "chr1".into(), + pos, + text: format!("chr1:{pos} A>T"), + genotype: "het".into(), + copies, + ploidy, + gt: gt.into(), + ps: ps.map(str::to_owned), + allele: None, + } + } + + #[test] + fn phase_reads_only_called_phased_copies() { + assert_eq!(phase_haps(Some("0|1")), Some(BTreeSet::from([1]))); + assert_eq!(phase_haps(Some("1|0")), Some(BTreeSet::from([0]))); + assert_eq!(phase_haps(Some("0/1")), None); + assert_eq!(phase_haps(Some(".|1")), None); + assert_eq!(phase_haps(None), None); + } + + #[test] + fn read_backed_phase_decides_outside_the_deletion() { + // Outside the deletion's span: phase alone decides. + let trans = pair_phase(&del(Some("0|1"), Some("7")), &var(500, "1|0", Some("7"))); + assert_eq!(trans.0, Verdict::Biallelic); + assert!(trans.1.starts_with("in trans")); + let cis = pair_phase(&del(Some("0|1"), Some("7")), &var(500, "0|1", Some("7"))); + assert_eq!(cis.0, Verdict::Cis); + // Another phase set links nothing. + let other = pair_phase(&del(Some("0|1"), Some("7")), &var(500, "1|0", Some("8"))); + assert_eq!(other.0, Verdict::PossiblyBiallelic); + assert!(other.1.starts_with("phase unknown")); + } + + #[test] + fn a_phased_cis_pair_is_no_second_hit_even_inside_the_deletion() { + let v = var(150, "0|1", Some("7")); + let cis = pair_phase(&del(Some("0|1"), Some("7")), &v); + assert_eq!(cis.0, Verdict::Cis); + let found = two_hits( + &[v], + &HashMap::from([("G".to_owned(), vec![del(Some("0|1"), Some("7"))])]), + ); + assert!(found.entries.is_empty()); + } + + #[test] + fn homozygous_inside_is_hemizygous() { + let v = var(150, "1/1", None); + let found = two_hits( + std::slice::from_ref(&v), + &HashMap::from([("G".to_owned(), vec![del(None, None)])]), + ); + let e = &found.entries["G"]; + assert_eq!(e.verdict, Verdict::Biallelic); + assert_eq!(e.genotype, "hemizygous"); + assert!(found.hemi.contains_key(&("G".to_owned(), v.key.clone()))); + } + + #[test] + fn the_best_pair_of_a_gene_wins() { + let v = var(500, "1|0", Some("7")); + let found = two_hits( + std::slice::from_ref(&v), + &HashMap::from([( + "G".to_owned(), + vec![del(None, None), del(Some("0|1"), Some("7"))], + )]), + ); + assert_eq!(found.entries["G"].verdict, Verdict::Biallelic); + assert!(found.pairs[&("G".to_owned(), v.key)].starts_with("biallelic with DEL: in trans")); + } +} diff --git a/rust/og/src/main.rs b/rust/og/src/main.rs index 35aecb25..c32cb0e6 100644 --- a/rust/og/src/main.rs +++ b/rust/og/src/main.rs @@ -50,7 +50,7 @@ const USAGE: &str = "usage: og [args] og prs format [args] og prs table [args] og mito annotate|deletions|copies [args] - og summary collect --only a [args] + og summary collect --only a|b [args] og chip prepare|report|check-table [args] og cn bed|call|intervals|check-bam|marker-regions|marker-read [args] og bloodtype report|check-markers|systems [args] diff --git a/tests/parity/cases.tsv b/tests/parity/cases.tsv index 0a9afd6b..191f8e4d 100644 --- a/tests/parity/cases.tsv +++ b/tests/parity/cases.tsv @@ -838,6 +838,20 @@ summary-collect empty gen:summary-collect.sh,summary-collect/oracle.sh --sample summary-collect declared-arg gen:summary-collect.sh,summary-collect/oracle.sh --sample S --sample-dir declared-status --assets assets --declared-sex female --out declared-arg.json --only a declared-arg.json json:paths=schema_version,sample,run,databases,manifest,sections.coverage,sections.callability,sections.sex_check,sections.sample_qc,sections.mosaic_sex_loss,sections.read_count,sections.variants PYTHON_IMAGE summary-collect no-assets gen:summary-collect.sh,summary-collect/oracle.sh --sample S --sample-dir full --assets nowhere --out no-assets.json --only a no-assets.json json:paths=schema_version,sample,run,databases,manifest,sections.coverage,sections.callability,sections.sex_check,sections.sample_qc,sections.mosaic_sex_loss,sections.read_count,sections.variants PYTHON_IMAGE summary-collect no-folder gen:summary-collect.sh,summary-collect/oracle.sh --sample S --sample-dir no-such-folder --assets assets --out no-folder.json --only a no-folder.json json:paths=schema_version,sample,run,databases,manifest,sections.coverage,sections.callability,sections.sex_check,sections.sample_qc,sections.mosaic_sex_loss,sections.read_count,sections.variants PYTHON_IMAGE +summary-collect cb-full gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-full --assets assets --out cb-full.json --only b cb-full.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-quoted-smn gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-quoted-smn --assets assets --out cb-quoted-smn.json --only b cb-quoted-smn.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-quoted-interval gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-quoted-interval --assets assets --out cb-quoted-interval.json --only b cb-quoted-interval.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-quoted-dosage gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-quoted-dosage --assets assets --out cb-quoted-dosage.json --only b cb-quoted-dosage.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-unicode-contig gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-unicode-contig --assets assets --out cb-unicode-contig.json --only b cb-unicode-contig.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-nf gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-nf --assets assets --out cb-nf.json --only b cb-nf.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-probable gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-probable --assets assets --out cb-probable.json --only b cb-probable.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-failed gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-failed --assets assets --out cb-failed.json --only b cb-failed.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-bad gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-bad --assets assets --out cb-bad.json --only b cb-bad.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-nfbad gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-nfbad --assets assets --out cb-nfbad.json --only b cb-nfbad.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-notbi gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-notbi --assets assets --out cb-notbi.json --only b cb-notbi.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-notgz gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-notgz --assets assets --out cb-notgz.json --only b cb-notgz.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-empty gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-empty --assets assets --out cb-empty.json --only b cb-empty.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE +summary-collect cb-no-assets gen:summary-clinical.sh,summary-collect/oracle.sh,summary-collect/cb_vep.vcf.gz,summary-collect/cb_vep.vcf.gz.tbi --sample S --sample-dir cb-full --assets nowhere --out cb-no-assets.json --only b cb-no-assets.json json:paths=schema_version,sample,run,databases,manifest,sections.clinvar,sections.carrier,sections.pharmcat,sections.cpic,sections.pypgx,sections.cyrius,sections.cpsr,sections.clinical,sections.slivar,cyp2d6 PYTHON_IMAGE y-tree-positions small gen:y-inputs.sh --tree-dir t pos.txt pos.txt bytes PYTHON_IMAGE y-tree-positions big gen:y-inputs.sh --tree-dir ./big/ pos.txt pos.txt bytes PYTHON_IMAGE y-tree-positions odd gen:y-inputs.sh --tree-dir odd pos.txt pos.txt bytes PYTHON_IMAGE diff --git a/tests/parity/gen/summary-clinical.sh b/tests/parity/gen/summary-clinical.sh new file mode 100644 index 00000000..6d5cfbac --- /dev/null +++ b/tests/parity/gen/summary-clinical.sh @@ -0,0 +1,262 @@ +#!/usr/bin/env bash +# summary-clinical.sh: writes the sample folders of the summary-collect +# parity cases of group b (og summary collect --only b) in the current +# directory, the same every time. The clinical sections: ClinVar hits with +# each review status, a call floor, combined alleles, half calls and +# phased calls, the cDNA check, scan and deep intronic tables, gnomAD AF +# from an indexed VEP output; the carrier list; PharmCAT, CPIC, pypgx and +# Cyrius and step 36's CYP2D6 row; CPSR with its own class beside +# ClinVar's; the clinical filter's table and VCF with the ACMG SF tier; +# slivar's table and VCF with the phased compound-het pairs; and the +# deletions the two-hit rule pairs with a variant (AnnotSV, step 35's +# interval, HBA and SMN1 tables), in trans, in cis and inside the deletion. +# +# The VEP output is the committed pair cb_vep.vcf.gz and its .tbi, staged +# beside the folders; vep/S_vep.vcf.gz links to it. assets/ is a copy of +# the repository's gene lists. pymods/ holds the helper modules the Python +# collector imports, at a pinned commit, so the oracle reads them instead +# of stubs. scripts/ci/oracle-diff.sh runs it. +set -euo pipefail +export LC_ALL=C TZ=UTC +here=$(cd "$(dirname "$0")" && pwd) +root=$(cd "${here}/../../.." && pwd) + +put() { mkdir -p "$(dirname "$1")"; printf '%b' "$2" > "$1"; } +putgz() { mkdir -p "$(dirname "$1")"; printf '%b' "$2" | gzip -n -c > "$1"; } +at() { touch -t "$1" "${@:2}"; } # [[CC]YY]MMDDhhmm[.SS], in UTC + +T0=1577836800 # 2020-01-01 00:00:00 UTC, the start of the run +OLD=201906011230.00 # before the run +NEW=202103040506.07 # after it +NEWER=202103040507.00 + +# The helper modules of bin/collect_summary.py, as they were when the +# collector was last the summary's writer. +HELPERS=7c1470bc50933f48a1a0dde3c0e0b2b87bc48718 +git -C "$root" cat-file -e "${HELPERS}^{commit}" 2> /dev/null \ + || git -C "$root" fetch -q --depth=1 origin "$HELPERS" +mkdir -p pymods +for m in pgx_outside_calls clip_check annotsv_postfilter paraphase_report traits_report; do + git -C "$root" show "${HELPERS}:bin/${m}.py" > "pymods/${m}.py" +done + +mkdir -p assets +for a in acmg_sf.tsv carrier_genes.tsv paralog_loci.tsv; do + cp "${root}/assets/${a}" assets/ +done + +VCF_HEAD='##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n' +PL='criteria_provided,_multiple_submitters,_no_conflicts' + +# --- the ClinVar hits of step 06 --------------------------------------------------- +HITS="##fileformat=VCFv4.2\n##PGP_CALL_FLOOR=GQ=20,DP=8\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" +HITS+="chr1\t45331556\t1\tC\tT\t.\tPASS\tGENEINFO=MUTYH:4595;CLNSIG=Pathogenic;CLNREVSTAT=${PL}\tGT:GQ:DP:PS\t0|1:50:30:45331000\n" +HITS+="chr6\t26092913\t2\tG\tA\t.\tPASS\tGENEINFO=HFE:3077;CLNSIG=Pathogenic;CLNREVSTAT=reviewed_by_expert_panel\tGT:GQ:DP\t0/1:99:40\n" +HITS+="chr7\t117559590\t3\tATCT\tA\t.\tPASS\tGENEINFO=CFTR:1080;CLNSIG=Pathogenic;CLNREVSTAT=practice_guideline\tGT:GQ:DP\t0/1:15:30\n" +HITS+="chr13\t20189473\t4\tC\tT\t.\tPASS\tGENEINFO=GJB2:2706;CLNSIG=Likely_pathogenic;CLNREVSTAT=criteria_provided,_single_submitter;PGP_JOINED_FROM=chr13:20189473:C:T|chr13:20189475:G:A\tGT:GQ:DP\t1/.:40:5\n" +HITS+="chr3\t15645186\t5\tG\tC\t.\tPASS\tGENEINFO=BTD:686;CLNSIG=Pathogenic/Likely_pathogenic;CLNREVSTAT=criteria_provided,_conflicting_classifications\tGT\t0/1\n" +HITS+="chr17\t43057100\t6\tC\tT\t.\tPASS\tGENEINFO=BRCA1:672|NBR2:10230;CLNSIG=Pathogenic;CLNREVSTAT=reviewed_by_expert_panel;PGP_PART_OF=chr17:43057100:C:T|chr17:43057102:A:G\tGT:GQ:DP\t0/1:60:35\n" +HITS+="chrX\t101398000\t7\tC\tT\t.\tPASS\tGENEINFO=GLA:2717;CLNSIG=Pathogenic;CLNREVSTAT=no_assertion_criteria_provided;PGP_ONE_ALLELE=chrX:101398004:G:A\tGT\t1\n" +HITS+="chr16\t176700\t8\tG\tA\t.\tPASS\tGENEINFO=HBA1:3039;CLNSIG=Pathogenic;CLNREVSTAT=${PL}\tGT:GQ:DP\t0/1:45:31\n" +HITS+="chr5\t70951946\t9\tC\tT\t.\tPASS\tGENEINFO=SMN1:6606;CLNSIG=Pathogenic;CLNREVSTAT=${PL}\tGT:GQ:DP\t1/1:45:31\n" +HITS+="chr2\t100\t10\tA\tG\t.\tPASS\tCLNREVSTAT=;CLNSIG=\tGT:DP\t0/1:x\n" +HITS+="chr2\t200\t11\tA\tG\t.\tPASS\n" + +CDNA='chrom\tpos\tref\talt\talt_reads\talt_clipped\tref_reads\tref_clipped\tcdna_check\tcdna_note\tconstruct_site\n' +CDNA+='chr1\t45331556\tC\tT\t12\t6\t10\t0\tflagged\tpossible cDNA/transgene contamination: check the reads\tyes\n' +CDNA+='chr6\t26092913\tG\tA\t20\t0\t20\t0\tclean\t\t\n' +CDNA+='chr6\tx\tG\tA\t1\t1\t1\t1\tclean\t\t\n' +SCAN='# cdna_scan: 2 genes with spliced reads\nchrom\tstart\tend\tjunctions\treads\tjunction_list\nchr1\t100\t200\t2\t5\tA:1-2,A:3-4\nchr2\tx\t1\t1\t1\t\n' +DEEP='chrom\tpos\tref\talt\tgene\thgvs\tclinvar_id\tcall\tdetail\nchr7\t117639961\tC\tT\tCFTR\tc.3718-2477C>T\t7129\t\t\nchr7\tNA\tC\tT\tCFTR\t\t\t\t\n' + +# --- pharmacogenomics ---------------------------------------------------------------- +CPIC='Gene\tDiplotype\tPhenotype\tStatus\nCYP2D6\t*1/*4\tIntermediate Metabolizer\tnon-normal\nCYP2C19\t*1/*17\tRapid\tambiguous\nCYP2C9\t*1/*1\tNormal\tnormal\nDPYD\n' +CYRIUS_HEADER='Sample\tGenotype\tFilter\n' +REC='CPIC lookup\n WARNING: CYP2C19 has two possible phenotypes\nnothing\nWARNING:last\n' +PYPGX='Gene\tDiplotype\tPhenotype\tStatus\tSource\tCNV_call\tQuality\tMin_GQ\tWeakest_variant\tDefining_filter\n' +PYPGX+='CYP2D6\t*1/*2\t Normal Metabolizer \tok\tBAM\tNormal\tok\t\t\tPASS\n' +PYPGX+='CYP2A6\tFAILED\t\tfailed\tBAM\t\t\t\t\t\n' +PYPGX+='CYP2B6\t*1/*6\tIM\tdepth check failed\tbam\tN/A\tok\t\t\t\n' +PYPGX+='CYP2E1\t.\t\tnot callable on this reference\tBAM\t.\t\t\t\t\n' +PYPGX+='UGT1A1\tIndeterminate\t\tok\tVCF\t\tlow-confidence: check the reads\t 12 \tchr2:1:A:G\t\n' +PYPGX+='SLCO1B1\t*1/*5\t\tok\tBAM\t.\tlow-confidence: check the reads\tx\t\t\n' +DEPTH='metric\tvalue\nstatus\tok\ngene_mapq1_share\t0.99\nbody_mapq20_share\t0.95\nother\t1\n' +COMP='Gene\tMatch\nCYP2D6\tYes\nCYP2C19\tNo\nCYP3A5\tpypgx only\nTPMT\tPharmCAT only\n' +CONS='Gene\tResult\tOutside_call\tReason\tEvidence\nHLA-A\tx\tno\t\t\nCYP2D6\t*1/*2\tyes\tpypgx and Cyrius agree\tpypgx *1/*2; Cyrius *2/*1\n' + +# --- CPSR (step 17) -------------------------------------------------------------------- +CPSR='VAR_ID\tGENOMIC_CHANGE\tSYMBOL\tCLASSIFICATION\tCPSR_CLASSIFICATION\tCLINVAR_CLASSIFICATION\tASSERTION_AUTHORITY\tACMG_CODE\tDP_CONTROL\n' +CPSR+='17_43057100_C_T\t17:g.43057100C>T\tBRCA1\tPathogenic\tLikely_Benign\tPathogenic\tClinVar\tPVS1\t6\n' +CPSR+='\t13:g.32338000G>A\tBRCA2\tLikely_Pathogenic\tVUS\t\tCPSR\t\t40\n' +CPSR+='2_47403000_A_G\t\tMSH2\tVUS\tLP\tB\t\t\t\n' +CPSR+='x\tbad\tATM\tBenign\tBenign\t\t\t\t3\n' +CPSR+='5_112839000_A_T\t\tAPC\t\tPathogenic\t\t\t\t30\n' +CPSR_CDNA="${CDNA}17\t43057100\tC\tT\t5\t5\t9\t0\tflagged\tsplice\t\n" + +# --- the clinical filter (step 23) ------------------------------------------------------ +CSQF='Allele|Consequence|IMPACT|SYMBOL|Feature|CANONICAL|MANE_SELECT|gnomADe_AF|gnomADg_AF' +CLIN="##fileformat=VCFv4.2\n##PGP_CALL_FLOOR=GQ=20,DP=8,HET_AB=0.25-0.75,HOM_AB=0.85\n" +CLIN+="##INFO=\n" +CLIN+='#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n' +CLIN+='chr1\t45332000\t.\tG\tA\t50\tPASS\tCSQ=A|stop_gained|HIGH|MUTYH|T1|YES|NM_1|0.0001|\tGT:GQ:DP:AD:PS\t1|0:40:30:15,15:45331000\n' +CLIN+='chr2\t47414000\t.\tA\tAT\t50\tPASS\tCSQ=AT|frameshift_variant|HIGH|MSH2|T2|YES||0.00002|0.00001\tGT:GQ:DP:AD\t0/1:50:30:28,2\n' +CLIN+='chr10\t100\t.\tC\tT\t50\tPASS\tCSQ=T|stop_gained|HIGH|PTEN|T3|YES||0.2|\tGT\t0/1\n' +CLIN+='chr6\t26093000\t.\tC\tT\t50\tPASS\tCSQ=T|stop_gained|HIGH|HFE|T4|YES||&|\tGT:GQ:DP:AD\t1/1:12:5:3,4\n' +CLIN+='chr17\t43057100\t.\tC\tT\t50\tPASS\tCSQ=T|missense_variant|MODERATE|BRCA1|T5|YES|NM_2||,T|stop_gained|HIGH|BRCA1|T6|||0.0001|\tGT\t0/1\n' +CLIN+='chr11\t47332000\t.\tG\tA\t50\tPASS\tCSQ=A|splice_donor_variant&intron_variant|HIGH|MYBPC3|T7|||0.001|,A|stop_gained|HIGH|MYBPC3|T8|||0.0001|\tGT:GQ\t0/1:.\n' +CLINT='GENE\tIMPACT\tFREQ_TIER\tCALL_CHECK\nMUTYH\tHIGH\trare\tok\nMSH2\tHIGH\trare\tweak\nHFE\tHIGH\tunknown\tweak\n.\tMODERATE\trare\tok\n\t\t\t\n' + +# --- slivar (step 31) ------------------------------------------------------------------- +SLIV='CHROM\tPOS\tGENE\tFREQ_TIER\nchr1\t1\tA\trare\nchr1\t2\tA\tunknown\n' +CHETS='GENE\tCHROM\nA\tchr1\nA\tchr1\nB\tchr2\n\tchr3\n' +PAIRS='GENE\tCHROM\tPOS_1\tREF_1\tALT_1\tPOS_2\tREF_2\tALT_2\tPHASE\tTRANS_READS\tCIS_READS\n' +for i in $(seq 1 23); do PAIRS+="G${i}\tchr1\t${i}0\tA\tG\t${i}5\tC\tT\ttrans\t${i}\t1\n"; done +PAIRS+='H\tchr2\t1\tA\tG\t2\tC\tT\tcis\t0\t4\nK\tchr2\t3\tA\tG\nL\tchr2\t4\tA\tG\t5\tC\tT\tunknown\tx\t\n' + +# --- the deletions of the two-hit rule -------------------------------------------------- +ANN='SV_chrom\tSV_start\tSV_end\tSV_length\tSV_type\tGene_name\tACMG_class\tCaller\tExon_edge\tExon_edge_detail\tPopulation_AF\tPopulation_source\tCommon\tAnnotSV_class\tFORMAT\tS\tAnnotation_mode\n' +ANN+='1\t45338000\t45340000\t-2000\tDEL\tMUTYH\t5\tManta\tyes\t\t\t\tno\t5\tGT:PS\t1|0:45331000\tfull\n' +ANN+='chr2\t47386000\t47388000\t-2000\tDEL\tEPCAM\t4\tDelly\tno\t\t\t\tno\t4\tGT\t0/1\tfull\n' +ANN+='17\t43000000\t43100000\t100000\tDUP\tBRCA1;NBR2\t4\tCNVpytor\tno\t\t0.001\tgnomAD\tno\t4\tGT\t0/1\tfull\n' +ANN+='13\t20180000\t20200000\t-20000\tDEL\tGJB2;NOTAGENE\t5\tManta\tno\t\t\t\tno\t5\tGT\t0/0\tfull\n' +ANN+='13\t20180000\t20200000\t-20000\tDEL\tGJB2;XYZ1\t4\tManta\tno\t\t\t\tno\t4\tGT\t0/1\tfull\n' +ANN+='3\t15600000\t15700000\t-100000\tDEL\tBTD\t3\tManta\tno\t\t\t\tno\t5\tGT\t0/1\tfull\n' +ANN+='3\t15600000\t15700000\t-100000\tDEL\tBTD\t5\tManta\tno\t\t0.2\tgnomAD\tyes\t5\tGT\t0/1\tfull\n' +CONSV='##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tmanta\tdelly\tcnvpytor\n' +CONSV+='chr1\t45338010\t.\tN\t\t.\tPASS\tSVTYPE=DEL;END=45340020;SUPP=2;SUPP_VEC=110\n' +CONSV+='chr1\t45338500\t.\tN\t\t.\tPASS\tSVTYPE=DEL;END=45340900;SUPP=3;SUPP_VEC=111\n' +CONSV+='chr17\t43000000\t.\tN\t\t.\tPASS\tSVTYPE=DUP:TANDEM;END=43100000;SUPP=1;SUPP_VEC=0011\n' +SMN='chrom\tstart\tend\tagCN\tagCN_filter\tagCN_qual\tpsCN\tpsCN_qual\tpsCN_filter\thomologous_regions\n' +SMN+='chr5\t70925000\t70953000\t4\tPASS\t40\t1,3\t30,30\tPASS\tchr5:70049000-70078000\n' +HBA='metric\tvalue\nstatus\tcalled\nalpha_globin_genes\t3\n' +IV='locus\tname\tchrom\tstart\tend\tstatus\tcopies\texpected\tchange\n' +IV+='CFTR_del\tCFTR region\tchr7\t117500000\t117600000\tcalled\t1\t2\tdeletion\n' +IV+='DUP1\tx\tchr1\t1\t100\tcalled\t3\t2\tduplication\n' +IV+='GJB6\tx\tchr13\tx\t2\tcalled\t1\t2\tdeletion\n' + +# full: every input, the bash layout, a run whose steps were all ok. +full() { + local d=$1 + put "$d/clinvar/S_clinvar_hits.vcf" "$HITS" + put "$d/clinvar/S_clinvar_cdna_check.tsv" "$CDNA" + put "$d/clinvar/S_clinvar_cdna_scan.tsv" "$SCAN" + put "$d/clinvar/S_clinvar_deep_intronic.tsv" "$DEEP" + mkdir -p "$d/vep" + ln -s ../../cb_vep.vcf.gz "$d/vep/S_vep.vcf.gz" + ln -s ../../cb_vep.vcf.gz.tbi "$d/vep/S_vep.vcf.gz.tbi" + put "$d/pharmcat/S.report.json" '{"pharmcatVersion": "3.4.0", "genes": {}}' + put "$d/pharmcat/S.report.html" '' + put "$d/cpic/S_phenotypes.tsv" "$CPIC" + put "$d/cpic/S_cpic_recommendations.txt" "$REC" + put "$d/pypgx/S_pypgx_summary.tsv" "$PYPGX" + put "$d/pypgx/S_cyp2d6_depth_check.tsv" "$DEPTH" + put "$d/pypgx/S_pharmcat_comparison.tsv" "$COMP" + put "$d/cyrius/S_cyp2d6.tsv" "${CYRIUS_HEADER}S\t*2/*1\tPASS\n" + put "$d/pgx_consensus/S_pgx_consensus.tsv" "$CONS" + put "$d/cpsr/S.cpsr.grch38.html" '' + putgz "$d/cpsr/S.cpsr.grch38.classification.tsv.gz" "$CPSR" + put "$d/cpsr/S_cpsr_cdna_check.tsv" "$CPSR_CDNA" + put "$d/clinical/S_clinical_summary.tsv" "$CLINT" + putgz "$d/clinical/S_clinical.vcf.gz" "$CLIN" + put "$d/slivar/S_slivar_summary.tsv" "$SLIV" + put "$d/slivar/S_compound_hets.tsv" "$CHETS" + put "$d/slivar/S_compound_het_pairs.tsv" "$PAIRS" + put "$d/annotsv/S_sv_class4_5_checked.tsv" "$ANN" + putgz "$d/sv_merged/S_sv_consensus.vcf.gz" "$CONSV" + put "$d/paralogs/S_smn_copy_number.tsv" "$SMN" + put "$d/paralogs/S_hba_copy_number.tsv" "$HBA" + put "$d/paralogs/S_cn_intervals.tsv" "$IV" + put "$d/logs/run_status.tsv" "meta\tstarted_epoch\t${T0}\nstep\t06\tok\nstep\t07\tok\nstep\t17\tok\nstep\t21\tok\nstep\t23\tok\nstep\t27\tok\nstep\t31\tok\nstep\t32\tok\n" +} +full cb-full + +# Change one header at a time. Quoted TSV fields have the same names as +# the unquoted baseline, including the dosage table's validation columns. +full cb-quoted-smn +put cb-quoted-smn/paralogs/S_smn_copy_number.tsv "\"chrom\"${SMN#chrom}" +full cb-quoted-interval +put cb-quoted-interval/paralogs/S_cn_intervals.tsv "\"locus\"${IV#locus}" +full cb-quoted-dosage +put cb-quoted-dosage/paralogs/S_cn_dosage.tsv '"locus"\tstatus\tcopies\testimate\treason\treading\n' + +# A valid UTF-8 contig with a code point crossing byte offset three must +# leave the consensus readable. The last chr17 record still matches. +full cb-unicode-contig +putgz cb-unicode-contig/sv_merged/S_sv_consensus.vcf.gz "${CONSV//chr1\\t/éé\\t}" + +# nf: the Nextflow names. The clinical VCF and slivar's VCF without their +# tables, two PharmCAT reports (the newer wins), Cyrius under another name, +# CPSR by its wildcard names, no VEP output, a consensus older than pypgx, +# PharmCAT's step failed, CPIC from an earlier run, pypgx failed. +put cb-nf/clinvar/S_clinvar_hits.vcf "$HITS" +putgz cb-nf/clinical/S_clinical.vcf.gz "$CLIN" +putgz cb-nf/slivar/S_prioritized.vcf.gz "${VCF_HEAD}chr1\t1\t.\tA\tG\t1\tPASS\t.\tGT\t0/1\nchr1\t2\t.\tA\tG\t1\tPASS\t.\tGT\t0/1\n" +put cb-nf/slivar/S_compound_het_pairs.tsv "$PAIRS" +put cb-nf/vcf/A.report.json '{"pharmcatVersion": 3.2}' +put cb-nf/vcf/B.report.json '[1, 2]' +put cb-nf/cpic/S_phenotypes.tsv 'Gene\tDiplotype\tPhenotype\tStatus\nCYP2D6\t*4/*4\tPoor\tnot called\n' +put cb-nf/cyrius/T_cyp2d6.tsv "${CYRIUS_HEADER}T\t*1/*2\tLowQ\n" +put cb-nf/pypgx/S_pypgx_summary.tsv "$PYPGX" +put cb-nf/pgx_consensus/S_pgx_consensus.tsv "$CONS" +put cb-nf/cpsr/T.cpsr.grch37.html '' +putgz cb-nf/cpsr/T.cpsr.grch37.classification.tsv.gz 'SYMBOL\tFINAL_CLASSIFICATION\tCPSR_CLASSIFICATION\nA\tPathogenic\tBenign\nB\t\t\n' +put cb-nf/annotsv/S_sv_class4_5_checked.tsv '' +put cb-nf/paralogs/S_cn_intervals.tsv "$IV" +put cb-nf/logs/run_status.tsv "meta\tstarted_epoch\t${T0}\nstep\t07\tfailed (exit 1)\nstep\t27\tskipped\nstep\t32\tok\nstep\t21\tok\n" + +# probable: pypgx alone, Cyrius not usable, the depth check ok; no consensus. +put cb-probable/pypgx/S_pypgx_summary.tsv "$PYPGX" +put cb-probable/pypgx/S_cyp2d6_depth_check.tsv "$DEPTH" +put cb-probable/cyrius/S_cyp2d6.tsv "${CYRIUS_HEADER}S\t*1/*3\tLowQ\n" +put cb-probable/cpic/S_phenotypes.tsv 'Gene\tDiplotype\tPhenotype\tStatus\nCYP2C9\t*1/*1\tNormal\tnormal\n' +put cb-probable/cpic/S_cpic_recommendations.txt 'PARSING FAILED for one gene\n' + +# failed: CPIC, pypgx and Cyrius failed in the run and wrote nothing. +put cb-failed/logs/run_status.tsv "meta\tstarted_utc\t2020-02-02 02:02 UTC\nstep\t21\tfailed (exit 1)\nstep\t27\tfailed\nstep\t32\tfailed (exit 3)\nstep\t06\tok\n" +put cb-failed/pgx_consensus/S_pgx_consensus.tsv 'Gene\tResult\nCYP2C19\tx\n' + +# bad: a hit whose POS is no number, a report that is no JSON, a CPSR +# table cut short, pypgx without a Gene column, an AnnotSV table without +# its columns, a Parascopy table without its columns, a clinical floor that +# is no number, slivar's VCF cut short. +put cb-bad/clinvar/S_clinvar_hits.vcf "${VCF_HEAD}chr1\tx\t.\tA\tG\t.\tPASS\tGENEINFO=MUTYH:1;CLNSIG=Pathogenic\tGT\t0/1\n" +put cb-bad/pharmcat/S.report.json '{"pharmcatVersion": ' +putgz cb-bad/cpsr/S.cpsr.grch38.classification.tsv.gz "$CPSR" +head -c 60 cb-bad/cpsr/S.cpsr.grch38.classification.tsv.gz > cb-bad.tmp && mv cb-bad.tmp cb-bad/cpsr/S.cpsr.grch38.classification.tsv.gz +put cb-bad/pypgx/S_pypgx_summary.tsv 'Name\tSource\tCNV_call\nCYP2D6\tBAM\tNormal\n' +put cb-bad/annotsv/S_sv_class4_5_checked.tsv 'SV_chrom\tSV_start\n1\t1\n' +put cb-bad/paralogs/S_smn_copy_number.tsv 'chrom\tstart\tend\nchr5\t1\t2\n' +putgz cb-bad/clinical/S_clinical.vcf.gz "##PGP_CALL_FLOOR=GQ=20,HOM_AB=0..8\n#CHROM\n" +put cb-bad/clinical/S_clinical_summary.tsv 'GENE\tIMPACT\nA\tHIGH\n' +put cb-bad/cyrius/S_cyp2d6.tsv '' +put cb-bad/cpic/S_phenotypes.tsv 'Gene\tDiplotype\n' +putgz cb-bad/slivar/S_prioritized.vcf.gz "$CLIN" +head -c 100 cb-bad/slivar/S_prioritized.vcf.gz > cb-bad.tmp && mv cb-bad.tmp cb-bad/slivar/S_prioritized.vcf.gz + +# nf-bad: the Nextflow clinical VCF cut short; Parascopy's catalog tables +# without their catalog; an HBA table without a status. +put cb-nfbad/clinvar/S_clinvar_hits.vcf "$HITS" +putgz cb-nfbad/clinical/S_clinical.vcf.gz "$CLIN" +head -c 100 cb-nfbad/clinical/S_clinical.vcf.gz > cb-nfbad.tmp && mv cb-nfbad.tmp cb-nfbad/clinical/S_clinical.vcf.gz +put cb-nfbad/paralogs/S_hba_copy_number.tsv 'metric\tvalue\nalpha_globin_genes\t2\n' +put cb-nfbad/paralogs/S_rhd_copy_number.tsv "$SMN" + +# notbi: the VEP output without its index; the HBA count of no gene left. +put cb-notbi/clinvar/S_clinvar_hits.vcf "$HITS" +mkdir -p cb-notbi/vep +ln -s ../../cb_vep.vcf.gz cb-notbi/vep/S_vep.vcf.gz +put cb-notbi/paralogs/S_hba_copy_number.tsv 'metric\tvalue\nstatus\tcalled\nalpha_globin_genes\t0\n' + +# notgz: the clinical and slivar VCFs as plain text under a .gz name. The +# Python tells gzip by the name and finds both unreadable; og tells it by +# the first bytes and reads them (a golden). +put cb-notgz/clinical/S_clinical.vcf.gz "$CLIN" +put cb-notgz/slivar/S_prioritized.vcf.gz "${VCF_HEAD}chr1\t1\t.\tA\tG\t1\tPASS\t.\tGT\t0/1\n" + +mkdir -p cb-empty + +# Every file after the run started; the CPIC table of nf from before it; +# nf's second PharmCAT report and its consensus' sources newer. +find . -type f -exec touch -t "$NEW" {} + +at "$OLD" cb-nf/cpic/S_phenotypes.tsv cb-nf/pgx_consensus/S_pgx_consensus.tsv +at "$NEWER" cb-nf/vcf/B.report.json diff --git a/tests/parity/golden/summary-collect/cb-notgz/cb-notgz.json b/tests/parity/golden/summary-collect/cb-notgz/cb-notgz.json new file mode 100644 index 00000000..e972fd82 --- /dev/null +++ b/tests/parity/golden/summary-collect/cb-notgz/cb-notgz.json @@ -0,0 +1,180 @@ +{ + "schema_version": 2, + "sample": "S", + "generated_utc": "2026-11-13 11:46 UTC", + "run": { + "status_file": false, + "started_utc": null, + "declared_sex": null + }, + "databases": { + "clinvar_release": null, + "cpsr_clinvar_release": null, + "hla_database": null + }, + "manifest": [], + "sections": { + "clinvar": { + "title": "ClinVar screen", + "step": "06", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "carrier": { + "title": "Carrier list", + "step": "06", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "pharmcat": { + "title": "PharmCAT", + "step": "07", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "cpic": { + "title": "CPIC drug-gene recommendations", + "step": "27", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "pypgx": { + "title": "pypgx", + "step": "32", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "cyrius": { + "title": "Cyrius (CYP2D6)", + "step": "21", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "cpsr": { + "title": "Cancer predisposition (CPSR)", + "step": "17", + "state": "missing", + "source": null, + "file_date": null, + "note": null, + "data": {} + }, + "clinical": { + "title": "Clinical variant filter", + "step": "23", + "state": "ok", + "source": "clinical/S_clinical.vcf.gz", + "file_date": "2021-03-04 05:06 UTC", + "note": null, + "data": { + "variants": 6, + "acmg_sf": { + "version": "ACMG SF v3.3", + "genes_on_list": 84, + "rules": "assets/acmg_sf.tsv", + "clinvar_hits": [], + "clinvar_checked": false, + "high_impact": [ + { + "gene": "MYBPC3", + "variant": "chr11:47332000 G>A", + "genotype": "het", + "consequence": "splice_donor_variant,intron_variant", + "transcript": "ANY", + "frequency": "rare", + "frequency_unknown": false, + "weak_call": false, + "call_note": null, + "gq": null, + "dp": null, + "ad": null, + "rule": "All P and LP", + "rule_note": null, + "severity": null + }, + { + "gene": "MSH2", + "variant": "chr2:47414000 A>AT", + "genotype": "het", + "consequence": "frameshift_variant", + "transcript": "CANONICAL", + "frequency": "rare", + "frequency_unknown": false, + "weak_call": true, + "call_note": "alt reads 0.07 of the depth on a heterozygote", + "gq": 50, + "dp": 30, + "ad": "28,2", + "rule": "All P and LP", + "rule_note": "A deletion of the 3' end of EPCAM silences MSH2 and counts as an MSH2 finding.", + "severity": null + } + ], + "high_impact_checked": true, + "high_impact_rare": 1, + "high_impact_weak": 1, + "high_impact_frequency_unknown": 0, + "structural": [], + "structural_checked": false, + "structural_note": null, + "not_reported": [ + { + "gene": "MUTYH", + "variant": "chr1:45332000 G>A", + "genotype": "het", + "source": "HIGH impact", + "reason": "one heterozygous variant; MUTYH is reported with P and LP (2 variants)" + }, + { + "gene": "HFE", + "variant": "chr6:26093000 C>T", + "genotype": "hom", + "source": "HIGH impact", + "reason": "HFE is reported only for: HFE p.C282Y homozygotes only" + } + ], + "two_hit_note": null + } + } + }, + "slivar": { + "title": "Variant prioritization (slivar)", + "step": "31", + "state": "ok", + "source": "slivar/S_prioritized.vcf.gz", + "file_date": "2021-03-04 05:06 UTC", + "note": null, + "data": { + "prioritized": 1 + } + } + }, + "cyp2d6": { + "calls": { + "PharmCAT": null, + "pypgx": null, + "Cyrius": null + }, + "agree": null, + "consensus": null, + "probable": null + } +} diff --git a/tests/parity/inputs/summary-collect/cb_vep.vcf.gz b/tests/parity/inputs/summary-collect/cb_vep.vcf.gz new file mode 100644 index 00000000..4a42ab23 Binary files /dev/null and b/tests/parity/inputs/summary-collect/cb_vep.vcf.gz differ diff --git a/tests/parity/inputs/summary-collect/cb_vep.vcf.gz.tbi b/tests/parity/inputs/summary-collect/cb_vep.vcf.gz.tbi new file mode 100644 index 00000000..1036d9b0 Binary files /dev/null and b/tests/parity/inputs/summary-collect/cb_vep.vcf.gz.tbi differ diff --git a/tests/parity/inputs/summary-collect/oracle.sh b/tests/parity/inputs/summary-collect/oracle.sh index 09d629c4..1a3e15a7 100644 --- a/tests/parity/inputs/summary-collect/oracle.sh +++ b/tests/parity/inputs/summary-collect/oracle.sh @@ -5,13 +5,16 @@ # --assets becomes ASSETS_DIR, the global render_report.py --assets sets; # --only is og's alone (the Python writes every section, and the case # compares the group's subtrees). The modules collect_summary imports are -# stubbed: none is reached for a sample folder that holds only the -# sections the cases cover, except run_delta's two readers of history/ and -# recheck/, which find nothing here. +# stubbed, except those a generator staged in pymods/ (the clinical cases +# need the real helpers); no other one is reached for a sample folder that +# holds only the sections the cases cover, except run_delta's two readers +# of history/ and recheck/, which find nothing here. prog=$1 shift d=$(mktemp -d) -sed -n 's/^import \([a-z_]*\) # noqa.*/\1/p' "$prog" | while read -r m; do : > "$d/$m.py"; done +sed -n 's/^import \([a-z_]*\) # noqa.*/\1/p' "$prog" | while read -r m; do + if [ -f "pymods/$m.py" ]; then cp "pymods/$m.py" "$d/"; else : > "$d/$m.py"; fi +done printf 'def read_recheck(path):\n return None\n\n\ndef newest_history(sample_dir):\n return None\n' > "$d/run_delta.py" cp "$prog" "$d/collect_summary.py" PYTHONPATH=$d exec python3 -c '