From 0f241f313473b599840d3e17bcd02ef4c9795764 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 12:27:46 +0100 Subject: [PATCH 1/7] feat(dev): og dev build rebuilds the three committed asset tables og dev build cn-intervals, callability and sv-rule-regions read the pinned downloads from a folder and write bin/cn_loci.tsv's interval rows, assets/callability_genes.tsv and assets/sv_rule_regions.tsv, byte for byte as the Python builders did. The GTF is read as bytes with a hand scanner of its attributes, the cut sets are merged once per chromosome and searched by binary search, and og never downloads. The SV rule builder now keeps the first GENCODE record per HGNC id and skips _PAR_Y copies, as the callability builder does; the committed table does not change. Bead: pgp-bs9u.35 --- rust/Cargo.lock | 1 + rust/crates/og-dev/Cargo.toml | 3 +- rust/crates/og-dev/src/build/callability.rs | 628 ++++++++++ rust/crates/og-dev/src/build/cn_intervals.rs | 1162 ++++++++++++++++++ rust/crates/og-dev/src/build/gtf.rs | 133 ++ rust/crates/og-dev/src/build/mod.rs | 298 +++++ rust/crates/og-dev/src/build/sv_rules.rs | 649 ++++++++++ rust/crates/og-dev/src/lib.rs | 6 +- rust/og/src/main.rs | 1 + 9 files changed, 2879 insertions(+), 2 deletions(-) create mode 100644 rust/crates/og-dev/src/build/callability.rs create mode 100644 rust/crates/og-dev/src/build/cn_intervals.rs create mode 100644 rust/crates/og-dev/src/build/gtf.rs create mode 100644 rust/crates/og-dev/src/build/mod.rs create mode 100644 rust/crates/og-dev/src/build/sv_rules.rs diff --git a/rust/Cargo.lock b/rust/Cargo.lock index 4d91a18ca..44bf8f81a 100644 --- a/rust/Cargo.lock +++ b/rust/Cargo.lock @@ -666,6 +666,7 @@ version = "0.1.0" dependencies = [ "flate2", "indexmap", + "og-intervals", "og-io", "og-nf", "png", diff --git a/rust/crates/og-dev/Cargo.toml b/rust/crates/og-dev/Cargo.toml index 20bcd7e4c..91e36986b 100644 --- a/rust/crates/og-dev/Cargo.toml +++ b/rust/crates/og-dev/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "og-dev" -description = "og dev: the developer and CI verbs: parity-compare, contract, schema, report-only" +description = "og dev: the developer and CI verbs: parity-compare, contract, schema, report-only, build" version.workspace = true edition.workspace = true license.workspace = true @@ -10,6 +10,7 @@ publish.workspace = true pycompat.workspace = true og-io.workspace = true og-nf.workspace = true +og-intervals.workspace = true flate2.workspace = true indexmap.workspace = true png = "0.18" diff --git a/rust/crates/og-dev/src/build/callability.rs b/rust/crates/og-dev/src/build/callability.rs new file mode 100644 index 000000000..d7ce4dece --- /dev/null +++ b/rust/crates/og-dev/src/build/callability.rs @@ -0,0 +1,628 @@ +//! `og dev build callability`: `assets/callability_genes.tsv`, the genes +//! where short reads are often not placed with confidence, so a "no +//! finding" there needs a reading of its own (step 16b measures each one). +//! +//! A gene is a GENCODE 50 gene with an HGNC id, on chr1-22, X or Y unless it +//! is ALT-only, and protein-coding unless it is a CMRG gene. It is listed +//! with every reason that applies: +//! +//! - `CMRG`: one of GIAB's 273 challenging medically relevant genes, +//! matched by symbol, or by the protein-coding gene that covers most of +//! its CMRG span when GENCODE has renamed it; +//! - `segdup`: half or more of the gene body in a segmental duplication; +//! - `false_duplication`: the gene body overlaps a false duplication of +//! GRCh38 (either copy); +//! - `collapsed`: the gene body overlaps a collapsed duplication; +//! - `ALT-only`: the gene is on an ALT scaffold and nowhere on chr1-22, X, +//! Y or M; +//! +//! and the step that calls it with a method of its own: `parascopy` for the +//! genes of `assets/paralog_loci.tsv`, `pypgx` for the bam genes of +//! `assets/pgx_genes.tsv` callable without ALT contigs, `cn_loci:` +//! for a gene a gene row of `bin/cn_loci.tsv` names; `.` for none. +//! +//! `--work` holds the seven downloads `scripts/ci/build-callability.sh` +//! names. The GTF is read once, line by line; each interval file is merged +//! once per chromosome, so a gene's coverage is a binary search. + +use super::gtf::Record; +use super::{Args, BuildError, CHROMS, Failure, chrom_index, each_line, each_raw_line, int}; +use og_intervals::sets::Set; +use og_io::Utf8; +use std::collections::{HashMap, HashSet}; +use std::fmt::Write as _; +use std::io::Write; +use std::path::Path; + +/// The reasons, in the order a row lists them. +const REASONS: [&str; 5] = [ + "CMRG", + "segdup", + "false_duplication", + "collapsed", + "ALT-only", +]; +const CMRG: u8 = 1; +const SEGDUP: u8 = 1 << 1; +const FALSE_DUP: u8 = 1 << 2; +const COLLAPSED: u8 = 1 << 3; +const ALT_ONLY: u8 = 1 << 4; + +/// A gene of GENCODE on chr1-22, X or Y: the first record of its HGNC id. +struct Gene { + symbol: String, + hgnc: String, + chrom: usize, + start: i64, + end: i64, + coding: bool, +} + +/// What the GTF gives. +#[derive(Default)] +struct Genes { + /// In the GTF's order. + primary: Vec, + by_hgnc: HashMap, + /// Every gene symbol on chr1-22, X, Y or M. + on_primary: HashSet, + /// HGNC id -> (symbol, UCSC ALT scaffold name) of the first + /// protein-coding ALT copy. + alt: HashMap, +} + +/// The ALT scaffolds of NCBI's assembly report: GenBank accession -> UCSC +/// name (GENCODE names them by accession). +fn read_alt(path: &Path) -> Result, BuildError> { + let mut alt = HashMap::new(); + each_line(path, Utf8::Strict, |line| { + let c: Vec<&str> = line.split('\t').collect(); + if !line.starts_with('#') && c.len() >= 10 && c[1] == "alt-scaffold" { + alt.insert(c[4].to_owned(), c[9].to_owned()); + } + Ok(()) + })?; + Ok(alt) +} + +fn read_gtf(path: &Path, alt: &HashMap) -> Result { + let mut g = Genes::default(); + each_raw_line(path, |line| { + let Some(r) = Record::parse(line)? else { + return Ok(()); + }; + if r.feature != b"gene" { + return Ok(()); + } + let chrom = chrom_index(r.seq); + if (chrom.is_some() || r.seq == b"chrM") + && let Some(n) = r.attr("gene_name")? + { + g.on_primary.insert(n.to_owned()); + } + let Some(h) = r.attr("hgnc_id")?.filter(|h| !h.is_empty()) else { + return Ok(()); + }; + let need = |k: &str| { + r.attr(k)?.ok_or_else(|| { + BuildError::Input(format!( + "GTF: a gene with an HGNC id and no {k}: {:?}", + String::from_utf8_lossy(line) + )) + }) + }; + if need("gene_id")?.ends_with("_PAR_Y") { + return Ok(()); + } + let coding = r.attr("gene_type")? == Some("protein_coding"); + match chrom { + Some(chrom) if !g.by_hgnc.contains_key(h) => { + let (start, end) = r.span()?; + g.by_hgnc.insert(h.to_owned(), g.primary.len()); + g.primary.push(Gene { + symbol: need("gene_name")?.to_owned(), + hgnc: h.to_owned(), + chrom, + start, + end, + coding, + }); + } + _ => { + if let Some(name) = std::str::from_utf8(r.seq).ok().and_then(|s| alt.get(s)) + && coding + && !g.alt.contains_key(h) + { + g.alt + .insert(h.to_owned(), (need("gene_name")?.to_owned(), name.clone())); + } + } + } + Ok(()) + })?; + Ok(g) +} + +/// A BED file as a [`Set`] per sequence. +fn read_sets(paths: &[&Path]) -> Result, BuildError> { + let mut rows: HashMap> = HashMap::new(); + for path in paths { + each_line(path, Utf8::Strict, |line| { + let c: Vec<&str> = line.split('\t').collect(); + if c.len() >= 3 + && !(line.starts_with('#') + || line.starts_with("track") + || line.starts_with("browser")) + { + let s = u64::try_from(int(c[1], "BED start")?); + let e = u64::try_from(int(c[2], "BED end")?); + let (Ok(s), Ok(e)) = (s, e) else { + return Err(BuildError::Input(format!( + "{}: a negative position: {line:?}", + path.display() + ))); + }; + rows.entry(c[0].to_owned()).or_default().push((s, e)); + } + Ok(()) + })?; + } + Ok(rows.into_iter().map(|(k, v)| (k, Set::new(v))).collect()) +} + +/// Bases of `[s, e)` on `chrom` inside `sets`. +fn covered(sets: &HashMap, chrom: usize, s: i64, e: i64) -> u64 { + let (s, e) = (u64::try_from(s).unwrap_or(0), u64::try_from(e).unwrap_or(0)); + sets.get(CHROMS[chrom]).map_or(0, |set| set.covered(s, e)) +} + +/// The rows of a tab-separated table after its comment lines, as +/// `csv.DictReader` with `QUOTE_NONE` reads them: the first line is the +/// header; each row is handed to `f` with a cell getter. +fn each_row( + path: &Path, + mut f: impl FnMut(&dyn Fn(&str) -> Result) -> Result<(), BuildError>, +) -> Result<(), BuildError> { + let mut header: Option> = None; + each_line(path, Utf8::Strict, |line| { + if pycompat::text::strip(line, None).is_empty() || line.starts_with('#') { + return Ok(()); + } + let cells: Vec<&str> = line.split('\t').collect(); + let Some(h) = &header else { + header = Some(cells.iter().map(|c| (*c).to_owned()).collect()); + return Ok(()); + }; + let get = |k: &str| { + h.iter() + .rposition(|n| n == k) + .and_then(|i| cells.get(i)) + .map(|c| (*c).to_owned()) + .ok_or_else(|| { + BuildError::Input(format!("{}: no {k} cell in {line:?}", path.display())) + }) + }; + f(&get) + }) +} + +/// Who calls a gene with a method of its own. +struct Owners { + parascopy: HashSet, + pypgx: HashSet, + /// A name a `bin/cn_loci.tsv` gene row reports -> its locus (the first + /// row that names it). + cn: HashMap, +} + +impl Owners { + fn read(repo: &Path) -> Result { + let mut parascopy = HashSet::new(); + each_row(&repo.join("assets/paralog_loci.tsv"), |get| { + parascopy.extend(get("genes")?.split(',').map(str::to_owned)); + Ok(()) + })?; + // A gene step 32 skips on the default no-ALT reference + // (callable_without_alt no) has no caller there, and step 32 + // refuses a BAM with ALT contigs, so it is not pypgx's. + let mut pypgx = HashSet::new(); + each_row(&repo.join("assets/pgx_genes.tsv"), |get| { + if get("mode")? == "bam" && get("callable_without_alt")? != "no" { + pypgx.insert(get("gene")?); + } + Ok(()) + })?; + // A gene row reports the genes it names (its locus, as RHD, or LPA + // for LPA_KIV2, or a target, as HBA1 and HBA2), not every gene whose + // span holds one of its targets. + let mut cn: HashMap = HashMap::new(); + each_row(&repo.join("bin/cn_loci.tsv"), |get| { + if get("kind")? != "gene" { + return Ok(()); + } + let locus = get("locus")?; + let targets = get("targets")?; + let mut names = vec![locus.as_str(), first(&locus, '_')]; + for part in targets.split(';') { + let n = first(part, ':'); + names.extend([n, first(n, '_')]); + } + for n in names { + cn.entry(n.to_owned()).or_insert_with(|| locus.clone()); + } + Ok(()) + })?; + Ok(Owners { + parascopy, + pypgx, + cn, + }) + } + + fn of(&self, symbol: &str) -> String { + if self.parascopy.contains(symbol) { + "parascopy".to_owned() + } else if self.pypgx.contains(symbol) { + "pypgx".to_owned() + } else if let Some(locus) = self.cn.get(symbol) { + format!("cn_loci:{locus}") + } else { + ".".to_owned() + } + } +} + +/// The text before the first `sep`. +fn first(s: &str, sep: char) -> &str { + s.split(sep).next().unwrap_or(s) +} + +/// The reasons of each gene, by HGNC id. +fn reasons(work: &Path, genes: &Genes) -> Result, BuildError> { + let mut why: HashMap = HashMap::new(); + let mut by_symbol: HashMap<&str, usize> = HashMap::new(); + for (i, g) in genes.primary.iter().enumerate() { + by_symbol.insert(&g.symbol, i); + } + let cmrg = work.join("cmrg.bed"); + each_line(&cmrg, Utf8::Strict, |line| { + let c: Vec<&str> = line.split('\t').collect(); + if c.len() < 4 { + return Err(BuildError::Input(format!( + "cmrg.bed: a row with under 4 columns: {line:?}" + ))); + } + let (chrom, sym) = (c[0], c[3]); + let (s, e) = (int(c[1], "cmrg.bed start")?, int(c[2], "cmrg.bed end")?); + let i = match by_symbol.get(sym) { + Some(&i) => i, + None => { + // Renamed since CMRG v1.00: the protein-coding gene that + // covers most of its span (the larger HGNC id on a tie). + let best = genes + .primary + .iter() + .enumerate() + .filter(|(_, g)| g.coding && CHROMS[g.chrom] == chrom) + .map(|(i, g)| (e.min(g.end) - s.max(g.start), &g.hgnc, i)) + .max_by(|a, b| (a.0, a.1).cmp(&(b.0, b.1))); + match best { + Some((overlap, _, i)) if overlap > 0 => i, + _ => { + return Err(BuildError::Data(format!( + "CMRG gene {sym} ({chrom}:{s}-{e}) matches no GENCODE 50 gene" + ))); + } + } + } + }; + *why.entry(genes.primary[i].hgnc.clone()).or_default() |= CMRG; + Ok(()) + })?; + let segdups = read_sets(&[&work.join("segdups.bed.gz")])?; + let collapsed = read_sets(&[&work.join("collapsed.bed.gz")])?; + let falsedup = read_sets(&[ + &work.join("falsedup_correct.bed.gz"), + &work.join("falsedup_incorrect.bed.gz"), + ])?; + for g in genes.primary.iter().filter(|g| g.coding) { + let mut bits = 0; + // Half or more: covered >= 0.5 * length, in whole numbers. + if 2 * covered(&segdups, g.chrom, g.start, g.end) + >= u64::try_from(g.end - g.start).unwrap_or(0) + { + bits |= SEGDUP; + } + if covered(&falsedup, g.chrom, g.start, g.end) > 0 { + bits |= FALSE_DUP; + } + if covered(&collapsed, g.chrom, g.start, g.end) > 0 { + bits |= COLLAPSED; + } + if bits != 0 { + *why.entry(g.hgnc.clone()).or_default() |= bits; + } + } + for (h, (sym, _)) in &genes.alt { + if !genes.by_hgnc.contains_key(h) && !genes.on_primary.contains(sym) { + *why.entry(h.clone()).or_default() |= ALT_ONLY; + } + } + Ok(why) +} + +/// A row of the catalog: its sort key (chrom order, start, symbol, HGNC +/// id; ALT-only genes last), its cells and its reasons. +type CatalogRow<'a> = ((usize, i64, &'a str, &'a str), [String; 7], u8); + +/// The catalog's text and the summary line, from the downloads in `work` +/// and the owner tables of `repo`. +pub fn build(work: &Path, repo: &Path) -> Result<(String, String), BuildError> { + let alt = read_alt(&work.join("report.txt"))?; + let genes = read_gtf(&work.join("gencode.gtf.gz"), &alt)?; + let why = reasons(work, &genes)?; + let owners = Owners::read(repo)?; + let mut rows: Vec = Vec::with_capacity(why.len()); + for (h, &bits) in &why { + let (key, chrom, start, end, sym) = match genes.by_hgnc.get(h) { + Some(&i) => { + let g = &genes.primary[i]; + let key = (g.chrom, g.start, g.symbol.as_str(), h.as_str()); + ( + key, + CHROMS[g.chrom].to_owned(), + g.start.to_string(), + g.end.to_string(), + &g.symbol, + ) + } + None => { + let (sym, name) = &genes.alt[h]; + ( + (CHROMS.len(), 0, sym.as_str(), h.as_str()), + name.clone(), + ".".to_owned(), + ".".to_owned(), + sym, + ) + } + }; + let listed: Vec<&str> = REASONS + .iter() + .enumerate() + .filter(|(i, _)| bits & (1 << i) != 0) + .map(|(_, r)| *r) + .collect(); + let cells = [ + sym.clone(), + h.clone(), + chrom, + start, + end, + listed.join(","), + owners.of(sym), + ]; + rows.push((key, cells, bits)); + } + rows.sort_by(|a, b| a.0.cmp(&b.0)); + let mut seen: HashMap<&str, usize> = HashMap::new(); + for r in &rows { + *seen.entry(r.0.2).or_default() += 1; + } + let mut dup: Vec<&str> = seen + .into_iter() + .filter(|&(_, n)| n > 1) + .map(|(s, _)| s) + .collect(); + if !dup.is_empty() { + dup.sort_unstable(); + return Err(BuildError::Data(format!( + "gene symbols listed twice: {}", + dup.join(", ") + ))); + } + let counts: Vec = REASONS + .iter() + .enumerate() + .map(|(i, r)| { + format!( + "{r} {}", + rows.iter().filter(|x| x.2 & (1 << i) != 0).count() + ) + }) + .collect(); + let owned = rows.iter().filter(|r| r.1[6] != ".").count(); + let summary = format!( + "Genes: {} ({}); with an owner: {owned}", + rows.len(), + counts.join(", ") + ); + let mut text = String::from( + "# assets/callability_genes.tsv: the genes where short reads are often not placed with confidence, so a\n\ + # 'no finding' there needs a reading of its own. Step 16b measures each gene on chr1-22, X and Y; the\n\ + # report lists the ones it could not read. Written by scripts/ci/build-callability.sh, which says how the\n\ + # rows and the owners are chosen; rebuild it there, do not edit it by hand. start is 0-based and end\n\ + # 1-based, as in a BED file; an ALT-only gene has its ALT scaffold as chrom and '.' as start and end.\n", + ); + let _ = writeln!(text, "# {summary}."); + text.push_str("gene\thgnc_id\tchrom\tstart\tend\treasons\towner\n"); + for r in &rows { + text.push_str(&r.1.join("\t")); + text.push('\n'); + } + Ok((text, summary)) +} + +pub(super) fn cli(args: &[String], err: &mut dyn Write) -> Result { + let a = Args::parse(args, &["--work", "--repo", "--out"], &[])?; + let (work, repo, out) = (a.need("--work")?, a.need("--repo")?, a.need("--out")?); + let (text, summary) = build(Path::new(work), Path::new(repo))?; + super::write_atomic(Path::new(out), text.as_bytes())?; + let _ = writeln!(err, "{summary}"); + Ok(0) +} + +#[cfg(test)] +mod tests { + use super::*; + use std::fs; + + /// A work folder and a repo folder with tiny inputs: A is a CMRG gene by + /// symbol, OLD a CMRG gene GENCODE renamed to B, C sits in a segmental + /// duplication, D over a false and a collapsed duplication, E is only on + /// an ALT scaffold, F is on an ALT scaffold and chr1, P is a CMRG gene whose + /// _PAR_Y copy comes first and must not move it to chrY, N is not protein-coding and has no reason. + fn fixture(dir: &Path) -> (std::path::PathBuf, std::path::PathBuf) { + let (work, repo) = (dir.join("work"), dir.join("repo")); + fs::create_dir_all(&work).unwrap(); + fs::create_dir_all(repo.join("assets")).unwrap(); + fs::create_dir_all(repo.join("bin")).unwrap(); + let gene = |seq: &str, s: u32, e: u32, id: &str, ty: &str, name: &str, h: &str| { + let hg = if h.is_empty() { + String::new() + } else { + format!(" hgnc_id \"{h}\";") + }; + format!( + "{seq}\tX\tgene\t{s}\t{e}\t.\t+\t.\tgene_id \"{id}\"; gene_type \"{ty}\"; gene_name \"{name}\";{hg}\n" + ) + }; + let gtf = [ + "##format: gtf\n".to_owned(), + gene("chr1", 101, 200, "G1", "protein_coding", "A", "HGNC:1"), + "chr1\tX\texon\t101\t150\t.\t+\t.\tgene_id \"G1\";\n".to_owned(), + gene("chr1", 1001, 2000, "G2", "protein_coding", "B", "HGNC:2"), + gene("chr2", 101, 1100, "G3", "protein_coding", "C", "HGNC:3"), + gene("chr2", 5001, 6000, "G4", "protein_coding", "D", "HGNC:4"), + gene("chr2", 9001, 9100, "G5", "lncRNA", "N", "HGNC:5"), + gene( + "chrY", + 100_001, + 100_500, + "G6_PAR_Y", + "protein_coding", + "P", + "HGNC:6", + ), + gene( + "chrX", + 100_001, + 100_500, + "G6", + "protein_coding", + "P", + "HGNC:6", + ), + gene("chr1", 3001, 3100, "G7", "protein_coding", "F", "HGNC:7"), + gene("KI1.1", 1, 100, "G8", "protein_coding", "E", "HGNC:8"), + gene("KI1.1", 1, 100, "G9", "protein_coding", "F", "HGNC:9"), + gene("chrM", 1, 50, "G10", "protein_coding", "MT", ""), + ] + .concat(); + let mut gz = flate2::write::GzEncoder::new(Vec::new(), flate2::Compression::fast()); + gz.write_all(gtf.as_bytes()).unwrap(); + fs::write(work.join("gencode.gtf.gz"), gz.finish().unwrap()).unwrap(); + fs::write( + work.join("report.txt"), + "# report\nKI1\talt-scaffold\t1\tc\tKI1.1\t=\tx\ty\tz\tchr1_KI1v1_alt\r\n", + ) + .unwrap(); + fs::write( + work.join("cmrg.bed"), + "chr1\t100\t200\tA\nchr1\t900\t1500\tOLD\nchrX\t100000\t100500\tP\n", + ) + .unwrap(); + fs::write( + work.join("segdups.bed.gz"), + "chr2\t100\t400\nchr2\t350\t700\n", + ) + .unwrap(); + fs::write(work.join("falsedup_correct.bed.gz"), "chr2\t5999\t6001\n").unwrap(); + fs::write(work.join("falsedup_incorrect.bed.gz"), "").unwrap(); + fs::write(work.join("collapsed.bed.gz"), "#c\nchr2\t5000\t5001\n").unwrap(); + fs::write( + repo.join("assets/paralog_loci.tsv"), + "# x\nlocus\tgenes\nl\tC,Z\n", + ) + .unwrap(); + fs::write( + repo.join("assets/pgx_genes.tsv"), + "gene\tmode\tcallable_without_alt\nD\tbam\tyes\nA\tbam\tno\n", + ) + .unwrap(); + fs::write( + repo.join("bin/cn_loci.tsv"), + "locus\tkind\ttargets\nA_x\tgene\tq_t1:1-2\nY\tinterval\tA:1-2\n", + ) + .unwrap(); + (work, repo) + } + + #[test] + fn the_tiny_catalog_has_every_reason_and_owner() { + let dir = std::env::temp_dir().join(format!("og-callability-{}", std::process::id())); + let (work, repo) = fixture(&dir); + let (text, summary) = build(&work, &repo).unwrap(); + fs::remove_dir_all(&dir).unwrap(); + let body: Vec<&str> = text.lines().filter(|l| !l.starts_with('#')).collect(); + assert_eq!( + body, + [ + "gene\thgnc_id\tchrom\tstart\tend\treasons\towner", + "A\tHGNC:1\tchr1\t100\t200\tCMRG\tcn_loci:A_x", + "B\tHGNC:2\tchr1\t1000\t2000\tCMRG\t.", + "C\tHGNC:3\tchr2\t100\t1100\tsegdup\tparascopy", + "D\tHGNC:4\tchr2\t5000\t6000\tfalse_duplication,collapsed\tpypgx", + "P\tHGNC:6\tchrX\t100000\t100500\tCMRG\t.", + "E\tHGNC:8\tchr1_KI1v1_alt\t.\t.\tALT-only\t.", + ] + ); + assert_eq!( + summary, + "Genes: 6 (CMRG 3, segdup 1, false_duplication 1, collapsed 1, ALT-only 1); with an owner: 3" + ); + assert!(text.contains(&format!("\n# {summary}.\n"))); + } + + #[test] + fn a_cmrg_gene_with_no_gencode_gene_stops_the_build() { + let dir = std::env::temp_dir().join(format!("og-callability-cmrg-{}", std::process::id())); + let (work, repo) = fixture(&dir); + fs::write(work.join("cmrg.bed"), "chr3\t1\t2\tGONE\n").unwrap(); + let e = build(&work, &repo).err().unwrap(); + fs::remove_dir_all(&dir).unwrap(); + assert_eq!( + e.to_string(), + "CMRG gene GONE (chr3:1-2) matches no GENCODE 50 gene" + ); + assert!(matches!(e, BuildError::Data(_))); + } + + #[test] + fn a_symbol_listed_twice_stops_the_build() { + let dir = std::env::temp_dir().join(format!("og-callability-dup-{}", std::process::id())); + let (work, repo) = fixture(&dir); + // Two genes named A, both CMRG (the second by its span). + let mut gz = flate2::write::GzEncoder::new(Vec::new(), flate2::Compression::fast()); + gz.write_all(b"chr1\tX\tgene\t101\t200\t.\t+\t.\tgene_id \"G1\"; gene_type \"protein_coding\"; gene_name \"A\"; hgnc_id \"HGNC:1\";\nchr1\tX\tgene\t901\t1500\t.\t+\t.\tgene_id \"G2\"; gene_type \"protein_coding\"; gene_name \"A\"; hgnc_id \"HGNC:2\";\n").unwrap(); + fs::write(work.join("gencode.gtf.gz"), gz.finish().unwrap()).unwrap(); + fs::write( + work.join("cmrg.bed"), + "chr1\t100\t200\tX1\nchr1\t900\t1500\tX2\n", + ) + .unwrap(); + let e = build(&work, &repo).err().unwrap(); + fs::remove_dir_all(&dir).unwrap(); + assert_eq!(e.to_string(), "gene symbols listed twice: A"); + } + + #[test] + fn coverage_is_union_coverage() { + let sets = HashMap::from([( + "chr2".to_owned(), + Set::new([(100, 400), (350, 700), (900, 950)]), + )]); + assert_eq!(covered(&sets, 1, 100, 1100), 650); + assert_eq!(covered(&sets, 0, 100, 1100), 0); + } +} diff --git a/rust/crates/og-dev/src/build/cn_intervals.rs b/rust/crates/og-dev/src/build/cn_intervals.rs new file mode 100644 index 000000000..8624daa9b --- /dev/null +++ b/rust/crates/og-dev/src/build/cn_intervals.rs @@ -0,0 +1,1162 @@ +//! `og dev build cn-intervals`: the interval rows of `bin/cn_loci.tsv`, +//! the copy-number checks step 35 runs with one `samtools bedcov` pass +//! (`og cn intervals`). +//! +//! The rows are: +//! +//! - every region row of ClinGen Dosage Sensitivity's recurrent CNV BED, +//! version 2.1, GRCh38 (`CLINGEN_RECURRENT_CNV_URL` in `versions.env`, +//! CC0). Its breakpoint and LCR rows (itemRgb `0,0,0`) are the segmental +//! duplications that bound the regions: they are cut out of the targets +//! and kept out of the controls. Its chrX rows are spelled `chrx` and are +//! read as chrX; any other name that is not a GRCh38 chromosome stops the +//! build; +//! - the rows the BED lacks (`HAND_ROWS`): two ClinGen curations (1p36, +//! Xq28 MECP2), five carrier deletions at their ClinVar or HbVar GRCh38 +//! spans, and the Y chromosome's AZF regions with two single-copy Y +//! controls. +//! +//! Targets are cut to pieces of at most 500 kb and lose every stretch inside +//! a breakpoint or LCR row, a GIAB v3.6 segmental duplication, an assembly +//! gap (GIAB's gaps with 15 kb on each side), a satellite or a region of a +//! model the pinned Parascopy data holds; the AZF targets keep their +//! segmental duplications, since the DAZ copies are one. Controls are two +//! 100 kb windows, one on each side of the cluster of ClinGen rows the +//! region belongs to, the nearest that stay clear of all of that and of +//! every target, within 3 Mb; a side with none (a chromosome end) is left +//! out. +//! +//! Every cut set is merged once per chromosome, so a subtraction or a +//! "clear of the cuts" test is a binary search. + +use super::{Args, BuildError, CHROMS, Failure, chrom_index, each_line, input_error, int}; +use og_io::Utf8; +use std::io::Write; +use std::path::Path; + +/// The line after which the catalog's rows are this builder's. +pub const MARKER: &str = + "# --- interval rows: written by scripts/ci/build-cn-intervals.sh; do not edit by hand ---"; +/// How a marker line starts, whatever builder it names. +const MARKER_START: &[u8] = b"# --- interval rows: "; + +const REGION_RGB: &str = "255,111,0"; +const MASK_RGB: &str = "0,0,0"; +/// GRCh38 chromosome lengths (the no-ALT analysis set's `.fai`), in +/// [`CHROMS`] order. +const LENGTHS: [i64; 24] = [ + 248_956_422, + 242_193_529, + 198_295_559, + 190_214_555, + 181_538_259, + 170_805_979, + 159_345_973, + 145_138_636, + 138_394_717, + 133_797_422, + 135_086_622, + 133_275_309, + 114_364_328, + 107_043_718, + 101_991_189, + 90_338_345, + 83_257_441, + 80_373_285, + 58_617_616, + 64_444_167, + 46_709_983, + 50_818_468, + 156_040_895, + 57_227_415, +]; +const PIECE: i64 = 500_000; +const FLANK: i64 = 100_000; +const GAP: i64 = 50_000; +const SEARCH: i64 = 3_000_000; +const STEP: i64 = 10_000; +const MIN_PIECE: i64 = 500; +const CHR_Y: usize = 23; + +/// A 0-based, half-open stretch. +type Span = (i64, i64); + +/// How a row's target depth is read. +#[derive(Clone, Copy, PartialEq, Eq)] +enum Depth { + /// Depth pooled over the targets. + Pooled, + /// The targets' depths added up (copies of one sequence). + Sum, +} + +impl Depth { + fn name(self) -> &'static str { + match self { + Depth::Pooled => "pooled", + Depth::Sum => "sum", + } + } +} + +/// A row the BED lacks. +struct Hand { + locus: &'static str, + chrom: usize, + /// (label, 0-based start, end); the first label names the row. + targets: &'static [(&'static str, i64, i64)], + /// (name, start, end, label), or `None` for flanks. + controls: Option<&'static [(&'static str, i64, i64, &'static str)]>, + depth: Depth, + /// Whether segmental duplications are cut out of the targets. + cut_segdups: bool, + pmids: &'static str, +} + +const Y_CONTROLS: &[(&str, i64, i64, &str)] = &[ + ( + "ZFY", + 2_935_200, + 2_982_579, + "ZFY, single-copy control on Yp", + ), + ( + "NLGN4Y", + 14_522_572, + 14_845_654, + "NLGN4Y, single-copy control on Yq", + ), +]; + +/// The rows the ClinGen BED lacks. 1p36 (ISCA-37434) and Xq28 (ISCA-46304) +/// from ClinGen_region_curation_list_GRCh38.tsv as read on 2026-11-04; +/// del(GJB6-D13S1830), CFTRdele2,3 and the 621 bp HBB deletion at their +/// ClinVar GRCh38 spans; the 619 bp HBB deletion (HbVar 979, Orkin 1979) +/// and the Sicilian (delta-beta)0 deletion (HbVar 1035, Henthorn 1990) at +/// their HbVar spans, placed on NC_000011.10 by 50 bp on each side of the +/// NG_000007.3 breakpoints on 2026-11-04; the AZF regions by the single-copy +/// genes each deletion removes (AZFc as the four DAZ copies added up), gene +/// spans from Ensembl's REST lookup, checked on 2026-11-04. The alpha-globin +/// --SEA deletion is read by step 35's HBA table, not here. +const HAND_ROWS: &[Hand] = &[ + Hand { + locus: "1p36_terminal", + chrom: 0, + targets: &[( + "1p36 terminal region (includes GABRD), ClinGen ISCA-37434", + 898_702, + 6_229_913, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: ".", + }, + Hand { + locus: "Xq28_MECP2", + chrom: 22, + targets: &[( + "Xq28 region (includes MECP2), ClinGen ISCA-46304: the MECP2 duplication", + 154_008_528, + 154_110_279, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: ".", + }, + Hand { + locus: "GJB6_D13S1830", + chrom: 12, + targets: &[("del(GJB6-D13S1830), ClinVar 5546", 20_223_037, 20_531_806)], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: "11807148", + }, + Hand { + locus: "CFTR_dele2_3", + chrom: 6, + targets: &[( + "CFTRdele2,3 (21 kb), ClinVar 66105", + 117_498_312, + 117_519_392, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: ".", + }, + Hand { + locus: "HBB_exon3_del", + chrom: 10, + targets: &[( + "HBB deletion of exon 3 (621 bp), ClinVar 3769040", + 5_225_254, + 5_225_875, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: ".", + }, + Hand { + locus: "HBB_619bp_del", + chrom: 10, + targets: &[( + "HBB 619 bp deletion (South Asian), HbVar 979", + 5_225_388, + 5_226_007, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: "287080", + }, + Hand { + locus: "HBD_HBB_Sicilian_del", + chrom: 10, + targets: &[( + "Sicilian (delta-beta)0 deletion (13.4 kb, HBD to HBB 3'), HbVar 1035", + 5_219_877, + 5_233_280, + )], + controls: None, + depth: Depth::Pooled, + cut_segdups: true, + pmids: "2307466", + }, + Hand { + locus: "AZFa", + chrom: CHR_Y, + targets: &[ + ("AZFa (USP9Y and DDX3Y)", 12_537_649, 12_860_852), + ("DDX3Y", 12_904_107, 12_920_478), + ], + controls: Some(Y_CONTROLS), + depth: Depth::Pooled, + cut_segdups: false, + pmids: "37674303", + }, + Hand { + locus: "AZFb", + chrom: CHR_Y, + targets: &[ + ("AZFb (KDM5D and EIF1AY)", 19_702_297, 19_745_148), + ("EIF1AY", 20_575_693, 20_593_154), + ], + controls: Some(Y_CONTROLS), + depth: Depth::Pooled, + cut_segdups: false, + pmids: "37674303", + }, + Hand { + locus: "AZFc_DAZ", + chrom: CHR_Y, + targets: &[ + ( + "AZFc (DAZ1 to DAZ4 added up; gr/gr and b2/b3 leave two)", + 23_129_354, + 23_199_010, + ), + ("DAZ2", 23_219_446, 23_291_356), + ("DAZ3", 24_763_068, 24_813_492), + ("DAZ4", 24_833_842, 24_907_040), + ], + controls: Some(Y_CONTROLS), + depth: Depth::Sum, + cut_segdups: false, + pmids: "37674303", + }, +]; + +/// A region row of the ClinGen BED. +#[derive(Debug)] +struct Region { + chrom: usize, + start: i64, + end: i64, + name: String, +} + +/// The stretches cut out of targets and controls, per chromosome. +#[derive(Default, Clone)] +struct Cuts { + /// ClinGen breakpoint and LCR rows. + mask: Vec, + segdup: Vec, + gap: Vec, + satellite: Vec, + parascopy: Vec, +} + +/// The downloads, read. +pub struct Inputs { + regions: Vec, + cuts: Vec, +} + +/// `chrx` and `chry` (the ClinGen BED's spelling) as chrX and chrY; any +/// other name that is not a GRCh38 chromosome is an error. +fn norm_chrom(c: &str) -> Result { + let fixed = match c { + "chrx" => "chrX", + "chry" => "chrY", + _ => c, + }; + chrom_index(fixed).ok_or_else(|| { + BuildError::Input(format!( + "chromosome {} is not a GRCh38 chromosome name", + py_repr(c) + )) + }) +} + +fn py_repr(s: &str) -> String { + pycompat::value::py_repr(&pycompat::json::Value::Str(s.to_owned())) + .unwrap_or_else(|_| format!("'{s}'")) +} + +/// The region rows and breakpoint rows of the ClinGen BED's lines. +fn read_clingen<'a>( + lines: impl IntoIterator, + cuts: &mut [Cuts], +) -> Result, BuildError> { + let mut regions = Vec::new(); + for line in lines { + let strip = pycompat::text::strip; + if strip(line, None).is_empty() + || line.starts_with('#') + || line.starts_with("track") + || line.starts_with("browser") + { + continue; + } + let c: Vec<&str> = line.split('\t').collect(); + if c.len() < 9 { + return Err(BuildError::Input(format!( + "ClinGen BED: a row with {} columns: {}", + c.len(), + py_repr(strip(line, None)) + ))); + } + let chrom = norm_chrom(c[0])?; + let (start, end) = ( + int(c[1], "ClinGen BED start")?, + int(c[2], "ClinGen BED end")?, + ); + let (name, rgb) = (strip(c[3], None), strip(c[8], None)); + if end <= start { + return Err(BuildError::Input(format!( + "ClinGen BED: {name}: empty span" + ))); + } + match rgb { + REGION_RGB => regions.push(Region { + chrom, + start, + end, + name: name.to_owned(), + }), + MASK_RGB => cuts[chrom].mask.push((start, end)), + _ => { + return Err(BuildError::Input(format!( + "ClinGen BED: {name}: itemRgb {} is neither a region nor a breakpoint row", + py_repr(rgb) + ))); + } + } + } + if regions.is_empty() { + return Err(BuildError::Input("ClinGen BED: no region row".to_owned())); + } + Ok(regions) +} + +/// The rows of a BED on chr1-22, X and Y, per chromosome. +fn read_bed( + path: &Path, + put: impl Fn(&mut Cuts) -> &mut Vec, + cuts: &mut [Cuts], +) -> Result<(), BuildError> { + each_line(path, Utf8::Strict, |line| { + if line.starts_with('#') || pycompat::text::strip(line, None).is_empty() { + return Ok(()); + } + let mut c = line.split('\t'); + let (seq, s, e) = (c.next(), c.next(), c.next()); + let (Some(seq), Some(s), Some(e)) = (seq, s, e) else { + return Err(BuildError::Input(format!( + "{}: a row with under 3 columns", + path.display() + ))); + }; + if let Some(i) = chrom_index(seq) { + let span = (int(s, "BED start")?, int(e, "BED end")?); + put(&mut cuts[i]).push(span); + } + Ok(()) + }) +} + +/// The regions of Parascopy's models: `NAME CHROM:START-END ...`, 1-based. +fn read_parascopy(path: &Path, cuts: &mut [Cuts]) -> Result<(), BuildError> { + each_line(path, Utf8::Strict, |line| { + if line.starts_with('#') || pycompat::text::strip(line, None).is_empty() { + return Ok(()); + } + let bad = || { + BuildError::Input(format!( + "{}: not NAME CHROM:START-END DETAIL: {line:?}", + path.display() + )) + }; + let f = pycompat::text::split(line, None, None); + let [_, reg, _] = f[..] else { + return Err(bad()); + }; + let (c, span) = reg.split_once(':').ok_or_else(bad)?; + let (s, e) = span.split_once('-').ok_or_else(bad)?; + if let Some(i) = chrom_index(c) { + cuts[i] + .parascopy + .push((int(s, "Parascopy start")? - 1, int(e, "Parascopy end")?)); + } + Ok(()) + }) +} + +/// Sorted, with overlapping and touching spans joined. +fn merge(mut spans: Vec) -> Vec { + spans.sort_unstable(); + let mut out: Vec = Vec::with_capacity(spans.len()); + for (s, e) in spans { + match out.last_mut() { + Some(last) if s <= last.1 => last.1 = last.1.max(e), + _ => out.push((s, e)), + } + } + out +} + +/// The parts of `[s, e)` outside the merged `cuts`. +fn subtract(s: i64, e: i64, cuts: &[Span]) -> Vec { + let mut out = Vec::new(); + let mut pos = s; + let first = cuts.partition_point(|c| c.1 <= s); + for &(cs, ce) in &cuts[first..] { + if cs >= e { + break; + } + if cs > pos { + out.push((pos, cs)); + } + pos = pos.max(ce); + } + if pos < e { + out.push((pos, e)); + } + out +} + +/// Each span of `MIN_PIECE` bp or more cut into equal pieces of at most +/// `PIECE` bp. +fn pieces(spans: &[Span]) -> Vec { + let mut out = Vec::new(); + for &(s, e) in spans { + let len = e - s; + if len < MIN_PIECE { + continue; + } + // Ceiling divisions of positive numbers. + let n = (len + PIECE - 1) / PIECE; + let step = (len + n - 1) / n; + out.extend((0..n).map(|i| (s + i * step, e.min(s + (i + 1) * step)))); + } + out +} + +/// Whether `[s, e)` touches none of the merged `cuts`. +fn clear(s: i64, e: i64, cuts: &[Span]) -> bool { + let i = cuts.partition_point(|c| c.1 <= s); + cuts.get(i).is_none_or(|c| c.0 >= e) +} + +/// The nearest `FLANK` window `GAP` beyond `edge` on the left (`left`) or +/// right, clear of `cuts`, within `SEARCH`. +fn flank(chrom: usize, edge: i64, left: bool, cuts: &[Span]) -> Option { + for k in 0..=SEARCH / STEP { + let (s, e) = if left { + let e = edge - GAP - k * STEP; + (e - FLANK, e) + } else { + let s = edge + GAP + k * STEP; + (s, s + FLANK) + }; + if s < 0 || e > LENGTHS[chrom] { + return None; + } + if clear(s, e, cuts) { + return Some((s, e)); + } + } + None +} + +/// The row id of a ClinGen region name: its band, then the words of its +/// parentheses that are not an "includes" note. +fn locus_id(name: &str) -> String { + let band = pycompat::text::split(name, None, None) + .first() + .copied() + .unwrap_or_default(); + let b = name.as_bytes(); + let mut groups: Vec<&str> = Vec::new(); + let mut i = 0; + while i < b.len() { + if b[i] == b'(' { + let j = b[i + 1..] + .iter() + .position(|&c| c == b'(' || c == b')') + .map(|k| i + 1 + k); + if let Some(j) = j + && b[j] == b')' + { + let g = &name[i + 1..j]; + if !pycompat::text::strip(g, None).starts_with("includes") { + groups.push(g); + } + i = j + 1; + continue; + } + } + i += 1; + } + let mut words = String::new(); + let mut in_run = false; + for ch in groups.join(" ").chars() { + if ch.is_ascii_alphanumeric() || ch == '.' || ch == '-' { + words.push(ch); + in_run = false; + } else if !in_run { + words.push('_'); + in_run = true; + } + } + let words = words.trim_matches('_'); + if words.is_empty() { + band.to_owned() + } else { + format!("{band}_{words}") + } +} + +/// A name made safe for a target or control cell. +fn label(text: &str) -> String { + let t: String = text + .chars() + .map(|c| { + if matches!(c, ';' | '\t' | '\n') { + ',' + } else { + c + } + }) + .collect(); + pycompat::text::strip(&t, None).to_owned() +} + +/// A control window: (name, start, end, label). +type Control = (String, i64, i64, String); + +/// The catalog rows, as Python wrote them, in their order. +struct Builder<'a> { + inputs: &'a Inputs, + /// Per chromosome, the ClinGen rows that touch or overlap, merged + /// (each row's end taken one base further). + clusters: Vec>, + /// Per chromosome, every cut (with and without segmental duplications) + /// and the cuts a control keeps out of. + cuts_all: Vec>, + cuts_no_segdups: Vec>, + flank_cuts: Vec>, + rows: Vec, + seen: std::collections::HashMap, +} + +impl<'a> Builder<'a> { + fn new(inputs: &'a Inputs) -> Builder<'a> { + let n = CHROMS.len(); + let mut clusters = vec![Vec::new(); n]; + let mut flank_cuts = vec![Vec::new(); n]; + for r in &inputs.regions { + clusters[r.chrom].push((r.start, r.end + 1)); + flank_cuts[r.chrom].push((r.start, r.end)); + } + for h in HAND_ROWS { + flank_cuts[h.chrom].extend(h.targets.iter().map(|&(_, s, e)| (s, e))); + } + let mut cuts_all = Vec::with_capacity(n); + let mut cuts_no_segdups = Vec::with_capacity(n); + for (i, c) in inputs.cuts.iter().enumerate() { + clusters[i].extend(c.mask.iter().map(|&(s, e)| (s, e + 1))); + let all: Vec = [&c.mask, &c.segdup, &c.gap, &c.satellite, &c.parascopy] + .into_iter() + .flatten() + .copied() + .collect(); + flank_cuts[i].extend(&all); + cuts_all.push(merge(all)); + cuts_no_segdups.push(merge( + [&c.mask, &c.gap, &c.parascopy] + .into_iter() + .flatten() + .copied() + .collect(), + )); + } + Builder { + inputs, + clusters: clusters.into_iter().map(merge).collect(), + cuts_all, + cuts_no_segdups, + flank_cuts: flank_cuts.into_iter().map(merge).collect(), + rows: Vec::new(), + seen: std::collections::HashMap::new(), + } + } + + /// The two flank controls of the cluster that holds `[s, e]`. + fn flanks(&self, chrom: usize, s: i64, e: i64) -> Result, BuildError> { + let (cs, ce) = self.clusters[chrom] + .iter() + .copied() + .find(|&(a, b)| a <= s && e <= b) + .unwrap_or((s, e)); + let mut out = Vec::new(); + for (left, edge, n) in [(true, cs, "L"), (false, ce, "R")] { + if let Some((a, b)) = flank(chrom, edge, left, &self.flank_cuts[chrom]) { + let side = if left { "left" } else { "right" }; + out.push((n.to_owned(), a, b, format!("flank {side} of the region"))); + } + } + if out.is_empty() { + return Err(BuildError::Input(format!( + "{}:{s}-{e}: no flank clear of the masks within {} Mb", + CHROMS[chrom], + SEARCH / 1_000_000 + ))); + } + Ok(out) + } + + fn add( + &mut self, + locus: &str, + chrom: usize, + targets: &[(i64, i64, String)], + controls: &[Control], + depth: Depth, + pmids: &str, + ) -> Result<(), BuildError> { + if targets.is_empty() { + return Err(BuildError::Input(format!( + "{locus}: no target stretch is left outside the masks" + ))); + } + let n = self.seen.entry(locus.to_owned()).or_insert(0); + *n += 1; + let locus = if *n > 1 { + format!("{locus}_{n}") + } else { + locus.to_owned() + }; + let tg: Vec = targets + .iter() + .enumerate() + .map(|(i, (s, e, lab))| { + let tail = if lab.is_empty() { + String::new() + } else { + format!(":{lab}") + }; + format!("{locus}_t{}:{s}-{e}{tail}", i + 1) + }) + .collect(); + let ct: Vec = controls + .iter() + .map(|(n, s, e, lab)| format!("{locus}_{n}:{s}-{e}:{lab}")) + .collect(); + let y = chrom == CHR_Y; + let cells = [ + locus.as_str(), + "interval", + CHROMS[chrom], + &tg.join(";"), + &ct.join(";"), + ".", + depth.name(), + if y { "pooled" } else { "flanks" }, + "2", + "0", + ".", + ".", + ".", + if y { "4" } else { "8" }, + if y { "." } else { "1.5" }, + "round:0.3", + ".", + ".", + "bedcov", + "cn-depth", + "cnvpytor:cnvs", + pmids, + ]; + self.rows.push(cells.join("\t")); + Ok(()) + } + + /// Targets with the first piece labelled. + fn labelled(segs: &[Span], name: &str) -> Vec<(i64, i64, String)> { + segs.iter() + .enumerate() + .map(|(i, &(a, b))| (a, b, if i == 0 { label(name) } else { String::new() })) + .collect() + } + + fn build(mut self) -> Result, BuildError> { + let mut regions: Vec<&Region> = self.inputs.regions.iter().collect(); + regions.sort_by_key(|r| (r.chrom, r.start, r.end)); + for r in regions { + let segs = pieces(&subtract(r.start, r.end, &self.cuts_all[r.chrom])); + let targets = Self::labelled(&segs, &r.name); + let controls = self.flanks(r.chrom, r.start, r.end)?; + self.add( + &locus_id(&r.name), + r.chrom, + &targets, + &controls, + Depth::Pooled, + ".", + )?; + } + for h in HAND_ROWS { + let cuts = if h.cut_segdups { + &self.cuts_all[h.chrom] + } else { + &self.cuts_no_segdups[h.chrom] + }; + let mut targets = Vec::new(); + for &(name, s, e) in h.targets { + let segs = match h.depth { + Depth::Pooled => pieces(&subtract(s, e, cuts)), + Depth::Sum => vec![(s, e)], + }; + if h.depth == Depth::Sum + && !clear(s, e, &merge(self.inputs.cuts[h.chrom].parascopy.clone())) + { + return Err(BuildError::Input(format!( + "{}: {name} overlaps a Parascopy model", + h.locus + ))); + } + targets.extend(Self::labelled(&segs, name)); + } + let controls = match h.controls { + None => self.flanks(h.chrom, h.targets[0].1, h.targets[h.targets.len() - 1].2)?, + Some(c) => c + .iter() + .enumerate() + .map(|(i, &(_, s, e, lab))| (format!("c{}", i + 1), s, e, lab.to_owned())) + .collect(), + }; + self.add(h.locus, h.chrom, &targets, &controls, h.depth, h.pmids)?; + } + Ok(self.rows) + } +} + +impl Inputs { + /// The downloads: the ClinGen BED's lines, the three GIAB BEDs and the + /// Parascopy model list. + pub fn read( + clingen: &Path, + segdups: &Path, + gaps: &Path, + satellites: &Path, + parascopy: &Path, + ) -> Result { + let mut cuts = vec![Cuts::default(); CHROMS.len()]; + let lines: Vec = super::lines(clingen, Utf8::Replace)? + .collect::>() + .map_err(|e| input_error(clingen, &e))?; + let regions = read_clingen(lines.iter().map(String::as_str), &mut cuts)?; + read_bed(segdups, |c| &mut c.segdup, &mut cuts)?; + read_bed(gaps, |c| &mut c.gap, &mut cuts)?; + read_bed(satellites, |c| &mut c.satellite, &mut cuts)?; + read_parascopy(parascopy, &mut cuts)?; + Ok(Inputs { regions, cuts }) + } + + /// The interval rows, without line endings, marker excluded. + pub fn rows(&self) -> Result, BuildError> { + Builder::new(self).build() + } +} + +/// The catalog with its rows from the marker line on replaced: the lines +/// before a marker line are kept (all of them when there is none). +fn splice(catalog: &[u8], rows: &[String]) -> Vec { + let mut out = Vec::with_capacity(catalog.len()); + for line in catalog.split_inclusive(|&c| c == b'\n') { + if line.starts_with(MARKER_START) && line.ends_with(b"\n") { + break; + } + out.extend_from_slice(line); + } + for line in std::iter::once(MARKER).chain(rows.iter().map(String::as_str)) { + out.extend_from_slice(line.as_bytes()); + out.push(b'\n'); + } + out +} + +/// How two catalogs differ: (lines only in `old`, lines only in `new`), +/// each line counted as often as it appears. +fn line_diff(old: &[u8], new: &[u8]) -> (usize, usize) { + let mut count: std::collections::HashMap<&[u8], i64> = std::collections::HashMap::new(); + for l in old.split_inclusive(|&c| c == b'\n') { + *count.entry(l).or_default() += 1; + } + for l in new.split_inclusive(|&c| c == b'\n') { + *count.entry(l).or_default() -= 1; + } + let only_old = count + .values() + .filter(|&&v| v > 0) + .map(|&v| v as usize) + .sum(); + let only_new = count + .values() + .filter(|&&v| v < 0) + .map(|&v| (-v) as usize) + .sum(); + (only_old, only_new) +} + +pub(super) fn cli( + args: &[String], + out: &mut dyn Write, + err: &mut dyn Write, +) -> Result { + let a = Args::parse( + args, + &[ + "--clingen", + "--segdups", + "--gaps", + "--satellites", + "--parascopy", + "--catalog", + ], + &["--check", "--stdout"], + )?; + let catalog = a.get("--catalog"); + if catalog.is_some() == a.has("--stdout") { + return Err(Failure::Usage( + "give one of --catalog FILE and --stdout".to_owned(), + )); + } + if a.has("--check") && catalog.is_none() { + return Err(Failure::Usage("--check needs --catalog".to_owned())); + } + let p = |f: &str| a.need(f).map(Path::new); + let inputs = Inputs::read( + p("--clingen")?, + p("--segdups")?, + p("--gaps")?, + p("--satellites")?, + p("--parascopy")?, + )?; + let rows = inputs.rows()?; + let Some(catalog) = catalog else { + let mut text = String::new(); + for line in std::iter::once(MARKER).chain(rows.iter().map(String::as_str)) { + text.push_str(line); + text.push('\n'); + } + let _ = out.write_all(text.as_bytes()); + return Ok(0); + }; + let path = Path::new(catalog); + let old = std::fs::read(path).map_err(|e| input_error(path, &e))?; + let new = splice(&old, &rows); + if a.has("--check") { + if new != old { + let (only_old, only_new) = line_diff(&old, &new); + let how = if only_old == 0 && only_new == 0 { + "the same lines in another order".to_owned() + } else { + format!("{only_old} line(s) only in the file, {only_new} only in the rebuild") + }; + let _ = writeln!( + err, + "{catalog} is not what the pinned downloads give ({how}); rerun scripts/ci/build-cn-intervals.sh" + ); + return Ok(1); + } + let _ = writeln!( + out, + "{catalog}: the {} interval rows are what the pinned downloads give", + rows.len() + ); + return Ok(0); + } + super::write_atomic(path, &new)?; + let _ = writeln!(out, "{catalog}: {} interval rows written", rows.len()); + Ok(0) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A small ClinGen BED: the 22q11.2 A-D region between its breakpoint + /// rows, and the Xp22.31 region with its breakpoint, spelled `chrx`. + const SMALL_BED: &str = "track name=\"test\"\n#chrom\tstart\tend\tname\tscore\tstrand\tthickStart\tthickEnd\titemRgb\nchr22\t18159724\t18924717\tBP-A\t0\t.\t18159724\t18924717\t0,0,0\nchr22\t18924718\t21111383\t22q11.2 recurrent (DGS) region (proximal, A-D) (includes TBX1)\t0\t.\t18924718\t21111383\t255,111,0\nchr22\t21111384\t21562827\tBP-D\t0\t.\t21111384\t21562827\t0,0,0\nchrx\t6527561\t6537770\tXp22.31 distal breakpoint\t0\t.\t6527561\t6537770\t0,0,0\nchrx\t6537771\t8156913\tXp22.31 recurrent region (includes STS)\t0\t.\t6537771\t8156913\t255,111,0\n"; + + fn small(bed: &str) -> Result, BuildError> { + let mut cuts = vec![Cuts::default(); CHROMS.len()]; + let regions = read_clingen(bed.lines(), &mut cuts)?; + Inputs { regions, cuts }.rows() + } + + fn cell(row: &str, i: usize) -> &str { + row.split('\t').nth(i).unwrap() + } + + /// (start, end) of each target of a row. + fn spans(row: &str) -> Vec { + cell(row, 3) + .split(';') + .map(|t| { + let (s, e) = t.split(':').nth(1).unwrap().split_once('-').unwrap(); + (s.parse().unwrap(), e.parse().unwrap()) + }) + .collect() + } + + #[test] + fn the_small_bed_gives_its_two_regions_then_the_hand_rows() { + let rows = small(SMALL_BED).unwrap(); + let ids: Vec<(&str, &str)> = rows.iter().map(|r| (cell(r, 0), cell(r, 2))).collect(); + assert_eq!( + ids[..2], + [("22q11.2_DGS_proximal_A-D", "chr22"), ("Xp22.31", "chrX")] + ); + let hand: Vec<&str> = ids[2..].iter().map(|x| x.0).collect(); + let want: Vec<&str> = HAND_ROWS.iter().map(|h| h.locus).collect(); + assert_eq!(hand, want, "the hand rows follow, in their order"); + assert!(rows.iter().all(|r| r.split('\t').count() == 22)); + assert!(rows.iter().all(|r| cell(r, 20) == "cnvpytor:cnvs")); + } + + #[test] + fn an_unknown_contig_stops_the_build() { + let e = small(&SMALL_BED.replace("chrx", "chrQ")).unwrap_err(); + assert!( + e.to_string().contains("'chrQ' is not a GRCh38 chromosome"), + "{e}" + ); + assert_eq!(e.code(), 2, "exit 2, as the Python builder did"); + } + + #[test] + fn no_target_piece_overlaps_a_breakpoint_row() { + let rows = small(SMALL_BED).unwrap(); + for (s, e) in spans(&rows[0]) { + assert!(!(s < 18_924_717 && 18_159_724 < e), "{s}-{e} on BP-A"); + assert!(!(s < 21_562_827 && 21_111_384 < e), "{s}-{e} on BP-D"); + } + // 2.19 Mb in pieces of at most 500 kb, the first labelled. + let t = cell(&rows[0], 3); + assert_eq!(t.split(';').count(), 5); + assert!(t.starts_with("22q11.2_DGS_proximal_A-D_t1:18924718-")); + assert!( + t.split(';') + .next() + .unwrap() + .ends_with(":22q11.2 recurrent (DGS) region (proximal, A-D) (includes TBX1)") + ); + } + + #[test] + fn controls_flank_the_cluster_and_never_sit_on_a_target() { + let rows = small(SMALL_BED).unwrap(); + // 22q11.2: the cluster runs from BP-A to BP-D, so the controls lie + // beyond those, 50 kb out. + assert_eq!( + cell(&rows[0], 4), + "22q11.2_DGS_proximal_A-D_L:18009724-18109724:flank left of the region;22q11.2_DGS_proximal_A-D_R:21612828-21712828:flank right of the region" + ); + for r in &rows { + if cell(r, 2) == "chrY" { + continue; + } + for c in cell(r, 4).split(';') { + let (s, e) = c.split(':').nth(1).unwrap().split_once('-').unwrap(); + let (cs, ce): (i64, i64) = (s.parse().unwrap(), e.parse().unwrap()); + for o in &rows { + if cell(o, 2) == cell(r, 2) { + assert!( + spans(o).iter().all(|&(ts, te)| !(ts < ce && cs < te)), + "{} control on {}", + cell(r, 0), + cell(o, 0) + ); + } + } + } + } + } + + #[test] + fn hand_rows_keep_their_spans_and_columns() { + let rows = small(SMALL_BED).unwrap(); + let get = |l: &str| rows.iter().find(|r| cell(r, 0) == l).unwrap().clone(); + assert_eq!(spans(&get("HBB_619bp_del")), [(5_225_388, 5_226_007)]); + assert_eq!( + spans(&get("HBD_HBB_Sicilian_del")), + [(5_219_877, 5_233_280)] + ); + let az = get("AZFc_DAZ"); + assert_eq!( + (cell(&az, 6), cell(&az, 7), cell(&az, 13), cell(&az, 14)), + ("sum", "pooled", "4", ".") + ); + assert_eq!(spans(&az).len(), 4, "the DAZ copies are never cut"); + assert!(cell(&az, 4).starts_with("AZFc_DAZ_c1:2935200-2982579:ZFY")); + let x = get("Xq28_MECP2"); + assert_eq!( + (cell(&x, 7), cell(&x, 13), cell(&x, 14)), + ("flanks", "8", "1.5") + ); + } + + #[test] + fn a_parascopy_model_is_cut_out_of_a_target() { + let mut cuts = vec![Cuts::default(); CHROMS.len()]; + let regions = read_clingen(SMALL_BED.lines(), &mut cuts).unwrap(); + cuts[21].parascopy.push((19_000_000, 19_100_000)); + let rows = Inputs { regions, cuts }.rows().unwrap(); + assert!( + spans(&rows[0]) + .iter() + .all(|&(s, e)| e <= 19_000_000 || s >= 19_100_000) + ); + assert!(spans(&rows[0]).iter().any(|&(_, e)| e == 19_000_000)); + } + + #[test] + fn a_parascopy_model_on_a_daz_copy_stops_the_build() { + let mut cuts = vec![Cuts::default(); CHROMS.len()]; + let regions = read_clingen(SMALL_BED.lines(), &mut cuts).unwrap(); + cuts[CHR_Y].parascopy.push((24_800_000, 24_800_100)); + let e = Inputs { regions, cuts }.rows().unwrap_err(); + assert!( + e.to_string() + .contains("AZFc_DAZ: DAZ3 overlaps a Parascopy model"), + "{e}" + ); + } + + #[test] + fn bad_rows_are_refused() { + let mut c = vec![Cuts::default(); CHROMS.len()]; + assert!( + read_clingen(["chr1\t1\t2"], &mut c) + .unwrap_err() + .to_string() + .contains("a row with 3 columns") + ); + let e = read_clingen(["chr1\t5\t5\tx\t0\t.\t5\t5\t255,111,0"], &mut c).unwrap_err(); + assert!(e.to_string().contains("x: empty span")); + let e = read_clingen(["chr1\t1\t5\tx\t0\t.\t1\t5\t1,2,3"], &mut c).unwrap_err(); + assert!(e.to_string().contains("itemRgb '1,2,3'")); + let e = read_clingen(["chr1\t1\t5\tx\t0\t.\t1\t5\t0,0,0"], &mut c).unwrap_err(); + assert!(e.to_string().contains("no region row")); + } + + #[test] + fn interval_helpers() { + assert_eq!( + merge(vec![(5, 8), (1, 3), (3, 4), (7, 10)]), + [(1, 4), (5, 10)] + ); + assert_eq!( + subtract(0, 20, &[(2, 4), (6, 8), (30, 40)]), + [(0, 2), (4, 6), (8, 20)] + ); + assert_eq!(subtract(3, 7, &[(0, 4), (6, 9)]), [(4, 6)]); + assert!(subtract(3, 7, &[(0, 9)]).is_empty()); + assert_eq!( + pieces(&[(0, 499), (0, 1_000_001)]), + [(0, 333_334), (333_334, 666_668), (666_668, 1_000_001)] + ); + assert_eq!(pieces(&[(0, 500)]), [(0, 500)]); + assert!(clear(4, 6, &[(0, 4), (6, 9)])); + assert!(!clear(3, 6, &[(0, 4), (6, 9)])); + assert_eq!(flank(0, 1_000_000, true, &[]), Some((850_000, 950_000))); + assert_eq!( + flank(0, 1_000_000, true, &[(900_000, 945_000)]), + Some((800_000, 900_000)) + ); + assert_eq!( + flank(0, 100_000, true, &[]), + None, + "past the chromosome start" + ); + } + + #[test] + fn locus_ids_and_labels() { + assert_eq!( + locus_id("22q11.2 recurrent (DGS) region (proximal, A-D) (includes TBX1)"), + "22q11.2_DGS_proximal_A-D" + ); + assert_eq!( + locus_id("Xp22.31 recurrent region (includes STS)"), + "Xp22.31" + ); + assert_eq!(locus_id("1q21.1 (a (b) c)"), "1q21.1_b"); + assert_eq!(label(" a;b\tc "), "a,b,c"); + } + + #[test] + fn a_repeated_locus_gets_a_number() { + let bed = "chr1\t10000000\t10100000\t1p (X) one\t0\t.\t0\t0\t255,111,0\nchr1\t20000000\t20100000\t1p (X) two\t0\t.\t0\t0\t255,111,0\n"; + let rows = small(bed).unwrap(); + assert_eq!((cell(&rows[0], 0), cell(&rows[1], 0)), ("1p_X", "1p_X_2")); + assert!(cell(&rows[1], 3).starts_with("1p_X_2_t1:")); + } + + #[test] + fn splice_keeps_the_head_and_check_counts_both_sides() { + let rows = vec!["a\t1".to_owned(), "b\t2".to_owned()]; + let old = b"# head\nx\n# --- interval rows: written by an older builder ---\nold\n"; + let new = splice(old, &rows); + assert_eq!( + new, + format!("# head\nx\n{MARKER}\na\t1\nb\t2\n").into_bytes() + ); + assert_eq!( + splice(b"x\n", &rows), + format!("x\n{MARKER}\na\t1\nb\t2\n").into_bytes() + ); + let added = [new.as_slice(), b"c\t3\n"].concat(); + assert_eq!(line_diff(&new, &added), (0, 1), "an added line is counted"); + assert_eq!(line_diff(&added, &new), (1, 0)); + let swapped = format!("# head\nx\n{MARKER}\nb\t2\na\t1\n").into_bytes(); + assert_eq!(line_diff(&new, &swapped), (0, 0)); + } +} diff --git a/rust/crates/og-dev/src/build/gtf.rs b/rust/crates/og-dev/src/build/gtf.rs new file mode 100644 index 000000000..89926e516 --- /dev/null +++ b/rust/crates/og-dev/src/build/gtf.rs @@ -0,0 +1,133 @@ +//! The GENCODE GTF reading the callability and SV rule builders share. +//! +//! A GTF line is read as bytes: its columns are found by their tabs, and an +//! attribute is found by a scan of the `key "value";` pairs of the ninth +//! column, only on the lines a builder keeps. GENCODE writes every pair as +//! `key "value"` or `key value` (a bare number), separated by `; `, and no +//! value holds a `;` or a `"`. + +use super::BuildError; + +/// One record of a GTF, by its columns. +pub(crate) struct Record<'a> { + pub seq: &'a [u8], + pub feature: &'a [u8], + start: &'a [u8], + end: &'a [u8], + pub attrs: &'a [u8], +} + +impl<'a> Record<'a> { + /// The record of a data line (without its line ending); `None` for a + /// comment. + pub fn parse(line: &'a [u8]) -> Result>, BuildError> { + if line.first() == Some(&b'#') { + return Ok(None); + } + let mut c = line.splitn(9, |&b| b == b'\t'); + let mut next = || { + c.next().ok_or_else(|| { + BuildError::Input(format!( + "GTF: a line with under 9 columns: {:?}", + String::from_utf8_lossy(line) + )) + }) + }; + let seq = next()?; + let _source = next()?; + let feature = next()?; + let start = next()?; + let end = next()?; + let (_score, _strand, _frame) = (next()?, next()?, next()?); + let attrs = next()?; + Ok(Some(Record { + seq, + feature, + start, + end, + attrs, + })) + } + + /// The 0-based start and the end. + pub fn span(&self) -> Result<(i64, i64), BuildError> { + Ok((number(self.start)? - 1, number(self.end)?)) + } + + /// The value of attribute `key`; the last one when it is given twice, + /// as a dict of the pairs keeps it. A value that is not UTF-8 is an + /// error. + pub fn attr(&self, key: &str) -> Result, BuildError> { + attr(self.attrs, key.as_bytes()) + .map(|v| { + std::str::from_utf8(v) + .map_err(|_| BuildError::Input(format!("GTF: {key} is not UTF-8: {v:?}"))) + }) + .transpose() + } +} + +/// A position column. +fn number(cell: &[u8]) -> Result { + std::str::from_utf8(cell) + .ok() + .and_then(|s| s.trim().parse().ok()) + .ok_or_else(|| { + BuildError::Input(format!( + "GTF: not a whole number: {:?}", + String::from_utf8_lossy(cell) + )) + }) +} + +/// The value of the last `key "value"` pair of an attribute column. +pub(crate) fn attr<'a>(attrs: &'a [u8], key: &[u8]) -> Option<&'a [u8]> { + let mut found = None; + for pair in attrs.split(|&b| b == b';') { + let pair = pair.trim_ascii_start(); + if let Some(rest) = pair.strip_prefix(key) + && let Some(quoted) = rest.strip_prefix(b" \"") + && let Some(end) = quoted.iter().position(|&b| b == b'"') + { + found = Some("ed[..end]); + } + } + found +} + +#[cfg(test)] +mod tests { + use super::*; + + const GENE: &[u8] = b"chr1\tHAVANA\tgene\t11\t20\t.\t+\t.\tgene_id \"ENSG1\"; gene_type \"protein_coding\"; gene_name \"A B\"; level 2; tag \"x\"; tag \"y\"; xgene_name \"no\";"; + + #[test] + fn a_gene_line_reads_its_columns_and_attributes() { + let r = Record::parse(GENE).unwrap().unwrap(); + assert_eq!((r.seq, r.feature), (&b"chr1"[..], &b"gene"[..])); + assert_eq!(r.span().unwrap(), (10, 20)); + assert_eq!(r.attr("gene_name").unwrap(), Some("A B"), "not xgene_name"); + assert_eq!( + r.attr("tag").unwrap(), + Some("y"), + "the last of a repeated key" + ); + assert_eq!( + r.attr("level").unwrap(), + None, + "an unquoted value is no pair" + ); + assert_eq!(r.attr("gene").unwrap(), None, "a key is matched whole"); + assert_eq!(r.attr("hgnc_id").unwrap(), None); + } + + #[test] + fn comments_pass_and_short_lines_stop() { + assert!(Record::parse(b"##description").unwrap().is_none()); + assert!(Record::parse(b"chr1\tx\tgene").is_err()); + let r = Record::parse(b"chr1\tx\tgene\tone\t2\t.\t+\t.\t") + .unwrap() + .unwrap(); + assert!(r.span().is_err()); + } +} diff --git a/rust/crates/og-dev/src/build/mod.rs b/rust/crates/og-dev/src/build/mod.rs new file mode 100644 index 000000000..b2098e3dc --- /dev/null +++ b/rust/crates/og-dev/src/build/mod.rs @@ -0,0 +1,298 @@ +//! `og dev build`: the builders of the asset tables OwnGenome commits. +//! +//! Each verb rebuilds one committed file from pinned public downloads that +//! a shell wrapper in `scripts/ci/` has already fetched and checked against +//! its sha256 (`scripts/ci/lib/fetch-pinned.sh`): `og` itself never +//! downloads. The committed file is the oracle, byte for byte. +//! +//! - `og dev build cn-intervals` ([`cn_intervals`]): the interval rows of +//! `bin/cn_loci.tsv` (`scripts/ci/build-cn-intervals.sh`). +//! - `og dev build callability` ([`callability`]): +//! `assets/callability_genes.tsv` (`scripts/ci/build-callability.sh`). +//! - `og dev build sv-rule-regions` ([`sv_rules`]): +//! `assets/sv_rule_regions.tsv` (`scripts/ci/build-sv-rule-regions.sh`). +//! +//! Exit 0 when the file was written (or `--check` found it right); 1 when +//! `--check` found it different, or when the callability or SV rule +//! downloads hold something their rules refuse (a CMRG gene with no GENCODE +//! gene, a symbol listed twice, not one MHC region), as the Python builders +//! exited; 2 for a usage error or an input that cannot be read, and for +//! every cn-intervals input the rules refuse, as its Python builder exited. + +pub mod callability; +pub mod cn_intervals; +mod gtf; +pub mod sv_rules; + +use og_io::{Lines, Utf8}; +use std::fmt; +use std::io::{self, BufRead, Write}; +use std::path::Path; + +const USAGE: &str = "usage: og dev build cn-intervals --clingen BED --segdups BED --gaps BED --satellites BED --parascopy TXT (--catalog FILE [--check] | --stdout)\n og dev build callability --work DIR --repo DIR --out FILE\n og dev build sv-rule-regions --work DIR --out FILE\n"; + +/// GRCh38's primary chromosomes in their sort order: chr1-22, X, Y. +pub(crate) const CHROMS: [&str; 24] = [ + "chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10", "chr11", + "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19", "chr20", "chr21", + "chr22", "chrX", "chrY", +]; + +/// The index of a primary chromosome in [`CHROMS`], from its name's bytes. +pub(crate) fn chrom_index(name: impl AsRef<[u8]>) -> Option { + let rest = name.as_ref().strip_prefix(b"chr")?; + match rest { + b"X" => Some(22), + b"Y" => Some(23), + [d @ b'1'..=b'9'] => Some(usize::from(d - b'1')), + [a @ b'1'..=b'2', b @ b'0'..=b'9'] => { + let n = usize::from(a - b'0') * 10 + usize::from(b - b'0'); + (n <= 22).then(|| n - 1) + } + _ => None, + } +} + +/// Why a build stopped: exit 1 for downloads the rules refuse, exit 2 for +/// an input that cannot be read. +#[derive(Debug)] +pub enum BuildError { + /// The inputs were read, but they hold what the rules refuse. + Data(String), + /// An input could not be read or parsed. + Input(String), +} + +impl BuildError { + fn code(&self) -> u8 { + match self { + BuildError::Data(_) => 1, + BuildError::Input(_) => 2, + } + } +} + +impl fmt::Display for BuildError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + BuildError::Data(m) | BuildError::Input(m) => f.write_str(m), + } + } +} + +/// An input that cannot be read, named by its path. +pub(crate) fn input_error(path: &Path, e: &io::Error) -> BuildError { + BuildError::Input(crate::py_os_error(e, &path.to_string_lossy())) +} + +/// The lines of a plain or gzip file, without their endings, as Python's +/// text mode reads them. +pub(crate) fn lines(path: &Path, utf8: Utf8) -> Result>, BuildError> { + let r = og_io::open(path).map_err(|e| input_error(path, &e))?; + Ok(Lines::new(r, utf8)) +} + +/// Reads each line of `path` into `f`, stopping at the first error. +pub(crate) fn each_line( + path: &Path, + utf8: Utf8, + mut f: impl FnMut(&str) -> Result<(), BuildError>, +) -> Result<(), BuildError> { + for line in lines(path, utf8)? { + f(&line.map_err(|e| input_error(path, &e))?)?; + } + Ok(()) +} + +/// Reads each line of `path` (plain or gzip) into `f` as bytes, without +/// its `\n` or `\r\n`, stopping at the first error. For the big files, +/// where a text line per line would cost more than the reading. +pub(crate) fn each_raw_line( + path: &Path, + mut f: impl FnMut(&[u8]) -> Result<(), BuildError>, +) -> Result<(), BuildError> { + let mut r = og_io::open(path).map_err(|e| input_error(path, &e))?; + let mut buf = Vec::with_capacity(1 << 12); + loop { + buf.clear(); + if r.read_until(b'\n', &mut buf) + .map_err(|e| input_error(path, &e))? + == 0 + { + return Ok(()); + } + let mut line = buf.as_slice(); + if let Some(l) = line.strip_suffix(b"\n") { + line = l.strip_suffix(b"\r").unwrap_or(l); + } + f(line)?; + } +} + +/// A whole number of a text column, as Python's `int()` reads one. +pub(crate) fn int(cell: &str, what: &str) -> Result { + pycompat::text::strip(cell, None) + .parse() + .map_err(|_| BuildError::Input(format!("{what}: not a whole number: {cell:?}"))) +} + +/// Takes `--flag VALUE` pairs and bare `--switch`es from `args`. +pub(crate) struct Args { + values: Vec<(String, String)>, + switches: Vec, +} + +impl Args { + pub(crate) fn parse( + args: &[String], + flags: &[&str], + switches: &[&str], + ) -> Result { + let mut out = Args { + values: Vec::new(), + switches: Vec::new(), + }; + let mut it = args.iter(); + while let Some(a) = it.next() { + if flags.contains(&a.as_str()) { + let v = it + .next() + .ok_or_else(|| format!("argument {a}: expected one argument"))?; + if out.values.iter().any(|(k, _)| k == a) { + return Err(format!("argument {a}: given twice")); + } + out.values.push((a.clone(), v.clone())); + } else if switches.contains(&a.as_str()) { + out.switches.push(a.clone()); + } else { + return Err(format!("unrecognized arguments: {a}")); + } + } + Ok(out) + } + + pub(crate) fn get(&self, flag: &str) -> Option<&str> { + self.values + .iter() + .find(|(k, _)| k == flag) + .map(|(_, v)| v.as_str()) + } + + pub(crate) fn need(&self, flag: &str) -> Result<&str, String> { + self.get(flag) + .ok_or_else(|| format!("the following arguments are required: {flag}")) + } + + pub(crate) fn has(&self, switch: &str) -> bool { + self.switches.iter().any(|s| s == switch) + } +} + +/// `og dev build VERB ARGS`; returns the exit code. +pub fn cli(args: &[String], out: &mut dyn Write, err: &mut dyn Write) -> u8 { + let result = match args.first().map(String::as_str) { + Some("cn-intervals") => cn_intervals::cli(&args[1..], out, err), + Some("callability") => callability::cli(&args[1..], err), + Some("sv-rule-regions") => sv_rules::cli(&args[1..], err), + Some("-h" | "--help") => { + let _ = out.write_all(USAGE.as_bytes()); + return 0; + } + Some(other) => Err(Failure::Usage(format!("unknown builder '{other}'"))), + None => Err(Failure::Usage( + "the following arguments are required: builder".to_owned(), + )), + }; + match result { + Ok(code) => code, + Err(Failure::Usage(m)) => { + let _ = writeln!(err, "{USAGE}og dev build: error: {m}"); + 2 + } + Err(Failure::Build(e)) => { + let _ = writeln!(err, "ERROR: {e}"); + e.code() + } + } +} + +/// How a builder's command line ended badly. +pub(crate) enum Failure { + Usage(String), + Build(BuildError), +} + +impl From for Failure { + fn from(e: BuildError) -> Self { + Failure::Build(e) + } +} + +impl From for Failure { + fn from(m: String) -> Self { + Failure::Usage(m) + } +} + +/// Writes `bytes` to `path` through a temporary file renamed over it. +pub(crate) fn write_atomic(path: &Path, bytes: &[u8]) -> Result<(), BuildError> { + let mut f = og_io::AtomicFile::create(path).map_err(|e| input_error(path, &e))?; + f.write_all(bytes).map_err(|e| input_error(path, &e))?; + f.commit().map_err(|e| input_error(path, &e)) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn s(v: &[&str]) -> Vec { + v.iter().map(|x| (*x).to_owned()).collect() + } + + #[test] + fn unknown_builder_is_a_usage_error() { + let (mut o, mut e) = (Vec::new(), Vec::new()); + assert_eq!(cli(&s(&["nope"]), &mut o, &mut e), 2); + assert!( + String::from_utf8(e) + .unwrap() + .contains("unknown builder 'nope'") + ); + assert_eq!(cli(&[], &mut o, &mut Vec::new()), 2); + } + + #[test] + fn args_refuse_a_repeated_or_unknown_flag() { + assert!(Args::parse(&s(&["--out", "a", "--out", "b"]), &["--out"], &[]).is_err()); + assert!(Args::parse(&s(&["--what"]), &["--out"], &[]).is_err()); + assert!(Args::parse(&s(&["--out"]), &["--out"], &[]).is_err()); + let a = Args::parse(&s(&["--out", "a", "--check"]), &["--out"], &["--check"]).unwrap(); + assert_eq!(a.get("--out"), Some("a")); + assert!(a.has("--check")); + } + + #[test] + fn chromosome_names() { + for (i, c) in CHROMS.iter().enumerate() { + assert_eq!(chrom_index(c), Some(i)); + } + for c in [ + "chr0", + "chr23", + "chr01", + "chrM", + "chr1_KI270706v1_random", + "1", + "chrx", + "", + ] { + assert_eq!(chrom_index(c), None, "{c}"); + } + } + + #[test] + fn int_reads_as_python_does() { + assert_eq!(int(" 12 ", "x").unwrap(), 12); + assert!(int("1.5", "x").is_err()); + } +} diff --git a/rust/crates/og-dev/src/build/sv_rules.rs b/rust/crates/og-dev/src/build/sv_rules.rs new file mode 100644 index 000000000..7004f02b7 --- /dev/null +++ b/rust/crates/og-dev/src/build/sv_rules.rs @@ -0,0 +1,649 @@ +//! `og dev build sv-rule-regions`: `assets/sv_rule_regions.tsv`, the GRCh38 +//! regions where `og sv postfilter` overrides AnnotSV's class 4 or 5. +//! +//! All on chr1-22, X and Y: +//! +//! - `tandem_repeat`: every TRExplorer v2.1 repeat of 50 bp or more (the +//! shortest SV) that overlaps a coding exon of a MANE Select transcript; +//! name `GENE:MOTIF` (a motif over 12 bases as its length); +//! - `ig_tr_locus`: the span of GENCODE's IG and TR genes of each of IGH, +//! IGK, IGL, TRA (with TRD inside it), TRB and TRG: per locus the largest +//! group of genes on one chromosome with gaps under 1 Mb, so the orphon +//! copies on other chromosomes stay out; +//! - `mhc`: the MHC region of NCBI's assembly regions report; +//! - `gene_alias`: the GENCODE span of each HGNC pseudogene whose alias or +//! previous symbol is the approved symbol of a gene with a GenCC +//! curation; name `ALIAS|PSEUDOGENE|its HGNC id|the HGNC id of the gene +//! that holds ALIAS as its symbol`. +//! +//! A gene's span is that of the first GENCODE record of its HGNC id, and the +//! `_PAR_Y` copies of the pseudoautosomal genes are skipped, as the +//! callability builder does: the Python this replaces kept the last record, +//! so a PAR pseudogene would have taken its chrY span (none does in the +//! pinned files, so the committed table is the same). +//! +//! The GTF and the catalog are streamed once. The MANE Select coding exons +//! of each chromosome are sorted by start with a running maximum of their +//! ends, so each repeat finds the exons it overlaps by binary search. + +use super::gtf::Record; +use super::{Args, BuildError, CHROMS, Failure, chrom_index, each_line, each_raw_line, int}; +use og_io::{Lines, Utf8}; +use pycompat::csv::{DictReader, Quoting}; +use std::collections::{BTreeMap, BTreeSet, HashMap}; +use std::io::{self, Write}; +use std::path::Path; + +/// The tag of a MANE Select transcript's records. +const MANE: &[u8] = b"tag \"MANE_Select\""; +const MIN_TR: i64 = 50; +const LOCUS_GAP: i64 = 1_000_000; +const SRC_TR: &str = "TRExplorer v2.1; GENCODE 50 MANE Select CDS"; +const SRC_IG: &str = "GENCODE 50 IG and TR genes"; +const SRC_MHC: &str = "GRCh38 assembly regions (GCA_000001405.15)"; +const SRC_ALIAS: &str = "HGNC 2026-10-06; GENCODE 50"; + +/// One row of the table. +struct Row { + chrom: usize, + start: i64, + end: i64, + rule: &'static str, + action: &'static str, + name: String, + source: &'static str, +} + +/// What the GTF gives. +#[derive(Default)] +struct Gtf { + /// HGNC id -> (chrom, start, end) of its first record. + span: HashMap, + /// IG or TR locus -> its genes as (chrom name, start, end). + ig: BTreeMap>, + /// Per chromosome, the MANE Select coding exons: (start, end, gene). + cds: Vec>, +} + +/// The IG or TR locus of a gene, from its type and symbol. +fn ig_locus(gene_type: &str, name: &str) -> Option { + let ig_type = gene_type.starts_with("IG_") || gene_type.starts_with("TR_"); + let b = name.as_bytes(); + let ig_name = b.len() >= 3 + && matches!( + (b[0], b[1], b[2]), + (b'I', b'G', b'H' | b'K' | b'L') | (b'T', b'R', b'A' | b'B' | b'D' | b'G') + ); + (ig_type && ig_name).then(|| { + if name.starts_with("TRD") { + "TRA".to_owned() + } else { + name[..3].to_owned() + } + }) +} + +fn read_gtf(path: &Path) -> Result { + let mut g = Gtf { + cds: vec![Vec::new(); CHROMS.len()], + ..Gtf::default() + }; + each_raw_line(path, |line| { + let Some(r) = Record::parse(line)? else { + return Ok(()); + }; + let Some(chrom) = chrom_index(r.seq) else { + return Ok(()); + }; + match r.feature { + b"gene" => { + let (start, end) = r.span()?; + if let Some(h) = r.attr("hgnc_id")?.filter(|h| !h.is_empty()) + && !r.attr("gene_id")?.is_some_and(|id| id.ends_with("_PAR_Y")) + { + g.span.entry(h.to_owned()).or_insert((chrom, start, end)); + } + let (name, ty) = ( + r.attr("gene_name")?.unwrap_or(""), + r.attr("gene_type")?.unwrap_or(""), + ); + if let Some(locus) = ig_locus(ty, name) { + g.ig.entry(locus) + .or_default() + .push((CHROMS[chrom], start, end)); + } + } + b"CDS" if r.attrs.windows(MANE.len()).any(|w| w == MANE) => { + let name = r.attr("gene_name")?.ok_or_else(|| { + BuildError::Input(format!( + "GTF: a MANE Select CDS with no gene_name: {:?}", + String::from_utf8_lossy(line) + )) + })?; + let (start, end) = r.span()?; + g.cds[chrom].push((start, end, name.to_owned())); + } + _ => {} + } + Ok(()) + })?; + Ok(g) +} + +/// The MANE Select coding exons of one chromosome, sorted, with the +/// running maximum of their ends. +struct Exons { + iv: Vec<(i64, i64, String)>, + max_end: Vec, +} + +impl Exons { + fn new(mut iv: Vec<(i64, i64, String)>) -> Exons { + iv.sort_unstable(); + let mut m = i64::MIN; + let max_end = iv + .iter() + .map(|x| { + m = m.max(x.1); + m + }) + .collect(); + Exons { iv, max_end } + } + + /// The genes of the exons that overlap `[s, e)`, sorted. + fn genes(&self, s: i64, e: i64) -> BTreeSet<&str> { + let mut out = BTreeSet::new(); + let mut j = self.iv.partition_point(|x| x.0 < e); + while j > 0 && self.max_end[j - 1] > s { + j -= 1; + if self.iv[j].1 > s { + out.insert(self.iv[j].2.as_str()); + } + } + out + } +} + +/// tandem_repeat rows: catalog repeats of `MIN_TR` bp or more on a MANE +/// Select coding exon. +fn tandem_repeats( + path: &Path, + cds: Vec>, + rows: &mut Vec, +) -> Result<(), BuildError> { + let index: Vec> = cds + .into_iter() + .map(|iv| (!iv.is_empty()).then(|| Exons::new(iv))) + .collect(); + each_raw_line(path, |line| { + let mut c = line.splitn(5, |&b| b == b'\t'); + let Some(ex) = c + .next() + .and_then(chrom_index) + .and_then(|i| index[i].as_ref().map(|x| (i, x))) + else { + return Ok(()); + }; + let short = || { + BuildError::Input(format!( + "{}: a row with under 4 columns: {:?}", + path.display(), + String::from_utf8_lossy(line) + )) + }; + let (s, e) = (c.next().ok_or_else(short)?, c.next().ok_or_else(short)?); + let (s, e) = (bytes_int(s)?, bytes_int(e)?); + if e - s < MIN_TR { + return Ok(()); + } + let genes = ex.1.genes(s, e); + if !genes.is_empty() { + let motif = std::str::from_utf8(c.next().ok_or_else(short)?).map_err(|_| { + BuildError::Input(format!("{}: a motif that is not UTF-8", path.display())) + })?; + let n = motif.chars().count(); + let motif = if n <= 12 { + motif.to_owned() + } else { + format!("{n}bp motif") + }; + rows.push(Row { + chrom: ex.0, + start: s, + end: e, + rule: "tandem_repeat", + action: "class_3", + name: format!( + "{}:{motif}", + genes.into_iter().collect::>().join("/") + ), + source: SRC_TR, + }); + } + Ok(()) + }) +} + +/// A position of the catalog. +fn bytes_int(cell: &[u8]) -> Result { + int(&String::from_utf8_lossy(cell), "catalog position") +} + +/// ig_tr_locus rows: per locus the largest group of genes on one +/// chromosome with gaps under `LOCUS_GAP` (the first such group on a tie). +fn ig_loci(ig: BTreeMap>, rows: &mut Vec) { + for (locus, mut genes) in ig { + genes.sort_unstable(); + // (chrom, start, end, genes) + let mut groups: Vec<(&str, i64, i64, usize)> = Vec::new(); + for (chrom, s, e) in genes { + match groups.last_mut() { + Some(g) if g.0 == chrom && s - g.2 < LOCUS_GAP => { + g.2 = g.2.max(e); + g.3 += 1; + } + _ => groups.push((chrom, s, e, 1)), + } + } + let mut best = groups[0]; + for g in &groups[1..] { + if g.3 > best.3 { + best = *g; + } + } + rows.push(Row { + chrom: chrom_index(best.0).expect("only chr1-22, X and Y are read"), + start: best.1, + end: best.2, + rule: "ig_tr_locus", + action: "class_3", + name: locus, + source: SRC_IG, + }); + } +} + +/// The mhc row: the assembly report's one MHC region (1-based, inclusive). +fn mhc(path: &Path) -> Result { + let mut found: BTreeSet> = BTreeSet::new(); + each_line(path, Utf8::Strict, |line| { + if line.starts_with("MHC\t") { + found.insert(line.split('\t').take(4).map(str::to_owned).collect()); + } + Ok(()) + })?; + let one = match found.len() { + 1 => found.pop_first().unwrap_or_default(), + n => { + return Err(BuildError::Data(format!( + "want one MHC region in the assembly report, got {n}: {found:?}" + ))); + } + }; + let bad = || { + BuildError::Input(format!( + "{}: an MHC row with under 4 columns: {one:?}", + path.display() + )) + }; + let (chrom, s, e) = ( + one.get(1).ok_or_else(bad)?, + one.get(2).ok_or_else(bad)?, + one.get(3).ok_or_else(bad)?, + ); + let name = format!("chr{chrom}"); + Ok(Row { + chrom: chrom_index(&name).ok_or_else(|| { + BuildError::Input(format!("{}: the MHC is on {name}", path.display())) + })?, + start: int(s, "MHC start")? - 1, + end: int(e, "MHC end")?, + rule: "mhc", + action: "class_3", + name: "MHC".to_owned(), + source: SRC_MHC, + }) +} + +/// The lines of a file with their endings, as a file opened with +/// `newline=""` hands them to `csv`. +struct RawLines(Lines>); + +impl Iterator for RawLines { + type Item = io::Result; + fn next(&mut self) -> Option { + let mut s = String::new(); + match self.0.read_line(&mut s) { + Ok(true) => Some(Ok(s)), + Ok(false) => None, + Err(e) => Some(Err(e)), + } + } +} + +/// gene_alias rows: a pseudogene that carries a disease gene's symbol as an +/// alias or previous symbol. +fn gene_aliases( + path: &Path, + span: &HashMap, + rows: &mut Vec, +) -> Result<(), BuildError> { + let csv_error = |e: pycompat::csv::Error| BuildError::Input(format!("{}: {e}", path.display())); + // symbol -> (HGNC id, has a GenCC curation); the last row of a symbol. + let mut approved: HashMap = HashMap::new(); + // (HGNC id, symbol, the alias and previous symbols) of the pseudogenes + // GENCODE places. + let mut pseudo: Vec<(String, String, BTreeSet)> = Vec::new(); + let lines = RawLines(super::lines(path, Utf8::Strict)?); + for row in DictReader::new(lines, '\t', Quoting::Minimal) { + let row = row.map_err(csv_error)?; + let get = |k: &str| { + row.get(k).flatten().ok_or_else(|| { + BuildError::Input(format!("{}: a row with no {k} cell", path.display())) + }) + }; + let (symbol, hgnc) = (get("symbol")?, get("hgnc_id")?); + approved.insert( + symbol.to_owned(), + (hgnc.to_owned(), !get("gencc")?.is_empty()), + ); + if get("locus_group")? == "pseudogene" && span.contains_key(hgnc) { + let names = format!("{}|{}", get("alias_symbol")?, get("prev_symbol")?); + let syms = names + .split('|') + .map(|x| x.trim_matches('"')) + .filter(|x| !x.is_empty()) + .map(str::to_owned) + .collect(); + pseudo.push((hgnc.to_owned(), symbol.to_owned(), syms)); + } + } + for (hgnc, symbol, syms) in pseudo { + for sym in syms { + if let Some((other, true)) = approved.get(&sym) + && *other != hgnc + { + let (chrom, start, end) = span[&hgnc]; + rows.push(Row { + chrom, + start, + end, + rule: "gene_alias", + action: "remap_gene", + name: format!("{sym}|{symbol}|{hgnc}|{other}"), + source: SRC_ALIAS, + }); + } + } + } + Ok(()) +} + +/// The table's text and its row counts line, from the downloads in `work`. +pub fn build(work: &Path) -> Result<(String, String), BuildError> { + let gtf = read_gtf(&work.join("gencode.gtf.gz"))?; + let mut rows = Vec::new(); + tandem_repeats(&work.join("trx.bed.gz"), gtf.cds, &mut rows)?; + ig_loci(gtf.ig, &mut rows); + rows.push(mhc(&work.join("regions.txt"))?); + gene_aliases(&work.join("hgnc.txt"), >f.span, &mut rows)?; + rows.sort_by(|a, b| { + (a.chrom, a.start, a.end, a.rule, &a.name).cmp(&(b.chrom, b.start, b.end, b.rule, &b.name)) + }); + let mut counts: BTreeMap<&str, usize> = BTreeMap::new(); + for r in &rows { + *counts.entry(r.rule).or_default() += 1; + } + let counts: Vec = counts.iter().map(|(k, n)| format!("{k} {n}")).collect(); + let summary = format!("Rows: {}", counts.join(", ")); + let mut text = String::from( + "# assets/sv_rule_regions.tsv: the GRCh38 regions where og sv postfilter overrides AnnotSV's class 4\n\ + # or 5. Written by scripts/ci/build-sv-rule-regions.sh, which says how each rule's rows are chosen;\n\ + # rebuild it there, do not edit it by hand. start is 0-based and end 1-based, as in a BED file.\n", + ); + text.push_str(&format!( + "# {summary}.\n#chrom\tstart\tend\trule\taction\tname\tsource\n" + )); + for r in &rows { + text.push_str(&format!( + "{}\t{}\t{}\t{}\t{}\t{}\t{}\n", + CHROMS[r.chrom], r.start, r.end, r.rule, r.action, r.name, r.source + )); + } + Ok((text, summary)) +} + +pub(super) fn cli(args: &[String], err: &mut dyn Write) -> Result { + let a = Args::parse(args, &["--work", "--out"], &[])?; + let (work, out) = (a.need("--work")?, a.need("--out")?); + let (text, summary) = build(Path::new(work))?; + super::write_atomic(Path::new(out), text.as_bytes())?; + let _ = writeln!(err, "{summary}"); + Ok(0) +} + +#[cfg(test)] +mod tests { + use super::*; + use std::fs; + + fn gz(text: &str) -> Vec { + let mut g = flate2::write::GzEncoder::new(Vec::new(), flate2::Compression::fast()); + g.write_all(text.as_bytes()).unwrap(); + g.finish().unwrap() + } + + fn gene(seq: &str, s: u32, e: u32, id: &str, ty: &str, name: &str, h: &str) -> String { + format!( + "{seq}\tX\tgene\t{s}\t{e}\t.\t+\t.\tgene_id \"{id}\"; gene_type \"{ty}\"; gene_name \"{name}\"; hgnc_id \"{h}\";\n" + ) + } + + fn cds(seq: &str, s: u32, e: u32, name: &str, mane: bool) -> String { + let tag = if mane { " tag \"MANE_Select\";" } else { "" }; + format!("{seq}\tX\tCDS\t{s}\t{e}\t.\t+\t0\tgene_id \"x\"; gene_name \"{name}\";{tag}\n") + } + + /// A work folder with a few rows of each rule; `par_first` puts the + /// _PAR_Y copy of the pseudogene PSP before its chrX record. + fn fixture(dir: &Path, par_first: bool) { + fs::create_dir_all(dir).unwrap(); + let px = gene( + "chrX", + 1001, + 2000, + "GP", + "unprocessed_pseudogene", + "PSP", + "HGNC:30", + ); + let py = gene( + "chrY", + 1001, + 2000, + "GP_PAR_Y", + "unprocessed_pseudogene", + "PSP", + "HGNC:30", + ); + let par = if par_first { + [py.clone(), px.clone()] + } else { + [px, py] + }; + let gtf = [ + "##description: test\n".to_owned(), + gene("chr1", 1001, 3000, "G1", "protein_coding", "AAA", "HGNC:1"), + cds("chr1", 1101, 1300, "AAA", true), + cds("chr1", 1251, 1400, "AAB", true), + cds("chr1", 2001, 2400, "AAA", false), + gene( + "chr14", + 100_001, + 100_500, + "I1", + "IG_V_gene", + "IGHV1", + "HGNC:10", + ), + gene( + "chr14", + 900_001, + 900_500, + "I2", + "IG_C_gene", + "IGHM", + "HGNC:11", + ), + gene( + "chr15", + 1, + 500, + "I3", + "IG_V_pseudogene", + "IGHV1OR15", + "HGNC:12", + ), + gene( + "chr14", + 20_001, + 20_500, + "T1", + "TR_V_gene", + "TRDV1", + "HGNC:13", + ), + gene( + "chr14", + 10_001, + 10_500, + "T2", + "TR_V_gene", + "TRAV1", + "HGNC:14", + ), + gene( + "chr2", + 5001, + 6000, + "G2", + "processed_pseudogene", + "PSA", + "HGNC:20", + ), + par[0].clone(), + par[1].clone(), + ] + .concat(); + fs::write(dir.join("gencode.gtf.gz"), gz(>f)).unwrap(); + fs::write( + dir.join("trx.bed.gz"), + gz("chr1\t1000\t1049\tCAG\nchr1\t1200\t1260\tCAG\nchr1\t1280\t1340\tAAAAAAAAAAAAAT\nchr1\t2100\t2200\tGC\nchrM\t1\t100\tA\n"), + ) + .unwrap(); + fs::write( + dir.join("regions.txt"), + "# Region-Name\tChromosome\tStart\tStop\tRole\nMHC\t6\t28510120\t33480577\tx\r\nMHC\t6\t28510120\t33480577\tx\r\nOTHER\t1\t1\t2\tx\n", + ) + .unwrap(); + fs::write( + dir.join("hgnc.txt"), + "hgnc_id\tsymbol\tlocus_group\talias_symbol\tprev_symbol\tgencc\n\ + HGNC:1\tAAA\tprotein-coding gene\t\t\t\"curated\"\n\ + HGNC:2\tBBB\tprotein-coding gene\t\t\t\n\ + HGNC:3\tCCC\tprotein-coding gene\t\t\tyes\n\ + HGNC:20\tPSA\tpseudogene\t\"AAA|BBB\"\tPSA\t\n\ + HGNC:30\tPSP\tpseudogene\tCCC\t\t\n\ + HGNC:40\tPSZ\tpseudogene\tAAA\t\t\n", + ) + .unwrap(); + } + + #[test] + fn each_rule_gives_its_rows() { + let dir = std::env::temp_dir().join(format!("og-svrules-{}", std::process::id())); + fixture(&dir, false); + let (text, summary) = build(&dir).unwrap(); + fs::remove_dir_all(&dir).unwrap(); + let body: Vec<&str> = text.lines().skip(4).collect(); + assert_eq!( + body, + [ + "#chrom\tstart\tend\trule\taction\tname\tsource", + "chr1\t1200\t1260\ttandem_repeat\tclass_3\tAAA/AAB:CAG\tTRExplorer v2.1; GENCODE 50 MANE Select CDS", + "chr1\t1280\t1340\ttandem_repeat\tclass_3\tAAA/AAB:14bp motif\tTRExplorer v2.1; GENCODE 50 MANE Select CDS", + "chr2\t5000\t6000\tgene_alias\tremap_gene\tAAA|PSA|HGNC:20|HGNC:1\tHGNC 2026-10-06; GENCODE 50", + "chr6\t28510119\t33480577\tmhc\tclass_3\tMHC\tGRCh38 assembly regions (GCA_000001405.15)", + "chr14\t10000\t20500\tig_tr_locus\tclass_3\tTRA\tGENCODE 50 IG and TR genes", + "chr14\t100000\t900500\tig_tr_locus\tclass_3\tIGH\tGENCODE 50 IG and TR genes", + "chrX\t1000\t2000\tgene_alias\tremap_gene\tCCC|PSP|HGNC:30|HGNC:3\tHGNC 2026-10-06; GENCODE 50", + ] + ); + assert_eq!( + summary, + "Rows: gene_alias 2, ig_tr_locus 2, mhc 1, tandem_repeat 2" + ); + assert!(text.contains("\n# Rows: gene_alias 2, ig_tr_locus 2, mhc 1, tandem_repeat 2.\n")); + } + + #[test] + fn a_par_pseudogene_keeps_its_chrx_span_whichever_record_comes_first() { + let dir = std::env::temp_dir().join(format!("og-svrules-par-{}", std::process::id())); + fixture(&dir, true); + let (text, _) = build(&dir).unwrap(); + fs::remove_dir_all(&dir).unwrap(); + assert!( + text.contains("\nchrX\t1000\t2000\tgene_alias\tremap_gene\tCCC|PSP|HGNC:30|HGNC:3\t"), + "{text}" + ); + assert!(!text.contains("\nchrY\t"), "{text}"); + } + + #[test] + fn two_mhc_regions_stop_the_build() { + let dir = std::env::temp_dir().join(format!("og-svrules-mhc-{}", std::process::id())); + fixture(&dir, false); + fs::write( + dir.join("regions.txt"), + "MHC\t6\t1\t2\tx\nMHC\t6\t1\t3\tx\n", + ) + .unwrap(); + let e = build(&dir).err().unwrap(); + fs::remove_dir_all(&dir).unwrap(); + assert!(matches!(e, BuildError::Data(_))); + assert!( + e.to_string() + .starts_with("want one MHC region in the assembly report, got 2"), + "{e}" + ); + } + + #[test] + fn exon_lookup_finds_every_overlap() { + let ex = Exons::new(vec![ + (0, 100, "A".into()), + (10, 20, "B".into()), + (50, 60, "C".into()), + (200, 300, "D".into()), + ]); + assert_eq!( + ex.genes(15, 55).into_iter().collect::>(), + ["A", "B", "C"] + ); + assert_eq!(ex.genes(100, 200).len(), 0, "half-open on both sides"); + assert_eq!( + ex.genes(99, 201).into_iter().collect::>(), + ["A", "D"] + ); + } + + #[test] + fn ig_loci_by_type_and_symbol() { + assert_eq!(ig_locus("IG_V_gene", "IGKV1").as_deref(), Some("IGK")); + assert_eq!(ig_locus("TR_D_gene", "TRDD1").as_deref(), Some("TRA")); + assert_eq!(ig_locus("protein_coding", "IGHM"), None); + assert_eq!(ig_locus("TR_V_gene", "TRX1"), None); + assert_eq!(ig_locus("IG_V_gene", "IG"), None); + } +} diff --git a/rust/crates/og-dev/src/lib.rs b/rust/crates/og-dev/src/lib.rs index 55c93fdf0..88c683e45 100644 --- a/rust/crates/og-dev/src/lib.rs +++ b/rust/crates/og-dev/src/lib.rs @@ -12,7 +12,10 @@ //! the report carries the report-only error strategy. //! - `og dev remnants` ([`remnants`]): does each Python file kept only for //! its importers (`tests/parity/remnants.tsv`) hold nothing else? +//! - `og dev build cn-intervals|callability|sv-rule-regions` ([`build`]): +//! rebuild a committed asset table from its pinned downloads. +pub mod build; pub mod contract; pub mod parity; pub mod remnants; @@ -22,7 +25,7 @@ mod zip; use std::io::Write; -const USAGE: &str = "usage: og dev parity-compare KIND ORACLE NEW\n og dev parity-compare --self-test\n og dev contract --static [--repo DIR] [--contract FILE]\n og dev contract --outputs DIR --leg bash|nextflow --sample S [--observe]\n og dev schema SCHEMA.json DOC.json\n og dev report-only [--repo DIR]\n og dev remnants [--list FILE] [--root DIR]\n"; +const USAGE: &str = "usage: og dev parity-compare KIND ORACLE NEW\n og dev parity-compare --self-test\n og dev contract --static [--repo DIR] [--contract FILE]\n og dev contract --outputs DIR --leg bash|nextflow --sample S [--observe]\n og dev schema SCHEMA.json DOC.json\n og dev report-only [--repo DIR]\n og dev remnants [--list FILE] [--root DIR]\n og dev build cn-intervals|callability|sv-rule-regions [args]\n"; /// Python's `str(OSError)`: `[Errno 2] No such file or directory: 'path'`. pub(crate) fn py_os_error(e: &std::io::Error, path: &str) -> String { @@ -57,6 +60,7 @@ pub fn main(args: &[String]) -> u8 { Some("contract") => run(contract::cli, &args[1..]), Some("schema") => run(schema::cli, &args[1..]), Some("report-only") => run(report_only::cli, &args[1..]), + Some("build") => run(build::cli, &args[1..]), Some("-h" | "--help") => { print!("{USAGE}"); 0 diff --git a/rust/og/src/main.rs b/rust/og/src/main.rs index 234ea7cac..a34aa302f 100644 --- a/rust/og/src/main.rs +++ b/rust/og/src/main.rs @@ -74,6 +74,7 @@ const USAGE: &str = "usage: og [args] og exomiser report [args] og lookup query|resolve|report [args] og dev parity-compare|contract|schema|report-only|remnants [args] + og dev build cn-intervals|callability|sv-rule-regions [args] og dev freshness [--out FILE] [--sections a,b] | --update-issue BODY --repo O/N | --self-test og version og --tree-hash From f47eb8a998025a300368b015e69ca5d59910bd7f Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 12:27:46 +0100 Subject: [PATCH 2/7] build(ci): switch the asset builders to og dev build, one fetch rule build-callability.sh and build-sv-rule-regions.sh keep their downloads and hand the folder to og dev build; their Python heredocs are gone. build-cn-intervals.py becomes build-cn-intervals.sh with the same --list, --cache and --check. The three share scripts/ci/lib/fetch-pinned.sh: a cached file with the right sha256 is read in place, a 403 or 404 fails at once, anything else is retried for about five minutes. The marker line of bin/cn_loci.tsv names the new script; no row changes. The cn-intervals case, the e2e data fill and the converter checks of tests/test_copy_number.py call the new script and og. Bead: pgp-bs9u.35 --- CHANGELOG.md | 6 + bin/cn_loci.tsv | 2 +- bin/copy_number.py | 2 +- docs/cn-loci.md | 6 +- docs/rust.md | 2 +- docs/testing.md | 2 +- scripts/ci/build-callability.sh | 218 +------------- scripts/ci/build-cn-intervals.py | 423 ---------------------------- scripts/ci/build-cn-intervals.sh | 97 +++++++ scripts/ci/build-sv-rule-regions.sh | 171 ++--------- scripts/ci/e2e-case-paths.tsv | 6 +- scripts/ci/e2e-run.sh | 2 +- scripts/ci/lib/fetch-pinned.sh | 61 ++++ tests/e2e/cn-intervals-1-catalog.sh | 10 +- tests/test_copy_number.py | 77 ++--- tests/test_e2e_data_cache.sh | 2 +- versions.env | 2 +- 17 files changed, 255 insertions(+), 834 deletions(-) delete mode 100755 scripts/ci/build-cn-intervals.py create mode 100755 scripts/ci/build-cn-intervals.sh create mode 100644 scripts/ci/lib/fetch-pinned.sh diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e21d1ac8..7ce7f1507 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,12 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n ## Unreleased +- 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; any other failure is tried for about five minutes, 30 s apart. `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. + - `build-cn-intervals.sh --check` counts the lines only in the file and the lines only in the rebuild. The Python counted only the lines the rebuild lacks, so a catalog with an added line said "0 line(s) differ". + - The SV rule builder takes each gene's span from the first GENCODE record of its HGNC id and skips the `_PAR_Y` copies, as the callability builder does. The Python kept the last record, so a pseudogene in a pseudoautosomal region would have taken its chrY span. No gene of the pinned files is one, so the table does not change. + - 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/bin/cn_loci.tsv b/bin/cn_loci.tsv index aa1b5d4c6..6f63e7cf7 100644 --- a/bin/cn_loci.tsv +++ b/bin/cn_loci.tsv @@ -66,7 +66,7 @@ HP gene chr16 HP_UNIT1:72056529-72057466:HP exons 3 and 4, the first copy of the LCE3BC interval chr1 LCE3BC_DEL:152600233-152614098:LCE3C to LCE3B, inside the LCE3C_LCE3B deletion LCE3_FLANK_L:152520000-152560000:flank on the LCE3D side, outside the deletion;LCE3_FLANK_R:152640000-152680000:flank on the LCE3A side, outside the deletion . sum flanks 2 0 .,-Q 1 . . 8 1.5 round:0.3 . 0=LCE3B and LCE3C deleted on both chromosomes (LCE3C_LCE3B-del on both), a common genotype that goes with a higher psoriasis risk|1=LCE3B and LCE3C deleted on one chromosome (LCE3C_LCE3B-del on one)|2=no LCE3B/LCE3C deletion seen|more={n} copies of the LCE3B/LCE3C stretch: more than two, a duplication bedcov cn-depth none 19169253 CCL3L1 gene chr17 CCL3L3:36194868-36196758:CCL3L3, the one copy of CCL3L1/CCL3L3 on the reference (ENST00000619989) CCL3L_FLANK_L:36112000-36148000:single-copy sequence between the CCL3 and CCL3L3 duplications;CCL3L_FLANK_R:36460000-36500000:single-copy sequence past the duplications . sum flanks 2 0 .,-Q 1 . . 8 1.5 round:0.3 . 0=no copy of CCL3L1 or CCL3L3|1=one copy of CCL3L1/CCL3L3 in all|more={n} copies of CCL3L1/CCL3L3 in all: the count differs widely between people, and no reading of it is given bedcov cn-depth none 15637236 CYP21A2 gene chr6 CYP21A2_BODY:32038950-32041644:CYP21A2 from the end of the C4A model to its last base (NM_000500.9);CYP21A1P_BODY:32006220-32008909:the same stretch of CYP21A1P (NR_040090.1) RCCX_FLANK_L:31945000-31975000:single-copy sequence before the RCCX duplication;RCCX_FLANK_R:32070000-32100000:single-copy sequence in TNXB past the duplication PSV_In2G:32039080-32039081:32006352 AC>G c.293-13A/C>G (In2G, intron 2);PSV_I172N:32039425-32039426:32006689 T>A c.518T>A (p.I172N, exon 4);PSV_V237E:32039809-32039810:32007073 T>A c.713T>A (p.V237E, exon 6 cluster);PSV_V281L:32040109-32040110:32007373 G>T c.844G>T (p.V281L, exon 7);PSV_Q318X:32040420-32040421:32007685 C>T c.955C>T (p.Q318X, exon 8);PSV_R356W:32040534-32040535:32007799 C>T c.1069C>T (p.R356W, exon 8) sum flanks 2 0 .,-Q 1 0.02 . 8 1.5 psv:0.4 . none=no sign of the 30 kb chimeric deletion or of a pseudogene conversion at the six sites read; short reads see neither a change elsewhere in CYP21A2 nor two copies on one chromosome and none on the other|deletion=one RCCX module fewer and one CYP21A2 copy fewer at the differing sites: the pattern of the 30 kb CYP21A1P/CYP21A2 chimeric deletion, a common cause of carrier status for 21-hydroxylase deficiency (congenital adrenal hyperplasia); unconfirmed: a clinical test (MLPA or long-range PCR) confirms it|conversion=fewer CYP21A2 copies carry the gene's own base at {sites} than at the other sites: a possible conversion of one CYP21A2 copy to the pseudogene sequence there, a carrier pattern for 21-hydroxylase deficiency (congenital adrenal hyperplasia) when it sits on one chromosome; which chromosome carries it is not known; unconfirmed: a clinical test confirms it|pseudogene_deletion=one RCCX module fewer with CYP21A2 kept at every site read: a CYP21A1P (pseudogene) deletion, common and not a carrier finding bedcov cn-depth paraphase:rccx 34285390,27228465 -# --- interval rows: written by scripts/ci/build-cn-intervals.py; do not edit by hand --- +# --- interval rows: written by scripts/ci/build-cn-intervals.sh; do not edit by hand --- 1q21.1_TAR_syndrome_proximal_BP2-BP3 interval chr1 1q21.1_TAR_syndrome_proximal_BP2-BP3_t1:145686999-145809382:1q21.1 recurrent (TAR syndrome) region (proximal, BP2-BP3) (includes RBM8A);1q21.1_TAR_syndrome_proximal_BP2-BP3_t2:145823828-145923649;1q21.1_TAR_syndrome_proximal_BP2-BP3_t3:145927476-146048495 1q21.1_TAR_syndrome_proximal_BP2-BP3_L:143034587-143134587:flank left of the region;1q21.1_TAR_syndrome_proximal_BP2-BP3_R:149926380-150026380:flank right of the region . pooled flanks 2 0 . . . 8 1.5 round:0.3 . . bedcov cn-depth cnvpytor:cnvs . 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4 interval chr1 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t1:145687000-145809382:1q21.1 recurrent (TAR syndrome) region (proximal and distal, BP1-BP4) (includes RBM8A and GJA5);1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t2:145823828-145923649;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t3:145927476-146048496;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t4:147105903-147473859;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t5:147473859-147841815;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t6:147842996-147917509 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_L:143034587-143134587:flank left of the region;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_R:149926380-150026380:flank right of the region . pooled flanks 2 0 . . . 8 1.5 round:0.3 . . bedcov cn-depth cnvpytor:cnvs . 1q21.1_distal_BP3-BP4 interval chr1 1q21.1_distal_BP3-BP4_t1:147105904-147473860:1q21.1 recurrent region (distal, BP3-BP4) (includes GJA5);1q21.1_distal_BP3-BP4_t2:147473860-147841815;1q21.1_distal_BP3-BP4_t3:147842996-147917509 1q21.1_distal_BP3-BP4_L:143034587-143134587:flank left of the region;1q21.1_distal_BP3-BP4_R:149926380-150026380:flank right of the region . pooled flanks 2 0 . . . 8 1.5 round:0.3 . . bedcov cn-depth cnvpytor:cnvs . diff --git a/bin/copy_number.py b/bin/copy_number.py index 2e1b95485..9412104b9 100755 --- a/bin/copy_number.py +++ b/bin/copy_number.py @@ -31,7 +31,7 @@ "call", "interval", "readings", "input", "owner", "cross_check", "pmids"] KINDS = ("gene", "interval", "organelle") # GRCh38 names of the primary chromosomes (the ClinGen BED spells chrX 'chrx': -# scripts/ci/build-cn-intervals.py normalises it, and a row with any other +# og dev build cn-intervals normalises it, and a row with any other # spelling is refused here). CHROMS = {f"chr{c}" for c in list(range(1, 23)) + ["X", "Y", "M"]} INPUTS = ("bedcov", "cnvpytor-level") diff --git a/docs/cn-loci.md b/docs/cn-loci.md index eb6f0a44e..a9b918f5f 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -89,7 +89,7 @@ A site is low when its estimate is more than 0.6 copies under the larger of the ### The interval rows -`scripts/ci/build-cn-intervals.py` writes these rows, everything in `bin/cn_loci.tsv` after its marker line; do not edit them by hand. Step 35 (and the PARASCOPY module) runs `og cn bed --kind interval`, one `samtools bedcov` over every segment of every row, then `og cn intervals`, which writes one line per row to `_cn_intervals.tsv`: locus, name, chrom, start, end, status, copies, expected, change (`deletion`, `duplication`, `no change`, or `.` when not called), estimate and reason. The step log lists the rows with a change. The rows give a total copy count only, never breakpoints. +`scripts/ci/build-cn-intervals.sh` writes these rows with `og dev build cn-intervals`, everything in `bin/cn_loci.tsv` after its marker line; do not edit them by hand. Step 35 (and the PARASCOPY module) runs `og cn bed --kind interval`, one `samtools bedcov` over every segment of every row, then `og cn intervals`, which writes one line per row to `_cn_intervals.tsv`: locus, name, chrom, start, end, status, copies, expected, change (`deletion`, `duplication`, `no change`, or `.` when not called), estimate and reason. The step log lists the rows with a change. The rows give a total copy count only, never breakpoints. - **The ClinGen recurrent regions**: every region row of ClinGen Dosage Sensitivity's recurrent CNV BED, version 2.1, GRCh38 (CC0, pinned by `CLINGEN_RECURRENT_CNV_URL` and its sha256 in `versions.env`). The file's breakpoint and LCR rows are the segmental duplications that bound the regions. They are cut out of the targets, as are GIAB v3.6 segmental duplications, assembly gaps, satellites and the regions of every pinned Parascopy model (CHRNA7 in 15q13.3 stays with Parascopy). Targets are cut to pieces of at most 500 kb, and their depth is pooled. Controls are two 100 kb windows, one on each side of the cluster of ClinGen rows the region belongs to, clear of all of that and of the targets of the hand rows below (so no row's control sits on another row's target); a side with no clear window within 3 Mb (a chromosome end, as next to 22q11.2 and 15q11.2) is left out. The BED spells its chrX rows `chrx`; the script normalises them, and the engine refuses any contig name that is not a GRCh38 chromosome. The classification of a CNV found here follows the ACMG and ClinGen technical standards (Riggs 2020, PMID 31690835); a call from this step is a lead for a clinical test, not a classification. - **Two more ClinGen curations**: the 1p36 terminal region (ISCA-37434, haploinsufficiency score 3) and the Xq28 region that holds MECP2 (ISCA-46304, triplosensitivity score 3: the MECP2 duplication), from ClinGen's region curation list, read on 2026-11-04. @@ -100,7 +100,7 @@ A row on chrX or chrY is read against the sample's own count there: its controls The couple screen (step 39) reads the table with `couple_report.py screen --cn-intervals ID=_cn_intervals.tsv`: each carrier allele is present, absent or not assessable for each partner, and a male partner's AZF deletions are listed the same way (present means the deletion is there). With `--cnvpytor ID=_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call is printed beside it. Step 39 and the Nextflow couple screen pass both when the partner has them. Both reports show the table in step 35's section: the counts by status and each row with a deletion or a duplication. -To change a row, change the script (its hand rows or a pinned download) and run `scripts/ci/build-cn-intervals.py`; `--check` fails when the catalog is not what the downloads give. +To change a row, change its rule or hand row in `rust/crates/og-dev/src/build/cn_intervals.rs`, or a pinned download in `scripts/ci/build-cn-intervals.sh`, rebuild og (`scripts/setup.sh --og`) and run `scripts/ci/build-cn-intervals.sh`; `--check` fails when the catalog is not what the downloads give. ## One owner per locus @@ -118,7 +118,7 @@ pypgx owns the genes it calls from the BAM, the `bam` rows of `assets/pgx_genes. 1. Check that no pinned caller owns the gene: Parascopy (`assets/parascopy_models_v1.7.txt`), pypgx, Paraphase. 2. Choose single-copy controls near the gene, outside any common deletion or duplication there and outside segmental duplications. -3. Add the row to `bin/cn_loci.tsv`, above its marker line, and cite the method in `pmids`. An interval row (everything after the marker) is never added by hand: add it to `scripts/ci/build-cn-intervals.py` and rerun the script. +3. Add the row to `bin/cn_loci.tsv`, above its marker line, and cite the method in `pmids`. An interval row (everything after the marker) is never added by hand: add it to the hand rows of `rust/crates/og-dev/src/build/cn_intervals.rs` and rerun `scripts/ci/build-cn-intervals.sh`. 4. Run it from a step: `og cn bed --locus NAME`, `samtools bedcov` with the row's passes, then `og cn call --locus NAME --bedcov ... --table ...`. 5. Add made-up bedcov tables for the row to `tests/test_copy_number.py`, at every count it can call and at each guard. diff --git a/docs/rust.md b/docs/rust.md index cb6924559..ae50a8042 100644 --- a/docs/rust.md +++ b/docs/rust.md @@ -37,7 +37,7 @@ bin/og the shim tasks call; runs bin/.og/og-linux- (og-da tests/rust/ the differential test against CPython, the image test ``` -Each ported tool becomes a subcommand named after its file in `bin/` without `.py`, in kebab-case. Today `og` has `og qc sex-loss` (step 33), `og version`, `og --tree-hash`, `og map`, `og setup mask-fasta`, `og run state|tasks|eta|headroom|delta`, `og dev parity-compare`, `og pgx cyp2d6-bed|cyp2d6-measure|cyp2d6-check`, `og str motif|catalog`, `og cdna check|scan` (were `clip_check.py` and `cdna_scan.py`; `clip_check.py` keeps only the names `collect_summary.py` imports), `og comphet phase` (was `comphet_phase.py`), `og bam contigs`, `og setup unzip|genetic-map` and two hidden ones the tests drive: `og pycompat-probe` and `og bam-fetch-probe`. The filters of steps 04, 05 and 07 are grouped by step: `og sv svlen-fix` and `og sv inv-filter` (were `sv_svlen_fix.py` and `manta_inv_filter.py`), `og pgx refblock-fill` and `og pgx mito-merge` (were `pgx_refblock_fill.py` and `pgx_mito_merge.py`). Step 05 builds AnnotSV's input with `og sv add-callers` and checks its class 4 and 5 calls with `og sv postfilter` (were `cnvpytor_to_vcf.py` and `annotsv_postfilter.py`); both apply one read-depth rule for long Delly calls (`og-steps::sv`). `cnvpytor_to_vcf.py` keeps only the Python copy of that rule, which step 39's couple screen still imports, and `annotsv_postfilter.py` only the column names and the consensus join the reports import. `og str catalog` (was `strchive_convert.py`) builds the STRchive catalog of steps 09 and 9b; `og setup unzip` and `og setup genetic-map` replace two small Python programs inside `setup.sh`. Step 25 has `og prs format` (a port of `collect_summary.py prs-format`, named by what it works on) and `og prs table` (was `prs-table`). Step 33 has `og qc sample` (was `sample-qc`), and step 20 has `og mito annotate|deletions|copies` (were `mito-annotate`, `mito-deletions` and `mito-copies`). Steps 23 and 31 use `og vcf csq-rank` and `og table constraint-join` (ports of the awk programs `bin/csq_rank.awk` and `bin/constraint_join.awk`, named by what they work on). Step 27 uses `og pgx cpic-report` and `og pgx step36-outside` (were the two commands of `pgx_parse.py`). Step 49 runs `og chip prepare`, `og chip report` and `og chip check-table` (the commands of `chip_report.py`, which keeps only the detection limit and the column names the reports import). Steps 32, 35, 41 and 46 run `og cn bed`, `og cn call`, `og cn intervals`, `og cn check-bam`, `og cn marker-regions` and `og cn marker-read` ([Copy number from depth](cn-loci.md)); they replace `hba_copy_number.py`, `lpa_kiv2.py`, the commands of `copy_number.py` (which keeps the catalog reader the catalog tests import) and the marker reader of `traits_report.py` (which keeps `MARKERS` for the summary). Step 41 runs `og bloodtype report` (was `bloodtype_report.py`), which reads the marker sites through the VCF's index. Step 26b's six commands are `og ancestry-fine panel-map|check-map|order|sample-bed|chroms|report`; `bin/ancestry_fine_report.py` keeps only the names `collect_summary.py` and the docs demo import. Step 40 runs `og consumer build-sites`, `og consumer export` and `og consumer check` (were the three commands of `export_consumer.py`). Step 38 and the Nextflow family workflow run `og family plan`, `prepare`, `check`, `report`, `segments` and `ibd` (were `family_report.py`, `family_segments.py` and `ibd_report.py`). `og family report` reads the joint VCF once, compresses each child's inheritance table on its own thread, and stops with exit 2 on a joint VCF that is not sorted. The opt-in steps read their tool's output with `og y str-bed|str-report` (45), `og viral report` (47), `og exomiser prepare|report` (48), `og mei report|annotsv-input` (50), `og mocha report` (51), `og vntyper report` (52) and `og locityper report` (53); they replaced `ystr_report.py`, `viral_report.py`, `exomiser_report.py`, `mei_report.py`, `mocha_report.py`, `vntyper_report.py` and `locityper_report.py`. Step 37 places the sample on the current YFull YTree with `og y tree-positions` and `og y tree-place` (were the two commands of `y_tree_place.py`), lists its private SNPs with `og y private-snps` (was `y_private_snps.py`) and writes Yleaf's marker positions with `og y yleaf-positions` (was `yleaf_run.py positions`; `yleaf_run.py predict` stays Python, as glue inside Yleaf's own code). Steps 08, 32 and 36 run `og pgx hla-concordance`, `og pgx pypgx-summary` and `og pgx outside-calls` (were the three commands of `pgx_outside_calls.py`, which keeps the CYP2D6 and HLA rules the reports import), and step 35b runs `og paraphase check-bam|table` (were `paraphase_report.py`, which keeps the SMN1 rule beside Parascopy). Both write their two output files together or not at all. `scripts/owngenome-lookup.sh` runs `og lookup query` ([genotype lookup](lookup.md)), which answers every query of the command line in one process; `og lookup resolve` and `og lookup report` are its two halves, with the command lines of `genome_lookup.py`, which keeps only the query reading and the names rule that `results_index.py` imports. `og recheck` replaced `owngenome_recheck.py` ([Newer databases since a run](recheck.md)); the file stays one release as a shim that runs it, for crontabs. `og dev freshness` replaced `scripts/ci/freshness.py`, the monthly freshness issue of `freshness.yml`. Both ask the network through the host's `curl` (the config on stdin, so no token reaches the process list), which keeps a TLS stack out of the static binary and its cross-builds. +Each ported tool becomes a subcommand named after its file in `bin/` without `.py`, in kebab-case. Today `og` has `og qc sex-loss` (step 33), `og version`, `og --tree-hash`, `og map`, `og setup mask-fasta`, `og run state|tasks|eta|headroom|delta`, `og dev parity-compare`, `og pgx cyp2d6-bed|cyp2d6-measure|cyp2d6-check`, `og str motif|catalog`, `og cdna check|scan` (were `clip_check.py` and `cdna_scan.py`; `clip_check.py` keeps only the names `collect_summary.py` imports), `og comphet phase` (was `comphet_phase.py`), `og bam contigs`, `og setup unzip|genetic-map` and two hidden ones the tests drive: `og pycompat-probe` and `og bam-fetch-probe`. The filters of steps 04, 05 and 07 are grouped by step: `og sv svlen-fix` and `og sv inv-filter` (were `sv_svlen_fix.py` and `manta_inv_filter.py`), `og pgx refblock-fill` and `og pgx mito-merge` (were `pgx_refblock_fill.py` and `pgx_mito_merge.py`). Step 05 builds AnnotSV's input with `og sv add-callers` and checks its class 4 and 5 calls with `og sv postfilter` (were `cnvpytor_to_vcf.py` and `annotsv_postfilter.py`); both apply one read-depth rule for long Delly calls (`og-steps::sv`). `cnvpytor_to_vcf.py` keeps only the Python copy of that rule, which step 39's couple screen still imports, and `annotsv_postfilter.py` only the column names and the consensus join the reports import. `og str catalog` (was `strchive_convert.py`) builds the STRchive catalog of steps 09 and 9b; `og setup unzip` and `og setup genetic-map` replace two small Python programs inside `setup.sh`. Step 25 has `og prs format` (a port of `collect_summary.py prs-format`, named by what it works on) and `og prs table` (was `prs-table`). Step 33 has `og qc sample` (was `sample-qc`), and step 20 has `og mito annotate|deletions|copies` (were `mito-annotate`, `mito-deletions` and `mito-copies`). Steps 23 and 31 use `og vcf csq-rank` and `og table constraint-join` (ports of the awk programs `bin/csq_rank.awk` and `bin/constraint_join.awk`, named by what they work on). Step 27 uses `og pgx cpic-report` and `og pgx step36-outside` (were the two commands of `pgx_parse.py`). Step 49 runs `og chip prepare`, `og chip report` and `og chip check-table` (the commands of `chip_report.py`, which keeps only the detection limit and the column names the reports import). Steps 32, 35, 41 and 46 run `og cn bed`, `og cn call`, `og cn intervals`, `og cn check-bam`, `og cn marker-regions` and `og cn marker-read` ([Copy number from depth](cn-loci.md)); they replace `hba_copy_number.py`, `lpa_kiv2.py`, the commands of `copy_number.py` (which keeps the catalog reader the catalog tests import) and the marker reader of `traits_report.py` (which keeps `MARKERS` for the summary). Step 41 runs `og bloodtype report` (was `bloodtype_report.py`), which reads the marker sites through the VCF's index. Step 26b's six commands are `og ancestry-fine panel-map|check-map|order|sample-bed|chroms|report`; `bin/ancestry_fine_report.py` keeps only the names `collect_summary.py` and the docs demo import. Step 40 runs `og consumer build-sites`, `og consumer export` and `og consumer check` (were the three commands of `export_consumer.py`). Step 38 and the Nextflow family workflow run `og family plan`, `prepare`, `check`, `report`, `segments` and `ibd` (were `family_report.py`, `family_segments.py` and `ibd_report.py`). `og family report` reads the joint VCF once, compresses each child's inheritance table on its own thread, and stops with exit 2 on a joint VCF that is not sorted. The opt-in steps read their tool's output with `og y str-bed|str-report` (45), `og viral report` (47), `og exomiser prepare|report` (48), `og mei report|annotsv-input` (50), `og mocha report` (51), `og vntyper report` (52) and `og locityper report` (53); they replaced `ystr_report.py`, `viral_report.py`, `exomiser_report.py`, `mei_report.py`, `mocha_report.py`, `vntyper_report.py` and `locityper_report.py`. Step 37 places the sample on the current YFull YTree with `og y tree-positions` and `og y tree-place` (were the two commands of `y_tree_place.py`), lists its private SNPs with `og y private-snps` (was `y_private_snps.py`) and writes Yleaf's marker positions with `og y yleaf-positions` (was `yleaf_run.py positions`; `yleaf_run.py predict` stays Python, as glue inside Yleaf's own code). Steps 08, 32 and 36 run `og pgx hla-concordance`, `og pgx pypgx-summary` and `og pgx outside-calls` (were the three commands of `pgx_outside_calls.py`, which keeps the CYP2D6 and HLA rules the reports import), and step 35b runs `og paraphase check-bam|table` (were `paraphase_report.py`, which keeps the SMN1 rule beside Parascopy). Both write their two output files together or not at all. `scripts/owngenome-lookup.sh` runs `og lookup query` ([genotype lookup](lookup.md)), which answers every query of the command line in one process; `og lookup resolve` and `og lookup report` are its two halves, with the command lines of `genome_lookup.py`, which keeps only the query reading and the names rule that `results_index.py` imports. `og recheck` replaced `owngenome_recheck.py` ([Newer databases since a run](recheck.md)); the file stays one release as a shim that runs it, for crontabs. `og dev freshness` replaced `scripts/ci/freshness.py`, the monthly freshness issue of `freshness.yml`. `og dev build cn-intervals|callability|sv-rule-regions` rebuild the committed tables `bin/cn_loci.tsv` (its interval rows), `assets/callability_genes.tsv` and `assets/sv_rule_regions.tsv`; they replaced `scripts/ci/build-cn-intervals.py` and the Python inside `build-callability.sh` and `build-sv-rule-regions.sh`. Those scripts keep the downloads (`curl` through `scripts/ci/lib/fetch-pinned.sh`, each file checked against its sha256) and hand `og` the folder, so `og dev build` reads local files only. Both ask the network through the host's `curl` (the config on stdin, so no token reaches the process list), which keeps a TLS stack out of the static binary and its cross-builds. A Python file whose tool moved to `og` but whose names other Python still imports stays as a *remnant*: a `# remnant:` header names its importers and the lane that deletes it, and `tests/parity/remnants.tsv` lists the names it may hold and the files that import it. `og dev remnants` (Guard, job `rust`) fails when the file defines a name the list does not allow, when an allowed name has no user left, when a listed importer no longer imports it, or when an unlisted file starts to. So a remnant only shrinks. `bin/family_report.py` is the first one: the pedigree, kinship and PAR rules that `couple_report.py` and `pedigree_svg.py` import. `bin/genome_lookup.py` keeps the query reading and the summary names rule that `results_index.py` imports. diff --git a/docs/testing.md b/docs/testing.md index 7cd0a5da2..7302f594c 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -102,7 +102,7 @@ The base leg runs every base case whenever any case runs, because every shard st The images the cases use are cached in four GitHub Actions cache entries, one per group. `scripts/ci/e2e-run.sh --image-groups` names each image's group, taken from its repository name without the tag, so a version bump changes the key of one group only; a group whose key missed restores the newest older entry of that group and pulls only what changed. A group keeps its smallest images up to 1500 MB of zstd output, and the larger ones are pulled on every run. Four groups of 1500 MB and the six data entries below stay under the repository's 10 GB cache limit even while main replaces one group: 5.9 GB of images, 1.5 GB more for the group being replaced, and about 2.4 GB of data. Images pinned by digest are not cached, because `docker load` gives them back without the digest, and they are left out of the group's key, so a digest bump misses no entry. Only the base leg saves, and only on main (or in a run by hand with `save_images`, on any branch), after it pulled every image the cases use and got through its cases, and only the groups whose key missed: a cache saved by a pull request can be read by that pull request alone and would only push main's entries out. The job summary lists every image with its size and what the cache did with it: restored on an exact key hit, restored from an older entry, saved, or pulled. Every merge in the e2e scope starts a run on main, so the cache follows main's images. GitHub drops an entry no run has read for 7 days, so after a quiet week the next run on main fills the cache again. Until this layout the cache was one entry, and by October 2026 it had grown to about 11 GB, over the limit, so every run missed it. -Six downloads are cached the same way. Two come from EBI's 1000 Genomes server, which is slow and at times stalls: MoChA's chr20 panel with its manifest (0.64 GB, case `mocha-1`), keyed on `MOCHA_PANEL_MANIFEST_SHA256` and the chromosome, and the 1000 Genomes phase 3 sites VCF (1.36 GB, case `archaic-ancestry-1`), keyed on `ARCHAIC_1KG_SITES_MD5`. On a miss, each of those two is tried first at `E2E_DATA_MIRROR`, the GitHub release `e2e-data-1` that holds copies of the pinned files, and then at EBI; the pin checks both. The release also holds the four PGS Catalog score files, the HAPNEST panel and the table builders' downloads below, so an EBI outage costs none of them on a cold cache. It holds the TRExplorer catalog of `build-sv-rule-regions.sh` too (the v2.1 file its `TRX_SHA256` pins), so a miss tries the mirror first and the catalog's own GitHub release after it. Two are small and failed whole runs on 2026-11-04, because their case is in the base leg or a shard that cannot start without it: `hla.dat` of the two IPD-IMGT/HLA releases `bash-step-08-hla` installs, from raw.githubusercontent.com (which served it under 1 KB/s), with GENCODE's gene and exon files from EBI's FTP server (curl exit 7 and 28), and the four real PGS Catalog score files of `prs-scores-and-run-options-1-strand-sex`, from the same server. `scripts/ci/e2e-run.sh --data-key` names their keys: the releases and a hash of the code in `scripts/lib/common.sh` that installs them, and a hash of the score ids and their server. `scripts/ci/e2e-run.sh` hands each case its files just before it runs, at the paths `setup.sh` and `install_data_file` read them from, so they download nothing; the score files the case copies from the cache itself. On a miss it downloads them first, checked against their pins or the checksum files the servers publish beside them (`MD5SUMS`, `md5checksum.txt`, the `.md5` of each score), trying up to 6 times 30 s apart. The four score files are pinned by md5 in `PGS_MD5` of `e2e-run.sh` (the md5 the PGS Catalog served beside each on 2026-11-06), tried on the mirror first and then at EBI; a restored or downloaded score must match its pin as well as its `.md5`. A file gets at most `E2E_MIRROR_BUDGET` seconds on the mirror (900, or half the fill's budget when that is less), so a mirror that trickles leaves EBI its turn. One fill may take at most `E2E_DATA_BUDGET` seconds (2400 in `e2e.yml`): a source that trickles above `fetch`'s stall limit would otherwise run on to the leg's 240-minute limit, as EBI's 1000 Genomes server did twice on 2026-11-05. A fill still running then is stopped and its case fails as `FAIL (data fill over 2400s)` without running, since its own download would be as slow. A restored entry is checked again before use (`--data-check`): the IPD-IMGT/HLA and GENCODE entry against the sha256 list written when it was filled, the score files against their `.md5`; one that fails is downloaded again. `bash-step-08-hla` then runs with every HTTP and HTTPS request sent to a closed port, after a control shows that a request under that block fails at the closed port, so a warm cache cannot reach either host. Only the leg whose cases include the case restores and saves its entry, and it saves only on main or in a run by hand with `save_images`, as with the images. The job summary says whether each key hit, and the case's log says where its files came from. The fifth is the PGS Catalog's synthetic HAPNEST panel (0.27 GB, cases `prs-2-ancestry` and `prs-scores-and-run-options-2-panel-chrx`) from EBI's FTP server, which failed cases with curl exit 7 on 2026-11-04. It is keyed on the sha256 `tests/e2e/lib.sh` pins (`HAPNEST_SHA256`). The cases handle it themselves (`setup_ancestry_panel` in `lib.sh`): with the panel in the cache, or on a miss from `E2E_DATA_MIRROR`, they check the sha256 and run `setup.sh --ancestry-panel` with every HTTP and HTTPS request sent to a closed port, so an EBI outage cannot fail them; when neither has it, `setup.sh` downloads the panel from EBI as a user's run does. The first case puts a panel that did not come from the cache in the cache for the save step. A red control in `prs-2-ancestry` runs `setup.sh` with the network blocked and no panel, which must fail. The sixth, entry `catalogs`, holds the pinned downloads of `scripts/ci/build-callability.sh`, `scripts/ci/build-cn-intervals.py` and `scripts/ci/build-sv-rule-regions.sh`, whose cases `callability-1-catalog`, `cn-intervals-1-catalog` and `followups-ci-3-sv-rule-regions` rebuild `assets/callability_genes.tsv`, the interval rows of `bin/cn_loci.tsv` and `assets/sv_rule_regions.tsv` and compare them with the committed files. Each builder lists its downloads with `--list`; the key is a hash of those lists, so a new pin is a new key. GENCODE's GTF among them comes from EBI's FTP server, which refused every connection for minutes on 2026-11-06 and failed `callability-1-catalog` on main (run 42259432681). The runner fills the entry as it does the others, mirror first and every file checked against its sha256, and each case points its builder's cache at its folder there; a builder that still has to download tries for about five minutes, 30 s apart. Every download, there and in `setup.sh`, gives up on a transfer that stalls (`FETCH_STALL`, see `fetch` in `scripts/lib/common.sh`), and the next try resumes it. +Six downloads are cached the same way. Two come from EBI's 1000 Genomes server, which is slow and at times stalls: MoChA's chr20 panel with its manifest (0.64 GB, case `mocha-1`), keyed on `MOCHA_PANEL_MANIFEST_SHA256` and the chromosome, and the 1000 Genomes phase 3 sites VCF (1.36 GB, case `archaic-ancestry-1`), keyed on `ARCHAIC_1KG_SITES_MD5`. On a miss, each of those two is tried first at `E2E_DATA_MIRROR`, the GitHub release `e2e-data-1` that holds copies of the pinned files, and then at EBI; the pin checks both. The release also holds the four PGS Catalog score files, the HAPNEST panel and the table builders' downloads below, so an EBI outage costs none of them on a cold cache. It holds the TRExplorer catalog of `build-sv-rule-regions.sh` too (the v2.1 file its `TRX_SHA256` pins), so a miss tries the mirror first and the catalog's own GitHub release after it. Two are small and failed whole runs on 2026-11-04, because their case is in the base leg or a shard that cannot start without it: `hla.dat` of the two IPD-IMGT/HLA releases `bash-step-08-hla` installs, from raw.githubusercontent.com (which served it under 1 KB/s), with GENCODE's gene and exon files from EBI's FTP server (curl exit 7 and 28), and the four real PGS Catalog score files of `prs-scores-and-run-options-1-strand-sex`, from the same server. `scripts/ci/e2e-run.sh --data-key` names their keys: the releases and a hash of the code in `scripts/lib/common.sh` that installs them, and a hash of the score ids and their server. `scripts/ci/e2e-run.sh` hands each case its files just before it runs, at the paths `setup.sh` and `install_data_file` read them from, so they download nothing; the score files the case copies from the cache itself. On a miss it downloads them first, checked against their pins or the checksum files the servers publish beside them (`MD5SUMS`, `md5checksum.txt`, the `.md5` of each score), trying up to 6 times 30 s apart. The four score files are pinned by md5 in `PGS_MD5` of `e2e-run.sh` (the md5 the PGS Catalog served beside each on 2026-11-06), tried on the mirror first and then at EBI; a restored or downloaded score must match its pin as well as its `.md5`. A file gets at most `E2E_MIRROR_BUDGET` seconds on the mirror (900, or half the fill's budget when that is less), so a mirror that trickles leaves EBI its turn. One fill may take at most `E2E_DATA_BUDGET` seconds (2400 in `e2e.yml`): a source that trickles above `fetch`'s stall limit would otherwise run on to the leg's 240-minute limit, as EBI's 1000 Genomes server did twice on 2026-11-05. A fill still running then is stopped and its case fails as `FAIL (data fill over 2400s)` without running, since its own download would be as slow. A restored entry is checked again before use (`--data-check`): the IPD-IMGT/HLA and GENCODE entry against the sha256 list written when it was filled, the score files against their `.md5`; one that fails is downloaded again. `bash-step-08-hla` then runs with every HTTP and HTTPS request sent to a closed port, after a control shows that a request under that block fails at the closed port, so a warm cache cannot reach either host. Only the leg whose cases include the case restores and saves its entry, and it saves only on main or in a run by hand with `save_images`, as with the images. The job summary says whether each key hit, and the case's log says where its files came from. The fifth is the PGS Catalog's synthetic HAPNEST panel (0.27 GB, cases `prs-2-ancestry` and `prs-scores-and-run-options-2-panel-chrx`) from EBI's FTP server, which failed cases with curl exit 7 on 2026-11-04. It is keyed on the sha256 `tests/e2e/lib.sh` pins (`HAPNEST_SHA256`). The cases handle it themselves (`setup_ancestry_panel` in `lib.sh`): with the panel in the cache, or on a miss from `E2E_DATA_MIRROR`, they check the sha256 and run `setup.sh --ancestry-panel` with every HTTP and HTTPS request sent to a closed port, so an EBI outage cannot fail them; when neither has it, `setup.sh` downloads the panel from EBI as a user's run does. The first case puts a panel that did not come from the cache in the cache for the save step. A red control in `prs-2-ancestry` runs `setup.sh` with the network blocked and no panel, which must fail. The sixth, entry `catalogs`, holds the pinned downloads of `scripts/ci/build-callability.sh`, `scripts/ci/build-cn-intervals.sh` and `scripts/ci/build-sv-rule-regions.sh`, whose cases `callability-1-catalog`, `cn-intervals-1-catalog` and `followups-ci-3-sv-rule-regions` rebuild `assets/callability_genes.tsv`, the interval rows of `bin/cn_loci.tsv` and `assets/sv_rule_regions.tsv` and compare them with the committed files. Each builder lists its downloads with `--list`; the key is a hash of those lists, so a new pin is a new key. GENCODE's GTF among them comes from EBI's FTP server, which refused every connection for minutes on 2026-11-06 and failed `callability-1-catalog` on main (run 42259432681). The runner fills the entry as it does the others, mirror first and every file checked against its sha256, and each case points its builder's cache at its folder there; a builder that still has to download tries for about five minutes, 30 s apart, and gives up at once on a 403 or 404 (`fetch_pinned` in `scripts/ci/lib/fetch-pinned.sh`, which the three builders share). The builders download with `curl` and hand the folder to `og dev build`, which reads local files only. Every download, there and in `setup.sh`, gives up on a transfer that stalls (`FETCH_STALL`, see `fetch` in `scripts/lib/common.sh`), and the next try resumes it. How long it took in one job, before the split: a passing e2e job took about 105 to 205 minutes with 84 cases, and 230 minutes of cases with 130 (run 41390434697, 2026-11-02). The runner decides much of it, not the cache: the same CPU model took 106 minutes in one run and 166 in another. In October 2026, with 84 cases, the same code took 106 and 166 minutes on an AMD EPYC 9V45 (runs 39197502832 and 39075176392, the second of which also saved the cache), 112 minutes on an Intel Xeon 6973P (run 39088086725), 156 and 191 minutes on an AMD EPYC 7763 (runs 39075764278 and 39075201465) and 205 minutes on an AMD EPYC 9V74 (run 39197512594, which saved three groups). The image cache saves about 7 minutes of that: pulling every image took 9 minutes, restoring and loading the four groups took 2.5 (run 39075201465 restored 39 images in 118 seconds), and saving them all took 12. Each leg has a limit of 240 minutes (355 for a traced shard leg); the base leg waits up to 75 minutes for a fixture a pull request bumps. Each leg's first steps print the runner's CPU model and CPU count to the log and the job summary, so a later timing can name the hardware. diff --git a/scripts/ci/build-callability.sh b/scripts/ci/build-callability.sh index e47058489..ab594545b 100755 --- a/scripts/ci/build-callability.sh +++ b/scripts/ci/build-callability.sh @@ -8,9 +8,8 @@ # scripts/ci/build-callability.sh --list the downloads, one "NAME URL SHA256" per line # # Downloads seven pinned files (about 85 MB, each checked against its -# sha256) into a temporary folder, or reads them from CALLABILITY_CACHE when -# that folder holds them with the right sha256, writes OUT and deletes the -# downloads: +# sha256) into CALLABILITY_CACHE, where a file with the right sha256 is read +# as it is, or else into a temporary folder that is deleted, and writes OUT: # # GIAB CMRG v1.00's GRCh38 gene coordinates (Wagner 2022, PMID 35132260): # the 273 challenging medically relevant genes @@ -53,7 +52,9 @@ # NEB, PIK3CA, PLG, VWF). Step 35's own report section reads all of them. # tests/test_paralog_loci.py holds that list and goes red when it changes. # -# Needs curl, python3 and sha256sum or shasum. The e2e case +# The rules are in og dev build callability +# (rust/crates/og-dev/src/build/callability.rs). Needs curl, og +# (scripts/setup.sh --og) and sha256sum or shasum. The e2e case # tests/e2e/callability-1-catalog.sh runs it and fails on any difference with # the committed file. set -euo pipefail @@ -81,208 +82,17 @@ if [ "${1:-}" = --list ]; then exit 0 fi -sha256() { - if command -v sha256sum >/dev/null 2>&1; then sha256sum "$1" | cut -d' ' -f1; else shasum -a 256 "$1" | cut -d' ' -f1; fi -} - -WORK=$(mktemp -d) -trap 'rm -rf "$WORK"' EXIT - -# get NAME URL SHA256: WORK/NAME from CALLABILITY_CACHE when it matches, else -# downloaded; stop unless it matches. -get() { - local cached="${CALLABILITY_CACHE:-}/$1" got - if [ -n "${CALLABILITY_CACHE:-}" ] && [ -f "$cached" ] && [ "$(sha256 "$cached")" = "$3" ]; then - cp "$cached" "${WORK}/$1" - return - fi - echo "Downloading $2" >&2 - # Tries spread over about five minutes: EBI's FTP server refused every - # connection for minutes at a time (curl exit 7, 2026-11-06). - curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "${WORK}/$1" "$2" - got=$(sha256 "${WORK}/$1") - if [ "$got" != "$3" ]; then - echo "ERROR: $1 has sha256 ${got}, want $3" >&2 - exit 1 - fi - if [ -n "${CALLABILITY_CACHE:-}" ]; then - mkdir -p "$CALLABILITY_CACHE" && cp "${WORK}/$1" "$cached" - fi -} +# shellcheck source=scripts/ci/lib/fetch-pinned.sh +. "${REPO}/scripts/ci/lib/fetch-pinned.sh" +DIR=${CALLABILITY_CACHE:-} +if [ -z "$DIR" ]; then + DIR=$(mktemp -d) + trap 'rm -rf "$DIR"' EXIT +fi for f in "${FILES[@]}"; do IFS='|' read -r name url sum <<<"$f" - get "$name" "$url" "$sum" + fetch_pinned "$DIR" "$name" "$url" "$sum" || exit 1 done -python3 - "$WORK" "$REPO" "${OUT}.part" <<'PY' -import csv -import gzip -import re -import sys -from collections import defaultdict - -work, repo, out = sys.argv[1:4] -CHROMS = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY"] -ORDER = {c: i for i, c in enumerate(CHROMS)} -PRIMARY = set(CHROMS) | {"chrM"} -ATT = re.compile(r'(\S+) "([^"]*)"') -SEGDUP_SHARE = 0.5 -REASONS = ("CMRG", "segdup", "false_duplication", "collapsed", "ALT-only") - - -def bed(name): - d = defaultdict(list) - with (gzip.open if name.endswith(".gz") else open)(f"{work}/{name}", "rt") as f: - for line in f: - c = line.rstrip("\n").split("\t") - if len(c) >= 3 and not line.startswith(("#", "track", "browser")): - d[c[0]].append((int(c[1]), int(c[2]))) - for k in d: - d[k].sort() - return d - - -def covered(chrom, s, e, d): - """Bases of [s, e) inside the intervals of d (merged as it goes).""" - n, cur = 0, s - for a, b in d.get(chrom, []): - if b <= cur: - continue - if a >= e: - break - a, b = max(a, cur), min(b, e) - if b > a: - n += b - a - cur = b - return n - - -# ALT scaffolds: GenBank accession -> UCSC name (GENCODE names them by accession) -alt = {} -with open(f"{work}/report.txt") as f: - for line in f: - c = line.rstrip("\n").split("\t") - if not line.startswith("#") and len(c) >= 10 and c[1] == "alt-scaffold": - alt[c[4]] = c[9] - -genes = {} # hgnc id -> [symbol, chrom, start0, end] (any gene type) -coding = set() # hgnc ids of the protein-coding ones -on_primary = set() # every gene symbol on chr1-22, X, Y or M -alt_genes = {} # hgnc id -> [symbol, UCSC ALT name] (first ALT copy) -with gzip.open(f"{work}/gencode.gtf.gz", "rt") as f: - for line in f: - if line.startswith("#"): - continue - c = line.split("\t", 9) - if c[2] != "gene": - continue - a = dict(ATT.findall(c[8])) - if c[0] in PRIMARY: - on_primary.add(a.get("gene_name")) - h = a.get("hgnc_id") - if not h or a["gene_id"].endswith("_PAR_Y"): - continue - if c[0] in ORDER and h not in genes: - genes[h] = [a["gene_name"], c[0], int(c[3]) - 1, int(c[4])] - if a.get("gene_type") == "protein_coding": - coding.add(h) - elif c[0] in alt and h not in alt_genes and a.get("gene_type") == "protein_coding": - alt_genes[h] = [a["gene_name"], alt[c[0]]] - -reasons = defaultdict(set) -by_symbol = {g[0]: h for h, g in genes.items()} -with open(f"{work}/cmrg.bed") as f: - for line in f: - chrom, s, e, sym = line.rstrip("\n").split("\t")[:4] - s, e = int(s), int(e) - h = by_symbol.get(sym) - if h is None: - # renamed since CMRG v1.00: the protein-coding gene that covers most of its span - best = max(((min(e, genes[k][3]) - max(s, genes[k][2]), k) for k in coding if genes[k][1] == chrom), - default=(0, None)) - if best[0] <= 0: - sys.exit(f"ERROR: CMRG gene {sym} ({chrom}:{s}-{e}) matches no GENCODE 50 gene") - h = best[1] - reasons[h].add("CMRG") - -segdups, collapsed = bed("segdups.bed.gz"), bed("collapsed.bed.gz") -falsedup = bed("falsedup_correct.bed.gz") -for k, v in bed("falsedup_incorrect.bed.gz").items(): - falsedup[k] = sorted(falsedup[k] + v) -for h in coding: - sym, chrom, s, e = genes[h] - if covered(chrom, s, e, segdups) >= SEGDUP_SHARE * (e - s): - reasons[h].add("segdup") - if covered(chrom, s, e, falsedup): - reasons[h].add("false_duplication") - if covered(chrom, s, e, collapsed): - reasons[h].add("collapsed") -for h, (sym, name) in alt_genes.items(): - if h not in genes and sym not in on_primary: - reasons[h].add("ALT-only") - - -def table(path): - with open(path) as f: - return list(csv.DictReader((x for x in f if x.strip() and not x.startswith("#")), delimiter="\t", - quoting=csv.QUOTE_NONE)) - - -parascopy = {g for r in table(f"{repo}/assets/paralog_loci.tsv") for g in r["genes"].split(",")} -# A gene step 32 skips on the default no-ALT reference (callable_without_alt no) has no -# caller there, and step 32 refuses a BAM with ALT contigs, so it is not pypgx's. -pypgx = {r["gene"] for r in table(f"{repo}/assets/pgx_genes.tsv") - if r["mode"] == "bam" and r["callable_without_alt"] != "no"} -# A cn_loci gene row reports the genes it names (its locus, as RHD, or LPA for LPA_KIV2, -# or a target, as HBA1 and HBA2), not every gene whose span holds one of its targets. -cn_names = {} -for r in table(f"{repo}/bin/cn_loci.tsv"): - if r["kind"] == "gene": - names = {r["locus"], r["locus"].split("_")[0]} - for part in r["targets"].split(";"): - names |= {part.split(":")[0], part.split(":")[0].split("_")[0]} - for n in names: - cn_names.setdefault(n, r["locus"]) - - -def owner(sym): - if sym in parascopy: - return "parascopy" - if sym in pypgx: - return "pypgx" - if sym in cn_names: - return f"cn_loci:{cn_names[sym]}" - return "." - - -rows = [] -for h, why in reasons.items(): - if h in genes: - sym, chrom, s, e = genes[h] - rows.append((ORDER[chrom], s, sym, h, chrom, s, e, why, owner(sym))) - else: - sym, name = alt_genes[h] - rows.append((len(ORDER), 0, sym, h, name, ".", ".", why, owner(sym))) -rows.sort() -syms = [r[2] for r in rows] -dup = sorted({x for x in syms if syms.count(x) > 1}) -if dup: - sys.exit(f"ERROR: gene symbols listed twice: {', '.join(dup)}") -counts = {k: sum(1 for r in rows if k in r[7]) for k in REASONS} -owned = sum(1 for r in rows if r[8] != ".") -with open(out, "w") as f: - f.write("# assets/callability_genes.tsv: the genes where short reads are often not placed with confidence, so a\n" - "# 'no finding' there needs a reading of its own. Step 16b measures each gene on chr1-22, X and Y; the\n" - "# report lists the ones it could not read. Written by scripts/ci/build-callability.sh, which says how the\n" - "# rows and the owners are chosen; rebuild it there, do not edit it by hand. start is 0-based and end\n" - "# 1-based, as in a BED file; an ALT-only gene has its ALT scaffold as chrom and '.' as start and end.\n" - f"# Genes: {len(rows)} (" + ", ".join(f"{k} {counts[k]}" for k in REASONS) + f"); with an owner: {owned}.\n") - f.write("gene\thgnc_id\tchrom\tstart\tend\treasons\towner\n") - for r in rows: - f.write("\t".join(map(str, (r[2], r[3], r[4], r[5], r[6], ",".join(k for k in REASONS if k in r[7]), r[8]))) - + "\n") -print(f"Genes: {len(rows)} (" + ", ".join(f"{k} {counts[k]}" for k in REASONS) + f"); with an owner: {owned}", - file=sys.stderr) -PY -mv "${OUT}.part" "$OUT" +"${REPO}/bin/og" dev build callability --work "$DIR" --repo "$REPO" --out "$OUT" echo "Wrote ${OUT}" >&2 diff --git a/scripts/ci/build-cn-intervals.py b/scripts/ci/build-cn-intervals.py deleted file mode 100755 index 5982debdf..000000000 --- a/scripts/ci/build-cn-intervals.py +++ /dev/null @@ -1,423 +0,0 @@ -#!/usr/bin/env python3 -"""build-cn-intervals.py: write the interval rows of bin/cn_loci.tsv. - -The interval rows are the copy-number checks step 35 runs with one -`samtools bedcov` pass (og cn intervals): - - - every region row of ClinGen Dosage Sensitivity's recurrent CNV BED, - version 2.1, GRCh38 (CLINGEN_RECURRENT_CNV_URL and _SHA256 in - versions.env, CC0): 22q11.2, 1q21.1, 15q11.2, 16p11.2, 17p12, 17p11.2, - 7q11.23, 3q29, Xp22.31 and the rest. The file's breakpoint and LCR rows - (itemRgb 0,0,0) are the segmental duplications that bound the regions: - they are cut out of the targets and kept out of the controls. Its chrX - rows are spelled 'chrx' and are normalised to chrX; any other name that is - not a GRCh38 chromosome stops the build; - - two ClinGen curations the BED lacks (HAND_ROWS): the 1p36 terminal - region (ISCA-37434, haploinsufficiency 3) and the Xq28 region that holds - MECP2 (ISCA-46304, triplosensitivity 3), from - ClinGen_region_curation_list_GRCh38.tsv as read on 2026-11-04; - - three carrier deletion alleles at their ClinVar GRCh38 spans (HAND_ROWS): - del(GJB6-D13S1830) (VCV 5546), CFTRdele2,3 (VCV 66105) and a 621 bp HBB - deletion that removes exon 3 (VCV 3769040); and two beta-globin - deletions at their HbVar GRCh38 spans, checked on 2026-11-04 by placing - 50 bp on each side of the NG_000007.3 breakpoints on NC_000011.10: the - 619 bp HBB deletion common in India and Pakistan (HbVar 979, - NG_000007.3:g.71609_72227del619, Orkin 1979) and the 13.4 kb Sicilian - (delta-beta)0 deletion from HBD to 3' of HBB (HbVar 1035, - NG_000007.3:g.64336_77738del13403, Henthorn 1990). The alpha-globin - --SEA deletion is read by the HBA row (step 35's HBA table), not here; - - the Y chromosome AZF regions (HAND_ROWS), read only on a sample with a Y: - AZFa as USP9Y and DDX3Y, AZFb as KDM5D and EIF1AY (single-copy genes each - deletion removes), and AZFc as the four DAZ copies added up (a complete - AZFc deletion leaves none, the partial gr/gr and b2/b3 deletions two; - which copies are lost, and the breakpoints, are not claimed). Gene spans - from Ensembl's REST lookup, checked on 2026-11-04. Controls: ZFY and - NLGN4Y, single-copy genes outside every AZF region. - -Targets are cut to pieces of at most 500 kb and lose every stretch inside a -breakpoint or LCR row of the BED, a GIAB v3.6 segmental duplication, an -assembly gap (GIAB's gaps with 15 kb on each side), a satellite or a region of -a model the pinned Parascopy data holds (assets/parascopy_models_v1.7.txt: -Parascopy owns those, CHRNA7 in 15q13.3 among them); the AZF targets keep -their segmental duplications, since the DAZ copies are one. Controls are two -100 kb windows, one on each side of the cluster of ClinGen rows the region -belongs to, the nearest that stay clear of all of that, within 3 Mb; a side -with none (a chromosome end) is left out. Every row: depth pooled over the -targets against depth pooled over the controls, scale 2, a whole number -within 0.3 (round:0.3), at least 8x of control depth (4x on chrY), the two -flanks within 1.5x of each other; cross_check cnvpytor:cnvs (step 18). - -Usage: - scripts/ci/build-cn-intervals.py [--catalog bin/cn_loci.tsv] [--cache DIR] [--check] - scripts/ci/build-cn-intervals.py --list the downloads, one "NAME URL SHA256" per line - -Rewrites the catalog's rows from the marker line on (MARKER), or with --check -exits 1 when the catalog is not what the downloads give. The downloads (about -150 KB) are checked against their sha256 and kept in --cache when given. The -e2e case tests/e2e/cn-intervals-1-catalog.sh runs --check. -""" -import argparse -import gzip -import hashlib -import os -import re -import sys -import tempfile -import time -import urllib.request - -REPO = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) -MARKER = "# --- interval rows: written by scripts/ci/build-cn-intervals.py; do not edit by hand ---" -GIAB = "https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/genome-stratifications/v3.6/GRCh38@all" -MASKS = [ - ("segdups.bed.gz", f"{GIAB}/SegmentalDuplications/GRCh38_segdups.bed.gz", - "09ac029a7d0d681a9f0ed6581a6bc10fc708a834bc3f91c46a9bef504f6055c5"), - ("gaps.bed.gz", f"{GIAB}/OtherDifficult/GRCh38_gaps_slop15kb.bed.gz", - "b5fdbfc7e5b9cf4fa356f23841fd7b8a74fa239fd8cbc5b5c3463c9e7bdc0902"), - ("satellites.bed.gz", f"{GIAB}/LowComplexity/GRCh38_satellites_slop5.bed.gz", - "10207aa2dac0e1189ea553a374abbdcb4de25a81378d88dc4ba03aea8718b5a2"), -] -REGION_RGB, MASK_RGB = "255,111,0", "0,0,0" -CHROMS = {f"chr{c}" for c in list(range(1, 23)) + ["X", "Y"]} -# GRCh38 chromosome lengths (the no-ALT analysis set's .fai). -LENGTHS = {"chr1": 248956422, "chr2": 242193529, "chr3": 198295559, "chr4": 190214555, "chr5": 181538259, - "chr6": 170805979, "chr7": 159345973, "chr8": 145138636, "chr9": 138394717, "chr10": 133797422, - "chr11": 135086622, "chr12": 133275309, "chr13": 114364328, "chr14": 107043718, "chr15": 101991189, - "chr16": 90338345, "chr17": 83257441, "chr18": 80373285, "chr19": 58617616, "chr20": 64444167, - "chr21": 46709983, "chr22": 50818468, "chrX": 156040895, "chrY": 57227415} -PIECE, FLANK, GAP, SEARCH, STEP, MIN_PIECE = 500_000, 100_000, 50_000, 3_000_000, 10_000, 500 - -# Rows the BED lacks: locus, chrom, [(name, 0-based start, end)] targets, [controls] or None (flanks), -# target_depth, whether segmental duplications are cut out of the targets, pmids. The label of the -# first target is the row's name in step 35's table. -Y_CONTROLS = [("ZFY", 2935200, 2982579, "ZFY, single-copy control on Yp"), - ("NLGN4Y", 14522572, 14845654, "NLGN4Y, single-copy control on Yq")] -HAND_ROWS = [ - ("1p36_terminal", "chr1", [("1p36 terminal region (includes GABRD), ClinGen ISCA-37434", 898702, 6229913)], - None, "pooled", True, "."), - ("Xq28_MECP2", "chrX", [("Xq28 region (includes MECP2), ClinGen ISCA-46304: the MECP2 duplication", - 154008528, 154110279)], None, "pooled", True, "."), - ("GJB6_D13S1830", "chr13", [("del(GJB6-D13S1830), ClinVar 5546", 20223037, 20531806)], None, "pooled", True, - "11807148"), - ("CFTR_dele2_3", "chr7", [("CFTRdele2,3 (21 kb), ClinVar 66105", 117498312, 117519392)], None, "pooled", True, - "."), - ("HBB_exon3_del", "chr11", [("HBB deletion of exon 3 (621 bp), ClinVar 3769040", 5225254, 5225875)], None, - "pooled", True, "."), - ("HBB_619bp_del", "chr11", [("HBB 619 bp deletion (South Asian), HbVar 979", 5225388, 5226007)], None, - "pooled", True, "287080"), - ("HBD_HBB_Sicilian_del", "chr11", [("Sicilian (delta-beta)0 deletion (13.4 kb, HBD to HBB 3'), HbVar 1035", - 5219877, 5233280)], None, "pooled", True, "2307466"), - ("AZFa", "chrY", [("AZFa (USP9Y and DDX3Y)", 12537649, 12860852), ("DDX3Y", 12904107, 12920478)], - Y_CONTROLS, "pooled", False, "37674303"), - ("AZFb", "chrY", [("AZFb (KDM5D and EIF1AY)", 19702297, 19745148), ("EIF1AY", 20575693, 20593154)], - Y_CONTROLS, "pooled", False, "37674303"), - ("AZFc_DAZ", "chrY", [("AZFc (DAZ1 to DAZ4 added up; gr/gr and b2/b3 leave two)", 23129354, 23199010), - ("DAZ2", 23219446, 23291356), ("DAZ3", 24763068, 24813492), ("DAZ4", 24833842, 24907040)], - Y_CONTROLS, "sum", False, "37674303"), -] -COLUMNS = ["locus", "kind", "chrom", "targets", "controls", "details", "target_depth", "control_depth", - "scale", "offset", "bedcov", "mapq0_min", "mapq0_max", "min_control_depth", "max_control_ratio", - "call", "interval", "readings", "input", "owner", "cross_check", "pmids"] - - -class BuildError(ValueError): - pass - - -def sha256(path): - h = hashlib.sha256() - with open(path, "rb") as f: - for b in iter(lambda: f.read(1 << 20), b""): - h.update(b) - return h.hexdigest() - - -def fetch(name, url, want, folder): - path = os.path.join(folder, name) - if not (os.path.exists(path) and sha256(path) == want): - for i in range(3): - try: - with urllib.request.urlopen(url, timeout=120) as r, open(path + ".part", "wb") as f: - f.write(r.read()) - os.replace(path + ".part", path) - break - except OSError as e: - # A missing or forbidden file stays so: no point waiting. - if i == 2 or getattr(e, "code", None) in (403, 404): - raise BuildError(f"cannot download {url}: {e}") - time.sleep(30) - got = sha256(path) - if got != want: - raise BuildError(f"{url}: sha256 {got}, pinned {want}") - return path - - -def versions_env(key): - with open(os.path.join(REPO, "versions.env")) as f: - for line in f: - m = re.match(rf'^{key}="([^"]*)"', line) - if m: - return m.group(1) - raise BuildError(f"no {key} in versions.env") - - -def norm_chrom(c): - """chrx -> chrX (the ClinGen BED's spelling); anything not a GRCh38 chromosome stops the build.""" - c2 = re.sub(r"^chr([xy])$", lambda m: "chr" + m.group(1).upper(), c) - if c2 not in CHROMS: - raise BuildError(f"chromosome {c!r} is not a GRCh38 chromosome name") - return c2 - - -def read_clingen(lines): - """(regions [(chrom, start, end, name)], masks [(chrom, start, end)]) of the BED's lines.""" - regions, masks = [], [] - for line in lines: - if not line.strip() or line.startswith(("#", "track", "browser")): - continue - c = line.rstrip("\n").split("\t") - if len(c) < 9: - raise BuildError(f"ClinGen BED: a row with {len(c)} columns: {line.strip()!r}") - chrom, s, e, name, rgb = norm_chrom(c[0]), int(c[1]), int(c[2]), c[3].strip(), c[8].strip() - if e <= s: - raise BuildError(f"ClinGen BED: {name}: empty span") - if rgb == REGION_RGB: - regions.append((chrom, s, e, name)) - elif rgb == MASK_RGB: - masks.append((chrom, s, e)) - else: - raise BuildError(f"ClinGen BED: {name}: itemRgb {rgb!r} is neither a region nor a breakpoint row") - if not regions: - raise BuildError("ClinGen BED: no region row") - return regions, masks - - -def read_bed(path): - out = [] - with gzip.open(path, "rt") as f: - for line in f: - if line.startswith("#") or not line.strip(): - continue - c = line.split("\t") - out.append((c[0], int(c[1]), int(c[2]))) - return out - - -def read_parascopy(path): - out = [] - with open(path) as f: - for line in f: - if line.startswith("#") or not line.strip(): - continue - _, reg, _ = line.split() - c, span = reg.split(":") - s, e = span.split("-") - out.append((c, int(s) - 1, int(e))) - return out - - -def merge(spans): - out = [] - for s, e in sorted(spans): - if out and s <= out[-1][1]: - out[-1][1] = max(out[-1][1], e) - else: - out.append([s, e]) - return [tuple(x) for x in out] - - -def subtract(s, e, cuts): - """[(s, e)] of s-e outside the merged spans cuts.""" - out, pos = [], s - for cs, ce in cuts: - if ce <= pos or cs >= e: - continue - if cs > pos: - out.append((pos, cs)) - pos = max(pos, ce) - if pos < e: - out.append((pos, e)) - return out - - -def pieces(spans): - out = [] - for s, e in spans: - if e - s < MIN_PIECE: - continue - n = -(-(e - s) // PIECE) - step = -(-(e - s) // n) - for i in range(n): - out.append((s + i * step, min(e, s + (i + 1) * step))) - return out - - -def clear(s, e, cuts): - return all(ce <= s or cs >= e for cs, ce in cuts) - - -def flank(chrom, edge, side, cuts): - """The nearest FLANK window GAP beyond edge on side (-1 left, +1 right), clear of cuts, within SEARCH.""" - for k in range(0, SEARCH // STEP + 1): - if side < 0: - e = edge - GAP - k * STEP - s = e - FLANK - else: - s = edge + GAP + k * STEP - e = s + FLANK - if s < 0 or e > LENGTHS[chrom]: - return None - if clear(s, e, cuts): - return (s, e) - return None - - -def locus_id(name): - band = name.split()[0] - groups = [g for g in re.findall(r"\(([^()]*)\)", name) if not g.strip().startswith("includes")] - words = re.sub(r"[^A-Za-z0-9.-]+", "_", " ".join(groups)).strip("_") - return f"{band}_{words}" if words else band - - -def label(text): - return re.sub(r"[;\t\n]", ",", text).strip() - - -def build(clingen_lines, segdups, gaps, sats, parascopy): - regions, masks = read_clingen(clingen_lines) - by_chrom = {} - for name, items in (("mask", masks), ("segdup", segdups), ("gap", gaps), ("sat", sats), ("para", parascopy)): - for c, s, e in items: - by_chrom.setdefault(c, {}).setdefault(name, []).append((s, e)) - # Clusters: ClinGen rows (regions and breakpoints) that touch or overlap, per chromosome. - clusters = {} - for c in CHROMS: - spans = [(s, e + 1) for cc, s, e, _ in regions if cc == c] + [(s, e + 1) for cc, s, e in masks if cc == c] - clusters[c] = merge(spans) - - def cuts(c, with_segdups=True): - d = by_chrom.get(c, {}) - keys = ("mask", "segdup", "gap", "sat", "para") if with_segdups else ("mask", "gap", "para") - return merge([x for k in keys for x in d.get(k, [])]) - - def flank_cuts(c): - # A control never sits on a target: the ClinGen regions and the HAND_ROWS targets are cut out too. - d = by_chrom.get(c, {}) - return merge([x for k in ("mask", "segdup", "gap", "sat", "para") for x in d.get(k, [])] - + [(s, e) for cc, s, e, _ in regions if cc == c] - + [(s, e) for _, cc, tl, *_ in HAND_ROWS if cc == c for _, s, e in tl]) - - rows, seen = [], {} - - def add(locus, chrom, targets, controls, target_depth, pmids): - if not targets: - raise BuildError(f"{locus}: no target stretch is left outside the masks") - if locus in seen: - seen[locus] += 1 - locus = f"{locus}_{seen[locus]}" - else: - seen[locus] = 1 - y = chrom == "chrY" - tg = ";".join(f"{locus}_t{i + 1}:{s}-{e}" + (f":{lab}" if lab else "") for i, (s, e, lab) in enumerate(targets)) - ct = ";".join(f"{locus}_{n}:{s}-{e}:{lab}" for n, s, e, lab in controls) - rows.append({"locus": locus, "kind": "interval", "chrom": chrom, "targets": tg, "controls": ct, - "details": ".", "target_depth": target_depth, "control_depth": "pooled" if y else "flanks", - "scale": "2", "offset": "0", "bedcov": ".", "mapq0_min": ".", "mapq0_max": ".", - "min_control_depth": "4" if y else "8", "max_control_ratio": "." if y else "1.5", - "call": "round:0.3", "interval": ".", "readings": ".", "input": "bedcov", "owner": "cn-depth", - "cross_check": "cnvpytor:cnvs", "pmids": pmids}) - - def flanks(chrom, s, e): - cs, ce = next(((a, b) for a, b in clusters[chrom] if a <= s and e <= b), (s, e)) - fc = flank_cuts(chrom) - out = [] - for side, edge, n in ((-1, cs, "L"), (1, ce, "R")): - w = flank(chrom, edge, side, fc) - if w: - out.append((n, w[0], w[1], "flank " + ("left" if side < 0 else "right") + " of the region")) - if not out: - raise BuildError(f"{chrom}:{s}-{e}: no flank clear of the masks within {SEARCH // 1_000_000} Mb") - return out - - for chrom, s, e, name in sorted(regions, key=lambda r: (sorted(CHROMS, key=_order).index(r[0]), r[1], r[2])): - segs = pieces(subtract(s, e, cuts(chrom))) - targets = [(a, b, label(name) if i == 0 else "") for i, (a, b) in enumerate(segs)] - add(locus_id(name), chrom, targets, flanks(chrom, s, e), "pooled", ".") - for locus, chrom, tlist, controls, td, with_sd, pmids in HAND_ROWS: - c = cuts(chrom, with_sd) - targets = [] - for name, s, e in tlist: - segs = pieces(subtract(s, e, c)) if td == "pooled" else [(s, e)] - if td == "sum" and not clear(s, e, merge(by_chrom.get(chrom, {}).get("para", []))): - raise BuildError(f"{locus}: {name} overlaps a Parascopy model") - targets += [(a, b, label(name) if i == 0 else "") for i, (a, b) in enumerate(segs)] - if controls is None: - ctl = flanks(chrom, tlist[0][1], tlist[-1][2]) - else: - ctl = [(f"c{i + 1}", s, e, lab) for i, (_, s, e, lab) in enumerate(controls)] - add(locus, chrom, targets, ctl, td, pmids) - return rows - - -def _order(c): - k = c[3:] - return (0, int(k)) if k.isdigit() else (1, k) - - -def render(rows): - return [MARKER + "\n"] + ["\t".join(r[c] for c in COLUMNS) + "\n" for r in rows] - - -def main(argv=None): - ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) - ap.add_argument("--catalog", default=os.path.join(REPO, "bin", "cn_loci.tsv")) - ap.add_argument("--cache") - ap.add_argument("--check", action="store_true") - ap.add_argument("--list", action="store_true") - a = ap.parse_args(argv) - if a.list and a.check: - ap.error("--list and --check do not go together") - if a.list: - try: - files = [("clingen_recurrent_cnv.bed", versions_env("CLINGEN_RECURRENT_CNV_URL"), - versions_env("CLINGEN_RECURRENT_CNV_SHA256"))] + MASKS - except BuildError as e: - print(f"ERROR: {e}", file=sys.stderr) - return 2 - for n, u, h in files: - print(n, u, h) - return 0 - try: - with tempfile.TemporaryDirectory() as tmp: - folder = a.cache or tmp - os.makedirs(folder, exist_ok=True) - bed = fetch("clingen_recurrent_cnv.bed", versions_env("CLINGEN_RECURRENT_CNV_URL"), - versions_env("CLINGEN_RECURRENT_CNV_SHA256"), folder) - sd, gp, st = (read_bed(fetch(n, u, h, folder)) for n, u, h in MASKS) - with open(bed, errors="replace") as f: - rows = build(f.readlines(), sd, gp, st, - read_parascopy(os.path.join(REPO, "assets", "parascopy_models_v1.7.txt"))) - except (OSError, BuildError) as e: - print(f"ERROR: {e}", file=sys.stderr) - return 2 - with open(a.catalog) as f: - lines = f.readlines() - keep = lines[:lines.index(MARKER + "\n")] if MARKER + "\n" in lines else lines - new = keep + render(rows) - if a.check: - if new != lines: - old = set(lines) - set(new) - print(f"{a.catalog} is not what the pinned downloads give ({len(old)} line(s) differ); " - "rerun scripts/ci/build-cn-intervals.py", file=sys.stderr) - return 1 - print(f"{a.catalog}: the {len(rows)} interval rows are what the pinned downloads give") - return 0 - with open(a.catalog + ".tmp", "w") as f: - f.writelines(new) - os.replace(a.catalog + ".tmp", a.catalog) - print(f"{a.catalog}: {len(rows)} interval rows written") - return 0 - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/scripts/ci/build-cn-intervals.sh b/scripts/ci/build-cn-intervals.sh new file mode 100755 index 000000000..ec40d052c --- /dev/null +++ b/scripts/ci/build-cn-intervals.sh @@ -0,0 +1,97 @@ +#!/usr/bin/env bash +# build-cn-intervals.sh: rewrite the interval rows of bin/cn_loci.tsv, the +# copy-number checks step 35 runs with one samtools bedcov pass +# (og cn intervals). +# +# Usage: scripts/ci/build-cn-intervals.sh [--catalog FILE] [--cache DIR] [--check] +# scripts/ci/build-cn-intervals.sh --list the downloads, one "NAME URL SHA256" per line +# +# Rewrites FILE (default bin/cn_loci.tsv) from its marker line on, or with +# --check exits 1 when it is not what the downloads give. The four downloads +# (about 150 KB, each checked against its sha256) are kept in --cache when +# given, else in a temporary folder that is deleted: +# +# ClinGen Dosage Sensitivity's recurrent CNV BED, version 2.1, GRCh38 +# (CLINGEN_RECURRENT_CNV_URL and _SHA256 in versions.env, CC0) +# GIAB genome stratifications v3.6, GRCh38: segmental duplications, +# assembly gaps with 15 kb on each side, satellites +# +# with assets/parascopy_models_v1.7.txt. Every region row of the BED is a +# row, then the rows the BED lacks: 1p36 and Xq28 (MECP2) from ClinGen's +# curations, five carrier deletions (GJB6-D13S1830, CFTRdele2,3, three +# beta-globin deletions) and the Y chromosome's AZF regions. Targets are cut +# to pieces of at most 500 kb and lose the BED's breakpoint and LCR rows, +# segmental duplications (not on the AZF rows), gaps, satellites and +# Parascopy models; controls are two 100 kb flanks clear of all of that. The +# rules and the hand rows are in og dev build cn-intervals +# (rust/crates/og-dev/src/build/cn_intervals.rs); docs/cn-loci.md says what +# each row reads. +# +# Needs curl, og (scripts/setup.sh --og) and sha256sum or shasum. The e2e +# case tests/e2e/cn-intervals-1-catalog.sh runs --check. +set -euo pipefail + +REPO=$(cd "$(dirname "$0")/../.." && pwd) +# shellcheck source=scripts/ci/lib/fetch-pinned.sh +. "${REPO}/scripts/ci/lib/fetch-pinned.sh" + +CATALOG=${REPO}/bin/cn_loci.tsv +CACHE="" +CHECK="" +LIST="" +while [ $# -gt 0 ]; do + case $1 in + --catalog) CATALOG=${2:?--catalog needs a file}; shift 2 ;; + --cache) CACHE=${2:?--cache needs a folder}; shift 2 ;; + --check) CHECK=1; shift ;; + --list) LIST=1; shift ;; + *) echo "Usage: $0 [--catalog FILE] [--cache DIR] [--check] | --list" >&2; exit 2 ;; + esac +done +if [ -n "$LIST" ] && [ -n "$CHECK" ]; then + echo "ERROR: --list and --check do not go together" >&2 + exit 2 +fi + +# versions_env KEY: the value of KEY in versions.env. +versions_env() { + local v + v=$(sed -n "s/^$1=\"\([^\"]*\)\".*/\1/p" "${REPO}/versions.env" | head -n 1) + [ -n "$v" ] || { echo "ERROR: no $1 in versions.env" >&2; return 1; } + echo "$v" +} + +GIAB=https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/genome-stratifications/v3.6/GRCh38@all +clingen_url=$(versions_env CLINGEN_RECURRENT_CNV_URL) || exit 2 +clingen_sha=$(versions_env CLINGEN_RECURRENT_CNV_SHA256) || exit 2 +FILES=( + "clingen_recurrent_cnv.bed|${clingen_url}|${clingen_sha}" + "segdups.bed.gz|${GIAB}/SegmentalDuplications/GRCh38_segdups.bed.gz|09ac029a7d0d681a9f0ed6581a6bc10fc708a834bc3f91c46a9bef504f6055c5" + "gaps.bed.gz|${GIAB}/OtherDifficult/GRCh38_gaps_slop15kb.bed.gz|b5fdbfc7e5b9cf4fa356f23841fd7b8a74fa239fd8cbc5b5c3463c9e7bdc0902" + "satellites.bed.gz|${GIAB}/LowComplexity/GRCh38_satellites_slop5.bed.gz|10207aa2dac0e1189ea553a374abbdcb4de25a81378d88dc4ba03aea8718b5a2" +) + +if [ -n "$LIST" ]; then + for f in "${FILES[@]}"; do + IFS='|' read -r name url sum <<< "$f" + echo "${name} ${url} ${sum}" + done + exit 0 +fi + +if [ -z "$CACHE" ]; then + CACHE=$(mktemp -d) + trap 'rm -rf "$CACHE"' EXIT +fi +for f in "${FILES[@]}"; do + IFS='|' read -r name url sum <<< "$f" + fetch_pinned "$CACHE" "$name" "$url" "$sum" || exit 2 +done + +# Not exec: the EXIT trap deletes a temporary download folder. +rc=0 +"${REPO}/bin/og" dev build cn-intervals \ + --clingen "${CACHE}/clingen_recurrent_cnv.bed" --segdups "${CACHE}/segdups.bed.gz" \ + --gaps "${CACHE}/gaps.bed.gz" --satellites "${CACHE}/satellites.bed.gz" \ + --parascopy "${REPO}/assets/parascopy_models_v1.7.txt" --catalog "$CATALOG" ${CHECK:+--check} || rc=$? +exit "$rc" diff --git a/scripts/ci/build-sv-rule-regions.sh b/scripts/ci/build-sv-rule-regions.sh index cb21651da..031b80fab 100755 --- a/scripts/ci/build-sv-rule-regions.sh +++ b/scripts/ci/build-sv-rule-regions.sh @@ -10,9 +10,8 @@ # "NAME URL SHA256" per line # # Downloads four pinned files (about 145 MB, each checked against its sha256) -# into a temporary folder, or reads them from SV_RULE_CACHE when that folder -# holds them with the right sha256 (and keeps what it downloads there), writes -# OUT and deletes the downloads: +# into SV_RULE_CACHE, where a file with the right sha256 is read as it is, or +# else into a temporary folder that is deleted, and writes OUT: # # TRExplorer v2.1's GRCh38 tandem repeat catalog (broadinstitute/trexplorer-catalog) # GENCODE 50's basic annotation (the release scripts/lib/common.sh pins) @@ -35,8 +34,12 @@ # GenCC disease curation; name ALIAS|PSEUDOGENE|its HGNC id| # the HGNC id of the gene that holds ALIAS as its symbol # -# Needs curl, python3 and sha256sum or shasum. It reads about 1 GB of text, -# so run it on a box with room, not in a unit test. +# A gene's span is that of the first GENCODE record of its HGNC id, and the +# _PAR_Y copies of the pseudoautosomal genes are skipped, as +# build-callability.sh does. The rules are in og dev build sv-rule-regions +# (rust/crates/og-dev/src/build/sv_rules.rs). Needs curl, og +# (scripts/setup.sh --og) and sha256sum or shasum. It reads about 3 GB of +# text in a few seconds and keeps little of it. set -euo pipefail REPO=$(cd "$(dirname "$0")/../.." && pwd) @@ -66,161 +69,19 @@ if [ "${1:-}" = --list ]; then exit 0 fi -sha256() { - if command -v sha256sum >/dev/null 2>&1; then sha256sum "$1" | cut -d' ' -f1; else shasum -a 256 "$1" | cut -d' ' -f1; fi -} - +# shellcheck source=scripts/ci/lib/fetch-pinned.sh +. "${REPO}/scripts/ci/lib/fetch-pinned.sh" WORK=$(mktemp -d) trap 'rm -rf "$WORK"' EXIT +DIR=${SV_RULE_CACHE:-$WORK} [ -z "$CHECK" ] || OUT="${WORK}/rebuilt.tsv" -# get NAME URL SHA256: WORK/NAME from SV_RULE_CACHE when it matches, else -# downloaded; stop unless it matches. -get() { - local cached="${SV_RULE_CACHE:-}/$1" got - if [ -n "${SV_RULE_CACHE:-}" ] && [ -f "$cached" ] && [ "$(sha256 "$cached")" = "$3" ]; then - cp "$cached" "${WORK}/$1" - return - fi - echo "Downloading $2" >&2 - # Tries spread over about five minutes: EBI's FTP server refused every - # connection for minutes at a time (curl exit 7, 2026-11-06). - curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "${WORK}/$1" "$2" - got=$(sha256 "${WORK}/$1") - if [ "$got" != "$3" ]; then - echo "ERROR: $1 has sha256 ${got}, want $3" >&2 - exit 1 - fi - if [ -n "${SV_RULE_CACHE:-}" ]; then - mkdir -p "$SV_RULE_CACHE" && cp "${WORK}/$1" "$cached" - fi -} -get trx.bed.gz "$TRX_URL" "$TRX_SHA256" -get gencode.gtf.gz "$GENCODE_URL" "$GENCODE_SHA256" -get hgnc.txt "$HGNC_URL" "$HGNC_SHA256" -get regions.txt "$REGIONS_URL" "$REGIONS_SHA256" - -python3 - "$WORK" "${OUT}.part" <<'PY' -import bisect -import csv -import gzip -import re -import sys -from collections import defaultdict - -work, out = sys.argv[1:3] -CHROMS = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY"] -ORDER = {c: i for i, c in enumerate(CHROMS)} -MIN_TR = 50 -LOCUS_GAP = 1_000_000 -ATT = re.compile(r'(\S+) "([^"]*)"') -SRC_TR = "TRExplorer v2.1; GENCODE 50 MANE Select CDS" -SRC_IG = "GENCODE 50 IG and TR genes" -SRC_MHC = "GRCh38 assembly regions (GCA_000001405.15)" -SRC_ALIAS = "HGNC 2026-10-06; GENCODE 50" - -rows = [] -gene_span = {} # HGNC id -> (chrom, start0, end) -ig = defaultdict(list) # locus -> [(chrom, start0, end)] -cds = defaultdict(list) # chrom -> [(start0, end, gene)] -with gzip.open(f"{work}/gencode.gtf.gz", "rt") as f: - for line in f: - if line.startswith("#"): - continue - c = line.rstrip("\n").split("\t") - if c[0] not in ORDER: - continue - start, end = int(c[3]) - 1, int(c[4]) - if c[2] == "gene": - a = dict(ATT.findall(c[8])) - if a.get("hgnc_id"): - gene_span[a["hgnc_id"]] = (c[0], start, end) - name, gtype = a.get("gene_name", ""), a.get("gene_type", "") - if re.match(r"(IG|TR)_", gtype) and re.match(r"(IG[HKL]|TR[ABDG])", name): - locus = "TRA" if name.startswith("TRD") else name[:3] - ig[locus].append((c[0], start, end)) - elif c[2] == "CDS" and 'tag "MANE_Select"' in c[8]: - cds[c[0]].append((start, end, dict(ATT.findall(c[8]))["gene_name"])) - -# tandem_repeat: catalog repeats of MIN_TR bp or more on a MANE Select coding exon -index = {} -for chrom, iv in cds.items(): - iv.sort() - starts = [s for s, _, _ in iv] - maxend, m = [], 0 - for _, e, _ in iv: - m = max(m, e) - maxend.append(m) - index[chrom] = (iv, starts, maxend) -with gzip.open(f"{work}/trx.bed.gz", "rt") as f: - for line in f: - c = line.rstrip("\n").split("\t") - if c[0] not in index: - continue - s, e = int(c[1]), int(c[2]) - if e - s < MIN_TR: - continue - iv, starts, maxend = index[c[0]] - j = bisect.bisect_left(starts, e) - 1 - genes = set() - while j >= 0 and maxend[j] > s: - if iv[j][1] > s: - genes.add(iv[j][2]) - j -= 1 - if genes: - motif = c[3] if len(c[3]) <= 12 else f"{len(c[3])}bp motif" - rows.append((c[0], s, e, "tandem_repeat", "class_3", f"{'/'.join(sorted(genes))}:{motif}", SRC_TR)) - -# ig_tr_locus: per locus the largest group of genes on one chromosome -for locus, genes in ig.items(): - groups = [] - for chrom, s, e in sorted(genes): - g = groups[-1] if groups else None - if g and g[0] == chrom and s - g[2] < LOCUS_GAP: - g[2] = max(g[2], e) - g[3] += 1 - else: - groups.append([chrom, s, e, 1]) - chrom, s, e, _ = max(groups, key=lambda g: g[3]) - rows.append((chrom, s, e, "ig_tr_locus", "class_3", locus, SRC_IG)) - -# mhc: the assembly report's region (1-based, inclusive) -with open(f"{work}/regions.txt") as f: - mhc = sorted({tuple(l.split("\t")[:4]) for l in f if l.startswith("MHC\t")}) -if len(mhc) != 1: - sys.exit(f"ERROR: want one MHC region in the assembly report, got {mhc}") -_, chrom, s, e = mhc[0] -rows.append((f"chr{chrom}", int(s) - 1, int(e), "mhc", "class_3", "MHC", SRC_MHC)) - -# gene_alias: a pseudogene that carries a disease gene's symbol as alias or previous symbol -with open(f"{work}/hgnc.txt", newline="") as f: - hgnc = list(csv.DictReader(f, delimiter="\t")) -approved = {r["symbol"]: r for r in hgnc} -for r in hgnc: - if r["locus_group"] != "pseudogene" or r["hgnc_id"] not in gene_span: - continue - for sym in sorted({x.strip('"') for x in (r["alias_symbol"] + "|" + r["prev_symbol"]).split("|")} - {""}): - d = approved.get(sym) - if d and d["hgnc_id"] != r["hgnc_id"] and d["gencc"]: - chrom, s, e = gene_span[r["hgnc_id"]] - rows.append((chrom, s, e, "gene_alias", "remap_gene", f"{sym}|{r['symbol']}|{r['hgnc_id']}|{d['hgnc_id']}", - SRC_ALIAS)) +fetch_pinned "$DIR" trx.bed.gz "$TRX_URL" "$TRX_SHA256" || exit 1 +fetch_pinned "$DIR" gencode.gtf.gz "$GENCODE_URL" "$GENCODE_SHA256" || exit 1 +fetch_pinned "$DIR" hgnc.txt "$HGNC_URL" "$HGNC_SHA256" || exit 1 +fetch_pinned "$DIR" regions.txt "$REGIONS_URL" "$REGIONS_SHA256" || exit 1 -rows.sort(key=lambda x: (ORDER[x[0]], x[1], x[2], x[3], x[5])) -counts = defaultdict(int) -for x in rows: - counts[x[3]] += 1 -with open(out, "w") as f: - f.write("# assets/sv_rule_regions.tsv: the GRCh38 regions where og sv postfilter overrides AnnotSV's class 4\n" - "# or 5. Written by scripts/ci/build-sv-rule-regions.sh, which says how each rule's rows are chosen;\n" - "# rebuild it there, do not edit it by hand. start is 0-based and end 1-based, as in a BED file.\n" - "# Rows: " + ", ".join(f"{k} {counts[k]}" for k in sorted(counts)) + ".\n") - f.write("#chrom\tstart\tend\trule\taction\tname\tsource\n") - for x in rows: - f.write("\t".join(map(str, x)) + "\n") -print("Rows: " + ", ".join(f"{k} {counts[k]}" for k in sorted(counts)), file=sys.stderr) -PY -mv "${OUT}.part" "$OUT" +"${REPO}/bin/og" dev build sv-rule-regions --work "$DIR" --out "$OUT" if [ -n "$CHECK" ]; then if cmp -s "$OUT" "$CHECK"; then echo "${CHECK} is what the pinned downloads give, byte for byte" >&2 diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index f76068dcf..fe64ffe70 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -62,7 +62,7 @@ bash-step-20-numt assets/callability_genes.tsv bin/collect_summary.py bin/traits bash-step-bash42 assets/construct_variants.tsv bin/clip_check.py bin/cpsr_sample_id bin/og scripts/03d-octopus.sh scripts/04a-tiddit.sh scripts/10-telomere-hunter.sh scripts/17-cpsr.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh bloodtype-1 assets/bloodtype_markers.tsv bin/cn_loci.tsv bin/og modules/local/bloodtype/* modules/local/vcf_precheck/* scripts/41-bloodtype.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh bumps-drop-in-baseline scripts/ci/fake-docker/lib.sh -callability-1-catalog assets/callability_genes.tsv assets/paralog_loci.tsv assets/pgx_genes.tsv bin/cn_loci.tsv bin/og scripts/ci/build-callability.sh scripts/ci/fake-docker/lib.sh +callability-1-catalog assets/callability_genes.tsv assets/paralog_loci.tsv assets/pgx_genes.tsv bin/cn_loci.tsv bin/og scripts/ci/build-callability.sh scripts/ci/fake-docker/lib.sh scripts/ci/lib/fetch-pinned.sh callability-2-genes assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/stub/sample.vcf.gz bin/collect_summary.py bin/render_report.py bin/results_index.py bin/traits_report.py modules/local/html_report/* modules/local/mosdepth/* modules/local/vcf_precheck/* scripts/16b-mosdepth.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh chip-clonal-1 assets/chip_genes.tsv bin/chip_report.py bin/og bin/render_report.py modules/local/chip/* modules/local/vcf_precheck/* scripts/49-chip.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh chip-mt-1-roundtrip assets/par_grch38.bed bin/og scripts/40-export-consumer.sh scripts/chip-to-vcf.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz @@ -75,7 +75,7 @@ clinvar-combined-allele-1 assets/construct_variants.tsv assets/deep_intronic_sit clinvar-cpsr-quality-1 assets/construct_variants.tsv assets/deep_intronic_sites.tsv bin/clip_check.py bin/og bin/render_report.py scripts/06-clinvar-screen.sh scripts/ci/fake-docker/lib.sh cn-engine-1-parity assets/bloodtype_markers.tsv bin/cn_loci.tsv bin/copy_number.py bin/og bin/traits_report.py scripts/41-bloodtype.sh scripts/46-lpa.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh cn-guard-1-mapq0 bin/cn_loci.tsv bin/copy_number.py bin/og bin/traits_report.py modules/local/parascopy/* scripts/35-paralogs.sh -cn-intervals-1-catalog bin/cn_loci.tsv bin/og scripts/ci/build-cn-intervals.py scripts/ci/fake-docker/lib.sh +cn-intervals-1-catalog assets/parascopy_models_v1.7.txt bin/cn_loci.tsv bin/og scripts/ci/build-cn-intervals.sh scripts/ci/fake-docker/lib.sh scripts/ci/lib/fetch-pinned.sh cn-intervals-2-hg002 assets/par_grch38.bed assets/paralog_loci.tsv bin/cn_loci.tsv bin/copy_number.py bin/og bin/traits_report.py modules/local/deepvariant/* modules/local/parascopy/* modules/local/vcf_precheck/* scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh cn-report-1-tables assets/carrier_genes.tsv assets/par_grch38.bed bin/collect_summary.py bin/couple_report.py bin/family_report.py bin/og bin/paraphase_report.py bin/render_report.py bin/results_index.py bin/traits_report.py scripts/24-html-report.sh scripts/39-couple.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh comphet-phasing-1-reads assets/construct_variants.tsv assets/stub/sample.bam assets/stub/sample.bam.bai bin/clip_check.py bin/og modules/local/slivar/* modules/local/vcf_precheck/* modules/local/vep/* scripts/13-vep-annotation.sh scripts/31-slivar.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/comphet/make_fixture.py @@ -91,7 +91,7 @@ expansions-and-rfc1-1-motif assets/str_extra_loci.json assets/stub/db/catalog.js families-file-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/families.example.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/pedigree.example.tsv assets/pgx_genes.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat assets/stub/pedigree.tsv assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz bin/clip_check.py bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/couple_report.py bin/family_report.py bin/og bin/pgx_outside_calls.py bin/render_report.py bin/results_index.py bin/traits_report.py modules/local/align_minimap2/* modules/local/clinvar_screen/* modules/local/couple/* modules/local/deepvariant/* modules/local/fastp/* modules/local/hla_typing/* modules/local/html_report/* modules/local/t1k_build/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/06-clinvar-screen.sh scripts/08-hla-typing.sh scripts/24-html-report.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/38-family.sh scripts/39-couple.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/families/stub_families.tsv family-joint-1 assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/families.example.tsv assets/par_grch38.bed assets/pedigree.example.tsv assets/stub/pedigree.tsv assets/stub/sample.bam assets/stub/sample.bam.bai bin/clip_check.py bin/collect_summary.py bin/family_report.py bin/og bin/traits_report.py modules/local/align_minimap2/* modules/local/clinvar_screen/* modules/local/deepvariant/* modules/local/fastp/* modules/local/roh/* modules/local/somalier/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/06-clinvar-screen.sh scripts/11-roh-analysis.sh scripts/33-sample-qc.sh scripts/38-family.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh fixitfelix-chr21-1 assets/callability_genes.tsv scripts/16b-mosdepth.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh -followups-ci-3-sv-rule-regions assets/sv_rule_regions.tsv scripts/ci/build-sv-rule-regions.sh scripts/ci/fake-docker/lib.sh +followups-ci-3-sv-rule-regions assets/sv_rule_regions.tsv bin/og scripts/ci/build-sv-rule-regions.sh scripts/ci/fake-docker/lib.sh scripts/ci/lib/fetch-pinned.sh followups-clinvar-traits-1-rules-html assets/genotype_rules.tsv assets/traits.tsv bin/render_report.py bin/traits_report.py scripts/42-traits.sh scripts/ci/check-rule-alleles.sh followups-cn-1-couple-report assets/carrier_genes.tsv assets/par_grch38.bed assets/stub/pedigree.tsv bin/collect_summary.py bin/couple_report.py bin/family_report.py bin/og bin/render_report.py bin/results_index.py bin/traits_report.py modules/local/couple/* modules/local/html_report/* scripts/39-couple.sh scripts/ci/fake-docker/lib.sh followups-nf-1-nf-wiring assets/stub/* assets/test_samplesheet.csv assets/steps.tsv bin/collect_summary.py bin/render_report.py bin/results_index.py bin/run_state.py bin/traits_report.py conf/test.config conf/test_all.config conf/containers.config modules/local/html_report/* modules/local/couple/* modules/local/parascopy/* modules/local/deepvariant/* modules/local/cpsr/* modules/local/clinvar_screen/* modules/local/locityper/* workflows/* tests/fixtures/consumer_export/sites.tsv.gz bin/og diff --git a/scripts/ci/e2e-run.sh b/scripts/ci/e2e-run.sh index 8d4bfe1c4..4a5a45800 100755 --- a/scripts/ci/e2e-run.sh +++ b/scripts/ci/e2e-run.sh @@ -196,7 +196,7 @@ PGS000066 b34f31b8b3a7dbbbb58882a70e14649d"} # EBI's FTP server, which refused every connection on 2026-11-06 (run # 42259432681). CATALOGS="callability bash ${REPO}/scripts/ci/build-callability.sh -cn-intervals python3 ${REPO}/scripts/ci/build-cn-intervals.py +cn-intervals bash ${REPO}/scripts/ci/build-cn-intervals.sh sv-rule-regions bash ${REPO}/scripts/ci/build-sv-rule-regions.sh" # catalog_files: "FOLDER NAME URL SHA256" of every download of CATALOGS. catalog_files() { diff --git a/scripts/ci/lib/fetch-pinned.sh b/scripts/ci/lib/fetch-pinned.sh new file mode 100644 index 000000000..ccdb82f01 --- /dev/null +++ b/scripts/ci/lib/fetch-pinned.sh @@ -0,0 +1,61 @@ +# shellcheck shell=bash +# fetch-pinned.sh: the one download rule of the asset builders in +# scripts/ci/ (build-callability.sh, build-cn-intervals.sh, +# build-sv-rule-regions.sh). Source it; it defines two functions. +# +# fetch_pinned DIR NAME URL SHA256 +# Makes DIR/NAME the file at URL with that sha256. A file already there +# with the right sha256 is kept, so DIR can be a cache that outlives the +# run. Otherwise it prints "Downloading URL" on stderr (the e2e cases +# count these lines), downloads to DIR/NAME.part and renames it when the +# sha256 matches. A 403 or 404 fails at once: a missing or forbidden file +# stays so. Any other failure is tried again 6 times, 30 s apart, for +# about five minutes: EBI's FTP server refused every connection for +# minutes at a time on 2026-11-06 (curl exit 7). Returns 1 with an ERROR +# line when the file cannot be had. +# +# pinned_sha256 FILE +# The sha256 of FILE, from sha256sum or shasum. +# +# og never downloads: the builders fetch here, then hand the folder to +# og dev build. + +pinned_sha256() { + if command -v sha256sum > /dev/null 2>&1; then + sha256sum "$1" | cut -d' ' -f1 + else + shasum -a 256 "$1" | cut -d' ' -f1 + fi +} + +fetch_pinned() { + local dir=$1 name=$2 url=$3 want=$4 file code got + file="${dir}/${name}" + if [ -f "$file" ] && [ "$(pinned_sha256 "$file")" = "$want" ]; then + return 0 + fi + mkdir -p "$dir" || return 1 + echo "Downloading ${url}" >&2 + if ! code=$(curl -fsSL -w '%{http_code}' -o "${file}.part" "$url"); then + case $code in + 403 | 404) + rm -f "${file}.part" + echo "ERROR: cannot download ${url}: HTTP ${code}" >&2 + return 1 + ;; + esac + if ! curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 \ + -o "${file}.part" "$url"; then + rm -f "${file}.part" + echo "ERROR: cannot download ${url}" >&2 + return 1 + fi + fi + got=$(pinned_sha256 "${file}.part") + if [ "$got" != "$want" ]; then + rm -f "${file}.part" + echo "ERROR: ${name} has sha256 ${got}, want ${want}" >&2 + return 1 + fi + mv "${file}.part" "$file" +} diff --git a/tests/e2e/cn-intervals-1-catalog.sh b/tests/e2e/cn-intervals-1-catalog.sh index 5b0db18cb..dd2a6262f 100755 --- a/tests/e2e/cn-intervals-1-catalog.sh +++ b/tests/e2e/cn-intervals-1-catalog.sh @@ -1,5 +1,5 @@ #!/usr/bin/env bash -# The interval rows of bin/cn_loci.tsv are what scripts/ci/build-cn-intervals.py +# The interval rows of bin/cn_loci.tsv are what scripts/ci/build-cn-intervals.sh # writes from its pinned downloads (ClinGen's recurrent CNV BED V2.1 and three # GIAB v3.6 stratifications, each checked against its sha256): a hand edit of # an interval row fails here. Control: a copy of the catalog with one target @@ -11,9 +11,9 @@ CACHE="${CASE_TMP}/downloads" [ -z "${E2E_DATA_CACHE:-}" ] || CACHE="${E2E_DATA_CACHE}/catalogs/cn-intervals" mkdir -p "$CACHE" -echo "+ scripts/ci/build-cn-intervals.py --check" -python3 "${REPO}/scripts/ci/build-cn-intervals.py" --cache "$CACHE" --check 2>&1 | tee "$STEP_LOG" -check_eq "build-cn-intervals.py --check exits 0 on the committed catalog" "${PIPESTATUS[0]}" 0 +echo "+ scripts/ci/build-cn-intervals.sh --check" +bash "${REPO}/scripts/ci/build-cn-intervals.sh" --cache "$CACHE" --check 2>&1 | tee "$STEP_LOG" +check_eq "build-cn-intervals.sh --check exits 0 on the committed catalog" "${PIPESTATUS[0]}" 0 check "the log says all 90 interval rows were checked" grep -q "the 90 interval rows are what the pinned downloads give" "$STEP_LOG" # Control: one interval target moved by 1 kb @@ -25,7 +25,7 @@ awk -F'\t' -v OFS='\t' '$1 == "22q11.2_DGS_proximal_A-B" && !done { "${REPO}/bin/cn_loci.tsv" > "$EDITED" # shellcheck disable=SC2016 # the positional arguments of bash -c check "control: the edited copy differs from the catalog" bash -c '! cmp -s "$1" "$2"' _ "$EDITED" "${REPO}/bin/cn_loci.tsv" -python3 "${REPO}/scripts/ci/build-cn-intervals.py" --cache "$CACHE" --catalog "$EDITED" --check > "${CASE_TMP}/edited.log" 2>&1 +bash "${REPO}/scripts/ci/build-cn-intervals.sh" --cache "$CACHE" --catalog "$EDITED" --check > "${CASE_TMP}/edited.log" 2>&1 check_eq "control: --check fails (exit 1) on a catalog with one interval target moved" "$?" 1 cat "${CASE_TMP}/edited.log" diff --git a/tests/test_copy_number.py b/tests/test_copy_number.py index b00a00302..337dff70a 100644 --- a/tests/test_copy_number.py +++ b/tests/test_copy_number.py @@ -36,7 +36,7 @@ each, and a row shifted by 1 kb no longer matches its fixture; the PARASCOPY module runs the same list; the CCL3L1 reading names no disease (a planted 'HIV' is caught); - 9. the interval rows (scripts/ci/build-cn-intervals.py): a small ClinGen + 9. the interval rows (og dev build cn-intervals, with og): a small ClinGen BED with a 'chrx' row converts to rows on chrX that load, the same rows left as 'chrx' are refused, and an unknown contig stops the build; the committed rows are the ClinGen regions plus the named hand rows, with @@ -684,7 +684,6 @@ def dosage(cat, read): # Section 9: the interval rows. -BUILD = os.path.join(REPO, "scripts", "ci", "build-cn-intervals.py") SMALL_BED = """track name="test" #chrom\tstart\tend\tname\tscore\tstrand\tthickStart\tthickEnd\titemRgb chr22\t18159724\t18924717\tBP-A\t0\t.\t18159724\t18924717\t0,0,0 @@ -724,39 +723,49 @@ def depths(r, target=30.0, control=30.0): return {n: target if n in tg else control for *_, n in cn.regions(r)} -def intervals(cat, read): - import importlib.util - spec = importlib.util.spec_from_file_location("build_cn_intervals", BUILD) - bld = importlib.util.module_from_spec(spec) - spec.loader.exec_module(bld) - rows = bld.build(SMALL_BED.splitlines(True), [], [], [], []) - check("converter: the two region rows of the small BED, its chrx row on chrX", - [(r["locus"], r["chrom"]) for r in rows[:2]] == [("22q11.2_DGS_proximal_A-D", "chr22"), ("Xp22.31", "chrX")], - [(r["locus"], r["chrom"]) for r in rows]) +def build_rows(bed): + """og dev build cn-intervals --stdout on a ClinGen BED and no masks: (exit code, rows as dicts, stderr).""" with tempfile.TemporaryDirectory() as t: - p = os.path.join(t, "c.tsv") - write_catalog(p, rows) - try: - got = len(cn.load_catalog(p)) - except cn.CatalogError as e: - got = str(e) - check("converter: its rows load", got == len(rows), got) - write_catalog(p, [dict(r, chrom="chrx") if r["chrom"] == "chrX" else r for r in rows]) - try: - cn.load_catalog(p) - got = "loaded" - except cn.CatalogError as e: - got = str(e) - check("converter: the same rows left as 'chrx' are refused", "chrom 'chrx'" in got, got) - try: - bld.build(SMALL_BED.replace("chrx", "chrQ").splitlines(True), [], [], [], []) - got = "built" - except bld.BuildError as e: - got = str(e) - check("converter: an unknown contig stops the build", "'chrQ' is not a GRCh38 chromosome" in got, got) - check("converter: no target piece overlaps a breakpoint row", - all(not (s < 18924717 and 18159724 < e) and not (s < 21562827 and 21111384 < e) for r in rows[:1] - for _, s, e, _ in cn.parse_row(r)["_targets"])) + p = os.path.join(t, "clingen.bed") + with open(p, "w") as f: + f.write(bed) + out = subprocess.run([OG, "dev", "build", "cn-intervals", "--clingen", p, "--segdups", os.devnull, + "--gaps", os.devnull, "--satellites", os.devnull, "--parascopy", os.devnull, + "--stdout"], capture_output=True, text=True) + rows = [dict(zip(cn.COLUMNS, line.split("\t"))) for line in out.stdout.splitlines() if not line.startswith("#")] + return out.returncode, rows, out.stderr + + +def intervals(cat, read): + if not OG: + print("skip: OG is not set: the converter checks of og dev build cn-intervals did not run") + else: + code, rows, err = build_rows(SMALL_BED) + check("converter: the two region rows of the small BED, its chrx row on chrX", + code == 0 and [(r["locus"], r["chrom"]) for r in rows[:2]] + == [("22q11.2_DGS_proximal_A-D", "chr22"), ("Xp22.31", "chrX")], + (code, err, [(r["locus"], r["chrom"]) for r in rows])) + with tempfile.TemporaryDirectory() as t: + p = os.path.join(t, "c.tsv") + write_catalog(p, rows) + try: + got = len(cn.load_catalog(p)) + except cn.CatalogError as e: + got = str(e) + check("converter: its rows load", got == len(rows), got) + write_catalog(p, [dict(r, chrom="chrx") if r["chrom"] == "chrX" else r for r in rows]) + try: + cn.load_catalog(p) + got = "loaded" + except cn.CatalogError as e: + got = str(e) + check("converter: the same rows left as 'chrx' are refused", "chrom 'chrx'" in got, got) + code, _, err = build_rows(SMALL_BED.replace("chrx", "chrQ")) + check("converter: an unknown contig stops the build", + code == 2 and "'chrQ' is not a GRCh38 chromosome" in err, (code, err)) + check("converter: no target piece overlaps a breakpoint row", + all(not (s < 18924717 and 18159724 < e) and not (s < 21562827 and 21111384 < e) for r in rows[:1] + for _, s, e, _ in cn.parse_row(r)["_targets"])) iv = [r for r in cat.values() if r["kind"] == "interval" and r["_passes"] == ["."]] names = [r["locus"] for r in iv] diff --git a/tests/test_e2e_data_cache.sh b/tests/test_e2e_data_cache.sh index bfaf97a2d..e403f51ca 100755 --- a/tests/test_e2e_data_cache.sh +++ b/tests/test_e2e_data_cache.sh @@ -74,7 +74,7 @@ else echo "FAIL: the pgs key"; failed=1; fi if "$RUN" --data-key catalogs | grep -qE '^e2e-data-catalogs-[0-9a-f]{16}$'; then echo "PASS: the catalogs key is a hash" else echo "FAIL: the catalogs key"; failed=1; fi # Every row must be "NAME URL SHA256": a malformed row would reach the fill. -list=$( { bash "${REPO}/scripts/ci/build-callability.sh" --list; python3 "${REPO}/scripts/ci/build-cn-intervals.py" --list; bash "${REPO}/scripts/ci/build-sv-rule-regions.sh" --list; } ) +list=$( { bash "${REPO}/scripts/ci/build-callability.sh" --list; bash "${REPO}/scripts/ci/build-cn-intervals.sh" --list; bash "${REPO}/scripts/ci/build-sv-rule-regions.sh" --list; } ) n=$(awk 'NF == 3 && $2 ~ /^https:/ && $3 ~ /^[0-9a-f]{64}$/' <<< "$list" | wc -l | tr -d ' ') all=$(wc -l <<< "$list" | tr -d ' ') if [ "$n" = 15 ] && [ "$all" = 15 ]; then echo "PASS: the three builders list 15 pinned downloads (7, 4 and 4), every row well formed" diff --git a/versions.env b/versions.env index b67c1e177..304024130 100644 --- a/versions.env +++ b/versions.env @@ -407,7 +407,7 @@ EXOMISER_PHENOTYPE_SHA256="aa923530ca43cf5051ed2ba7b221bb73deb2b0dd66695b2ba6046 # The recurrent CNV regions of the interval rows of bin/cn_loci.tsv (step 35): # ClinGen Dosage Sensitivity's recurrent CNV BED, version 2.1, GRCh38 (CC0), # checked against the sha256 of the file as downloaded on 2026-11-04. -# scripts/ci/build-cn-intervals.py turns it into the rows; the e2e case +# scripts/ci/build-cn-intervals.sh turns it into the rows; the e2e case # tests/e2e/cn-intervals-1-catalog.sh rebuilds them and fails on any # difference. Change the two together, then rerun the script. CLINGEN_RECURRENT_CNV_URL="https://ftp.clinicalgenome.org/ClinGen_recurrent_CNV_V2.1-hg38.bed" From 1f2b7fcadaf904c7d42df9284c78f1e9b679dfe2 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 12:29:25 +0100 Subject: [PATCH 3/7] test(dev): plant three wrong builder rules to show the tests go red Red-first control, reverted in the next commit: callability keeps the _PAR_Y copies, the SV rule builder keeps the last record per HGNC id, and each cn-intervals piece ends one base short. Bead: pgp-bs9u.35 --- rust/crates/og-dev/src/build/callability.rs | 3 ++- rust/crates/og-dev/src/build/cn_intervals.rs | 3 ++- rust/crates/og-dev/src/build/sv_rules.rs | 5 +++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/rust/crates/og-dev/src/build/callability.rs b/rust/crates/og-dev/src/build/callability.rs index d7ce4dece..556c65865 100644 --- a/rust/crates/og-dev/src/build/callability.rs +++ b/rust/crates/og-dev/src/build/callability.rs @@ -111,7 +111,8 @@ fn read_gtf(path: &Path, alt: &HashMap) -> Result Vec { // Ceiling divisions of positive numbers. let n = (len + PIECE - 1) / PIECE; let step = (len + n - 1) / n; - out.extend((0..n).map(|i| (s + i * step, e.min(s + (i + 1) * step)))); + // PLANTED, reverted in the next commit: each piece ends one base short. + out.extend((0..n).map(|i| (s + i * step, e.min(s + (i + 1) * step) - 1))); } out } diff --git a/rust/crates/og-dev/src/build/sv_rules.rs b/rust/crates/og-dev/src/build/sv_rules.rs index 7004f02b7..1007fc81b 100644 --- a/rust/crates/og-dev/src/build/sv_rules.rs +++ b/rust/crates/og-dev/src/build/sv_rules.rs @@ -99,9 +99,10 @@ fn read_gtf(path: &Path) -> Result { b"gene" => { let (start, end) = r.span()?; if let Some(h) = r.attr("hgnc_id")?.filter(|h| !h.is_empty()) - && !r.attr("gene_id")?.is_some_and(|id| id.ends_with("_PAR_Y")) + && !r.attr("gene_id")?.is_some_and(|id| id.ends_with("_PAR_X")) { - g.span.entry(h.to_owned()).or_insert((chrom, start, end)); + // PLANTED, reverted in the next commit: the last record wins. + g.span.insert(h.to_owned(), (chrom, start, end)); } let (name, ty) = ( r.attr("gene_name")?.unwrap_or(""), From 1b8998fff3c1710a5194bf6b8bbee9e58fc0113b Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 12:59:37 +0100 Subject: [PATCH 4/7] Revert "test(dev): plant three wrong builder rules to show the tests go red" The plant turned five builder tests red in Guard's rust job (run 43417543846); back to the right rules. Bead: pgp-bs9u.35 --- rust/crates/og-dev/src/build/callability.rs | 3 +-- rust/crates/og-dev/src/build/cn_intervals.rs | 3 +-- rust/crates/og-dev/src/build/sv_rules.rs | 5 ++--- 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/rust/crates/og-dev/src/build/callability.rs b/rust/crates/og-dev/src/build/callability.rs index 556c65865..d7ce4dece 100644 --- a/rust/crates/og-dev/src/build/callability.rs +++ b/rust/crates/og-dev/src/build/callability.rs @@ -111,8 +111,7 @@ fn read_gtf(path: &Path, alt: &HashMap) -> Result Vec { // Ceiling divisions of positive numbers. let n = (len + PIECE - 1) / PIECE; let step = (len + n - 1) / n; - // PLANTED, reverted in the next commit: each piece ends one base short. - out.extend((0..n).map(|i| (s + i * step, e.min(s + (i + 1) * step) - 1))); + out.extend((0..n).map(|i| (s + i * step, e.min(s + (i + 1) * step)))); } out } diff --git a/rust/crates/og-dev/src/build/sv_rules.rs b/rust/crates/og-dev/src/build/sv_rules.rs index 1007fc81b..7004f02b7 100644 --- a/rust/crates/og-dev/src/build/sv_rules.rs +++ b/rust/crates/og-dev/src/build/sv_rules.rs @@ -99,10 +99,9 @@ fn read_gtf(path: &Path) -> Result { b"gene" => { let (start, end) = r.span()?; if let Some(h) = r.attr("hgnc_id")?.filter(|h| !h.is_empty()) - && !r.attr("gene_id")?.is_some_and(|id| id.ends_with("_PAR_X")) + && !r.attr("gene_id")?.is_some_and(|id| id.ends_with("_PAR_Y")) { - // PLANTED, reverted in the next commit: the last record wins. - g.span.insert(h.to_owned(), (chrom, start, end)); + g.span.entry(h.to_owned()).or_insert((chrom, start, end)); } let (name, ty) = ( r.attr("gene_name")?.unwrap_or(""), From cebf6b944a8c14b5d8deac6cfd689a63b24b212e Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 13:02:06 +0100 Subject: [PATCH 5/7] fix(ci): a fetch of its own per run, and a stalled transfer counts as failed fetch_pinned downloads to a mktemp file beside the target, so two runs on one cache never write the same partial file, and gives curl a 30 s connect timeout and a 1 KB/s for 60 s floor, so a stalled server reaches the retry or the error instead of holding the build. Bead: pgp-bs9u.35 --- CHANGELOG.md | 2 +- scripts/ci/lib/fetch-pinned.sh | 36 ++++++++++++++++++++-------------- 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7ce7f1507..20ebc7ff0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,7 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n ## Unreleased - 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; any other failure is tried for about five minutes, 30 s apart. `build-cn-intervals.py` tried three times, 30 s apart. + - 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; any other failure, a stalled transfer included, is tried for about five minutes, 30 s apart. `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. - `build-cn-intervals.sh --check` counts the lines only in the file and the lines only in the rebuild. The Python counted only the lines the rebuild lacks, so a catalog with an added line said "0 line(s) differ". - The SV rule builder takes each gene's span from the first GENCODE record of its HGNC id and skips the `_PAR_Y` copies, as the callability builder does. The Python kept the last record, so a pseudogene in a pseudoautosomal region would have taken its chrY span. No gene of the pinned files is one, so the table does not change. diff --git a/scripts/ci/lib/fetch-pinned.sh b/scripts/ci/lib/fetch-pinned.sh index ccdb82f01..f05b1c724 100644 --- a/scripts/ci/lib/fetch-pinned.sh +++ b/scripts/ci/lib/fetch-pinned.sh @@ -7,12 +7,15 @@ # Makes DIR/NAME the file at URL with that sha256. A file already there # with the right sha256 is kept, so DIR can be a cache that outlives the # run. Otherwise it prints "Downloading URL" on stderr (the e2e cases -# count these lines), downloads to DIR/NAME.part and renames it when the -# sha256 matches. A 403 or 404 fails at once: a missing or forbidden file -# stays so. Any other failure is tried again 6 times, 30 s apart, for -# about five minutes: EBI's FTP server refused every connection for -# minutes at a time on 2026-11-06 (curl exit 7). Returns 1 with an ERROR -# line when the file cannot be had. +# count these lines), downloads to a temporary file of its own beside +# DIR/NAME (so two runs on one cache never share one) and renames it +# when the sha256 matches. A 403 or 404 fails at once: a missing or +# forbidden file stays so. Any other failure is tried again 6 times, +# 30 s apart, for about five minutes: EBI's FTP server refused every +# connection for minutes at a time on 2026-11-06 (curl exit 7). A +# connection that does not open in 30 s, or a transfer under 1 KB/s for +# 60 s, counts as a failure, so a stalled server cannot hold the build. +# Returns 1 with an ERROR line when the file cannot be had. # # pinned_sha256 FILE # The sha256 of FILE, from sha256sum or shasum. @@ -29,33 +32,36 @@ pinned_sha256() { } fetch_pinned() { - local dir=$1 name=$2 url=$3 want=$4 file code got + local dir=$1 name=$2 url=$3 want=$4 file part code got file="${dir}/${name}" if [ -f "$file" ] && [ "$(pinned_sha256 "$file")" = "$want" ]; then return 0 fi mkdir -p "$dir" || return 1 + part=$(mktemp "${file}.part.XXXXXX") || return 1 echo "Downloading ${url}" >&2 - if ! code=$(curl -fsSL -w '%{http_code}' -o "${file}.part" "$url"); then + if ! code=$(curl -fsSL --connect-timeout 30 --speed-limit 1024 --speed-time 60 \ + -w '%{http_code}' -o "$part" "$url"); then case $code in 403 | 404) - rm -f "${file}.part" + rm -f "$part" echo "ERROR: cannot download ${url}: HTTP ${code}" >&2 return 1 ;; esac - if ! curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 \ - -o "${file}.part" "$url"; then - rm -f "${file}.part" + if ! curl -fsSL --connect-timeout 30 --speed-limit 1024 --speed-time 60 \ + --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "$part" "$url"; then + rm -f "$part" echo "ERROR: cannot download ${url}" >&2 return 1 fi fi - got=$(pinned_sha256 "${file}.part") + got=$(pinned_sha256 "$part") if [ "$got" != "$want" ]; then - rm -f "${file}.part" + rm -f "$part" echo "ERROR: ${name} has sha256 ${got}, want ${want}" >&2 return 1 fi - mv "${file}.part" "$file" + chmod 644 "$part" + mv "$part" "$file" } From 0948aae258d8b6c383ccb2a49133d8e1bc119480 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 13:34:54 +0100 Subject: [PATCH 6/7] fix(ci): stop pinned downloads on every forbidden or missing response --- CHANGELOG.md | 2 +- rust/crates/og-dev/src/build/mod.rs | 130 ++++++++++++++++++++++++++++ scripts/ci/lib/fetch-pinned.sh | 14 +-- 3 files changed, 140 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 20ebc7ff0..6eee801af 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,7 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n ## Unreleased - 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; any other failure, a stalled transfer included, is tried for about five minutes, 30 s apart. `build-cn-intervals.py` tried three times, 30 s apart. + - 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. - `build-cn-intervals.sh --check` counts the lines only in the file and the lines only in the rebuild. The Python counted only the lines the rebuild lacks, so a catalog with an added line said "0 line(s) differ". - The SV rule builder takes each gene's span from the first GENCODE record of its HGNC id and skips the `_PAR_Y` copies, as the callability builder does. The Python kept the last record, so a pseudogene in a pseudoautosomal region would have taken its chrY span. No gene of the pinned files is one, so the table does not change. diff --git a/rust/crates/og-dev/src/build/mod.rs b/rust/crates/og-dev/src/build/mod.rs index b2098e3dc..d8e9c5931 100644 --- a/rust/crates/og-dev/src/build/mod.rs +++ b/rust/crates/og-dev/src/build/mod.rs @@ -296,3 +296,133 @@ mod tests { assert!(int("1.5", "x").is_err()); } } + +#[cfg(all(test, unix))] +mod fetch_tests { + use std::io::{BufRead, BufReader, Write}; + use std::net::TcpListener; + use std::process::Command; + use std::sync::atomic::{AtomicBool, AtomicU64, Ordering}; + use std::time::Duration; + + static NEXT: AtomicU64 = AtomicU64::new(0); + const BODY: &str = "hello\n"; + const SHA: &str = "5891b5b522d5df086d0ff0b110fbd9d21bb4fc7163af34d08286a2e846f6be03"; + + // Real HTTP responses exercise curl's own retries too. Shorten only the + // test delays; the production script still supplies its stall limits. + fn fetch(statuses: &[u16], sha: &str, cached: bool, succeeds: bool, expected: &[u16]) { + let dir = std::env::temp_dir().join(format!( + "og-fetch-{}-{}", + std::process::id(), + NEXT.fetch_add(1, Ordering::Relaxed) + )); + std::fs::create_dir_all(&dir).unwrap(); + if cached { + std::fs::write(dir.join("asset"), BODY).unwrap(); + } + let server = TcpListener::bind("127.0.0.1:0").unwrap(); + server.set_nonblocking(true).unwrap(); + let url = format!("http://{}/asset", server.local_addr().unwrap()); + let stop = AtomicBool::new(false); + let (output, seen) = std::thread::scope(|scope| { + let worker = scope.spawn(|| { + let mut seen = Vec::new(); + while !stop.load(Ordering::Relaxed) { + match server.accept() { + Ok((mut stream, _)) => { + stream.set_read_timeout(Some(Duration::from_secs(5))).unwrap(); + let mut reader = BufReader::new(&stream); + let mut line = String::new(); + loop { + line.clear(); + if reader.read_line(&mut line).unwrap() == 0 || line == "\r\n" { + break; + } + } + let code = statuses.get(seen.len()).copied().unwrap_or(200); + seen.push(code); + write!(stream, + "HTTP/1.1 {code} Test\r\nContent-Length: {}\r\nConnection: close\r\n\r\n{BODY}", + BODY.len()).unwrap(); + } + Err(e) if e.kind() == std::io::ErrorKind::WouldBlock => { + std::thread::sleep(Duration::from_millis(5)); + } + Err(e) => panic!("HTTP test server: {e}"), + } + } + seen + }); + let output = Command::new("bash") + .args([ + "-c", + r#" + source "$1" + sleep() { :; } + curl() { command curl "$@" --noproxy '*' --retry-delay 1 --max-time 5; } + fetch_pinned "$2" asset "$3" "$4" + "#, + "_", + ]) + .arg(concat!( + env!("CARGO_MANIFEST_DIR"), + "/../../../scripts/ci/lib/fetch-pinned.sh" + )) + .arg(&dir) + .arg(url) + .arg(sha) + .output() + .unwrap(); + stop.store(true, Ordering::Relaxed); + (output, worker.join().unwrap()) + }); + let file = std::fs::read_to_string(dir.join("asset")).ok(); + let files: Vec<_> = std::fs::read_dir(&dir) + .unwrap() + .map(|e| e.unwrap().file_name()) + .collect(); + std::fs::remove_dir_all(&dir).unwrap(); + assert_eq!(seen, expected, "HTTP sequence {statuses:?}"); + assert_eq!( + output.status.success(), + succeeds, + "{}", + String::from_utf8_lossy(&output.stderr) + ); + assert_eq!(file.as_deref(), succeeds.then_some(BODY)); + assert_eq!( + files.len(), + usize::from(succeeds), + "temporary download leaked: {files:?}" + ); + if let Some(code @ (403 | 404)) = expected.last() { + assert!(String::from_utf8_lossy(&output.stderr).contains(&format!("HTTP {code}"))); + } + } + + #[test] + fn http_403_stops_after_any_attempt() { + fetch(&[403], SHA, false, false, &[403]); + fetch(&[503, 403], SHA, false, false, &[503, 403]); + } + + #[test] + fn http_404_stops_after_any_attempt() { + fetch(&[404], SHA, false, false, &[404]); + fetch(&[503, 404], SHA, false, false, &[503, 404]); + } + + #[test] + fn other_failures_retry_but_stop_at_the_bound() { + fetch(&[503, 503, 200], SHA, false, true, &[503, 503, 200]); + fetch(&[503; 7], SHA, false, false, &[503; 7]); + } + + #[test] + fn downloads_are_pinned_and_cached() { + fetch(&[200], SHA, false, true, &[200]); + fetch(&[200], "wrong-sha256", false, false, &[200]); + fetch(&[], SHA, true, true, &[]); + } +} diff --git a/scripts/ci/lib/fetch-pinned.sh b/scripts/ci/lib/fetch-pinned.sh index f05b1c724..53a1fdc9b 100644 --- a/scripts/ci/lib/fetch-pinned.sh +++ b/scripts/ci/lib/fetch-pinned.sh @@ -33,6 +33,7 @@ pinned_sha256() { fetch_pinned() { local dir=$1 name=$2 url=$3 want=$4 file part code got + local attempt=0 started=$SECONDS file="${dir}/${name}" if [ -f "$file" ] && [ "$(pinned_sha256 "$file")" = "$want" ]; then return 0 @@ -40,8 +41,8 @@ fetch_pinned() { mkdir -p "$dir" || return 1 part=$(mktemp "${file}.part.XXXXXX") || return 1 echo "Downloading ${url}" >&2 - if ! code=$(curl -fsSL --connect-timeout 30 --speed-limit 1024 --speed-time 60 \ - -w '%{http_code}' -o "$part" "$url"); then + while ! code=$(curl -fsSL --connect-timeout 30 --speed-limit 1024 --speed-time 60 \ + -w '%{http_code}' -o "$part" "$url"); do case $code in 403 | 404) rm -f "$part" @@ -49,13 +50,16 @@ fetch_pinned() { return 1 ;; esac - if ! curl -fsSL --connect-timeout 30 --speed-limit 1024 --speed-time 60 \ - --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "$part" "$url"; then + attempt=$((attempt + 1)) + # Seven attempts at most. Like curl's retry-max-time, this bounds when + # another attempt starts; an active transfer keeps its stall limits. + if [ "$attempt" -ge 7 ] || [ "$((SECONDS - started + 30))" -ge 300 ]; then rm -f "$part" echo "ERROR: cannot download ${url}" >&2 return 1 fi - fi + sleep 30 + done got=$(pinned_sha256 "$part") if [ "$got" != "$want" ]; then rm -f "$part" From 030fbdfec608cb0cee5a9722bb7c0572dab0e257 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 13 Nov 2026 14:51:23 +0100 Subject: [PATCH 7/7] fix(ci): stop E2E legs dying at the spare-disk mount on runners that list zram0 Some runner images list /dev/zram0 as a disk of 0 bytes with no file system. The spare-disk step took it for the blank disk, and mkfs.ext4 stopped the leg before any case ran. The step now skips zram, ram and loop devices and any disk of 0 bytes. --- .github/workflows/e2e.yml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/.github/workflows/e2e.yml b/.github/workflows/e2e.yml index e7c02c998..05cf8d80e 100644 --- a/.github/workflows/e2e.yml +++ b/.github/workflows/e2e.yml @@ -500,13 +500,17 @@ jobs: run: | if mountpoint -q /mnt; then echo "/mnt is already a mount"; df -h /mnt; exit 0; fi disk="" - while read -r name type; do + while read -r name type size; do [ "$type" = disk ] || continue + # A real disk: some runner images list /dev/zram0 as a disk of + # 0 bytes, and mkfs.ext4 stops on it. + case $name in /dev/zram*|/dev/ram*|/dev/loop*) continue ;; esac + [ "${size:-0}" -gt 0 ] || continue # Blank: no partition, no file system, nothing mounted. [ "$(lsblk -nro NAME "$name" | wc -l)" -eq 1 ] || continue [ -z "$(lsblk -ndro FSTYPE,MOUNTPOINT "$name" | tr -d ' ')" ] || continue disk=$name - done < <(lsblk -dnpo NAME,TYPE) + done < <(lsblk -dnbpo NAME,TYPE,SIZE) if [ -z "$disk" ]; then echo "No blank disk on this runner"; exit 0; fi echo "Formatting ${disk} ($(lsblk -ndo SIZE "$disk")) and mounting it at /mnt" sudo mkfs.ext4 -q -F -E nodiscard "$disk"