feat(summary): copy-number, traits, SV, ancestry and the other sections in Rust

docs/rust.md

@@ -21,9 +21,7 @@ rust/
   crates/og-qc/         og qc sex-loss: mosaic loss of Y or X for step 33
   crates/og-setup/      what setup.sh runs once: the STRchive catalog, zips, the genetic map
   crates/og-steps/      the step helpers: constraint-join, prs format, step 40's consumer export; og family (step 38); og pgx outside-calls|hla-concordance|pypgx-summary (steps 08, 32, 36); og paraphase (step 35b)
-  crates/og-summary/    the step tables of steps 20, 25 and 33: og mito, og prs table, og qc sample;
-                        og summary collect --only a, the first seven summary sections (no caller yet);
-                        --only b, the clinical sections: ClinVar, carrier, PGx, CPSR, ACMG SF, two hits
+  crates/og-summary/    the step tables of steps 20, 25 and 33: og mito, og prs table, og qc sample; og summary collect --only a|b|c, QC, clinical and the remaining summary sections (no caller yet); og roh summary (step 11)
   crates/og-pedigree/   the pedigree rules: families files, members, kinship, PARs
   crates/og-chip/       step 49, clonal haematopoiesis: og chip prepare, report, check-table
   crates/og-cn/         copy number from depth (steps 32, 35, 41, 46): og cn bed, call, intervals, check-bam, marker-regions, marker-read

modules/local/roh/main.nf

@@ -45,24 +45,11 @@ process ROH {
         | bcftools roh \${ROH_FLAGS} -o ${meta.id}_roh.txt -
 
     # Summary: autosomal segments with 5 Mb or more outside GRCh38's
-    # centromere and heterochromatin bands (the acen, gvar and stalk rows of
-    # UCSC's cytoBand.hg38.txt, merged; the line of scripts/11-roh-analysis.sh
-    # and bin/collect_summary.py), the threshold of docs/11-roh-analysis.md;
-    # length_bp and length_mb are of that part. The bash step writes the same
-    # file.
-    MASKED_BANDS="chr1:121700000-143200000 chr2:91800000-96000000 chr3:87800000-98600000 chr4:48200000-51800000 chr5:46100000-51400000 chr6:58500000-62600000 chr7:58100000-62100000 chr8:43200000-47200000 chr9:42200000-61500000 chr10:38000000-41600000 chr11:51000000-55800000 chr12:33200000-37800000 chr13:0-18900000 chr14:0-18200000 chr15:0-20500000 chr16:35300000-47000000 chr17:22700000-27400000 chr18:15400000-21500000 chr19:19900000-31900000 chr20:25700000-30400000 chr21:0-13000000 chr22:0-17400000 chrX:58100000-63800000 chrY:10300000-10600000 chrY:26600000-57227415"
-    echo "# ROH Summary for ${meta.id}" > ${meta.id}_roh_summary.txt
-    echo "# Segments >=5MB on autosomes (excludes chrX/chrY)" >> ${meta.id}_roh_summary.txt
-    printf 'chrom\\tstart\\tend\\tlength_bp\\tlength_mb\\n' >> ${meta.id}_roh_summary.txt
-    awk -v bands="\$MASKED_BANDS" '
-        BEGIN { n = split(bands, b, " "); for (i = 1; i <= n; i++) { split(b[i], x, /[:-]/); c[i] = x[1]; s[i] = x[2]; e[i] = x[3] } }
-        \$1 == "RG" && \$3 ~ /^chr([1-9]|1[0-9]|2[0-2])\$/ {
-          m = 0
-          for (i = 1; i <= n; i++) if (c[i] == \$3) { o = (\$5 < e[i] ? \$5 : e[i]) - (\$4 - 1 > s[i] ? \$4 - 1 : s[i]); if (o > 0) m += o }
-          if (m > \$6) m = \$6
-          l = \$6 - m
-          if (l >= 5000000) printf "%s\\t%s\\t%s\\t%d\\t%.1f\\n", \$3, \$4, \$5, l, l / 1e6
-        }' ${meta.id}_roh.txt >> ${meta.id}_roh_summary.txt
+    # centromere and heterochromatin bands (MASKED_BANDS in og's step 11
+    # code, the one copy), the threshold of docs/11-roh-analysis.md;
+    # length_bp and length_mb are of that part. The bash step runs the same
+    # command.
+    og roh summary --sample ${meta.id} --roh ${meta.id}_roh.txt --out ${meta.id}_roh_summary.txt
 
     cat <<-END_VERSIONS > versions.yml
     "${task.process}":

rust/README.md

@@ -10,7 +10,7 @@ OwnGenome's own tools, ported from `bin/*.py` one at a time into one static bina
 | `crates/og-bam` | BAM reading as pysam does it (fetch, `reference_end`, `cigarstring`, `get_blocks`), on htslib linked statically, and the tools that read BAMs: the reads at a variant (`reads`: one CIGAR walk per read, SA pieces, filters per tool), `og cdna check`, `og cdna scan`, `og comphet phase`, `og bam contigs`, the CYP2D6 depth check and the RFC1 motif count |
 | `crates/og-pedigree` | the pedigree rules of step 38: families files, members, kinship, declared relationships, PARs |
 | `crates/og-qc` | `og qc sex-loss`: step 33's mosaic loss of Y or X, with CPython's medians, `fmean` and exact `variance` |
-| `crates/og-summary` | the step tables: `og qc sample` (step 33), `og prs table` (step 25), `og mito annotate`, `og mito deletions`, `og mito copies` (step 20); the sections of summary.json, `og summary collect --only a` (QC) and `--only b` (ClinVar, carrier list, PGx, CPSR, the clinical filter with ACMG SF and two hits, slivar), beside the Python collector, which still writes it |
+| `crates/og-summary` | the step tables: `og qc sample` (step 33), `og prs table` (step 25), `og mito annotate`, `og mito deletions`, `og mito copies` (step 20); the sections of summary.json, `og summary collect --only a` (QC), `--only b` (ClinVar, carrier list, PGx, CPSR, the clinical filter with ACMG SF and two hits, slivar) and `--only c` (the remaining summary sections), beside the Python collector, which still writes summary.json; `og roh summary` (step 11) |
 | `crates/og-chip` | step 49, clonal haematopoiesis: `og chip prepare`, `og chip report`, `og chip check-table` |
 | `crates/og-cn` | copy number from depth, steps 32, 35, 41 and 46: `og cn bed`, `og cn call`, `og cn intervals`, `og cn check-bam`, `og cn marker-regions`, `og cn marker-read` |
 | `crates/og-bloodtype` | step 41, blood groups: `og bloodtype report`, `og bloodtype check-markers`, `og bloodtype systems` |

rust/crates/og-summary/Cargo.toml

@@ -12,6 +12,9 @@ og-io.workspace = true
 og-bam.workspace = true
 og-clinvar.workspace = true
 og-steps.workspace = true
+og-ancestry.workspace = true
+og-chip.workspace = true
+og-cn.workspace = true
 indexmap.workspace = true
 memchr.workspace = true
 

rust/crates/og-summary/src/collect.rs

@@ -10,6 +10,10 @@
 //! section's state. Group `b` is the clinical sections ([`group_b`]):
 //! clinvar, carrier, pharmcat, cpic, pypgx, cyrius, cpsr, clinical (with
 //! the ACMG SF tier and the two-hit rule) and slivar, and the CYP2D6 block.
+//! Group `c` is copy number, HLA, the polygenic scores,
+//! ancestry, structural variants, repeats, telomeres, runs of
+//! homozygosity, the haplogroups, the opt-in step tables, chrM, traits with
+//! the APOL1 G2 deletion, and the couple screen ([`group_c`]).
 //!
 //! How it reads a folder: one walk of the folder ([`scan::FolderIndex`])
 //! instead of a `glob` per candidate path; each input read once (the
@@ -25,6 +29,7 @@ pub mod clinical;
 pub mod clinvar;
 pub mod cpsr;
 pub mod group_b;
+pub mod group_c;
 pub mod model;
 pub mod pgx;
 pub mod scan;
@@ -41,7 +46,8 @@ use sections::{Callability, Coverage, Data, Mosaic, ReadCount, SampleQc, SexChec
 use std::io::Write;
 use std::path::Path;
 
-const USAGE: &str = "usage: og summary collect --sample SAMPLE --sample-dir SAMPLE_DIR --only {a,b}
+const USAGE: &str =
+    "usage: og summary collect --sample SAMPLE --sample-dir SAMPLE_DIR --only {a,b,c}
                           --out OUT [--declared-sex DECLARED_SEX] [--assets ASSETS]
 ";
 
@@ -51,7 +57,7 @@ const OPTS: &[Opt] = &[
     Opt("--declared-sex", Kind::Value, false),
     Opt("--out", Kind::Value, true),
     Opt("--assets", Kind::Value, false),
-    Opt("--only", Kind::Choice(&["a", "b"]), true),
+    Opt("--only", Kind::Choice(&["a", "b", "c"]), true),
 ];
 
 /// What one collect reads.
@@ -64,6 +70,18 @@ pub struct Input {
     pub declared_sex: Option<String>,
     /// The repository's `assets/` folder.
     pub assets: String,
+    /// The group of sections to read: `a`, `b` or `c`.
+    pub group: Group,
+    /// The folder of the two tables group c reads from `bin/`.
+    pub tables: String,
+}
+
+/// The group of sections a collect reads.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Group {
+    A,
+    B,
+    C,
 }
 
 /// `og summary collect`: exit 0 when summary.json is written (whatever
@@ -83,6 +101,12 @@ pub fn cli(argv: &[String], out: &mut impl Write, err: &mut impl Write) -> u8 {
         sample_dir: a.text("--sample-dir").to_owned(),
         declared_sex: a.get("--declared-sex").map(str::to_owned),
         assets: a.get("--assets").map_or_else(default_assets, str::to_owned),
+        tables: group_c::util::default_tables(),
+        group: match a.text("--only") {
+            "b" => Group::B,
+            "c" => Group::C,
+            _ => Group::A,
+        },
     };
     // A part of the sections must never replace the sample's complete
     // summary.json, so there is no default path while `--only` exists.
@@ -149,6 +173,34 @@ pub fn collect(input: &Input) -> Result<Summary, String> {
     let s = pycompat::mtime::escape(&input.sample);
     let find = |pats: &[String]| index.first(pats).map(str::to_owned);
     let file = |rel: String| -> (String, String) { (crate::join(&input.sample_dir, &rel), rel) };
+    if input.group == Group::C {
+        // The declared sex the chrY filter of CNVpytor reads comes from
+        // step 33's table when the run status and manifest have none.
+        let qc = find(&[format!("qc/{s}_sample_qc.tsv")]).map(file);
+        let qc_declared = qc.and_then(|(p, _)| {
+            let rows = sections::read_tsv(&p).ok()?;
+            SampleQc::read(&p, &rows)
+                .ok()?
+                .declared_sex()
+                .map(str::to_owned)
+        });
+        let declared_sex = declared(
+            input,
+            run.as_ref(),
+            manifest.as_ref(),
+            qc_declared.as_deref(),
+        );
+        let sections = read_group_c(input, &index, run.as_ref(), declared_sex.as_deref());
+        return Ok(Summary {
+            sample: input.sample.clone(),
+            generated_utc: time::utc_minute(now()).unwrap_or_default(),
+            run,
+            manifest,
+            declared_sex,
+            sections,
+            cyp2d6: None,
+        });
+    }
 
     let variants_at = find(&[format!("vcf/{s}.vcf.gz")]);
     let (variants, rest) = std::thread::scope(|scope| {
@@ -194,25 +246,10 @@ pub fn collect(input: &Input) -> Result<Summary, String> {
             Some(Data::SampleQc(q)) => q.declared_sex(),
             _ => None,
         });
-    let declared_sex = [
-        input.declared_sex.as_deref(),
-        run.as_ref().and_then(|r| r.declared_sex.as_deref()),
-        manifest
-            .as_ref()
-            .and_then(|m| m.value("run", "declared_sex")),
-        qc_declared,
-    ]
-    .into_iter()
-    .flatten()
-    .find(|v| !v.is_empty())
-    .map(str::to_owned);
-
-    let now = std::time::SystemTime::now()
-        .duration_since(std::time::UNIX_EPOCH)
-        .map_or(0.0, |d| d.as_secs_f64());
+    let declared_sex = declared(input, run.as_ref(), manifest.as_ref(), qc_declared);
     Ok(Summary {
         sample: input.sample.clone(),
-        generated_utc: time::utc_minute(now).unwrap_or_default(),
+        generated_utc: time::utc_minute(now()).unwrap_or_default(),
         run,
         manifest,
         declared_sex,
@@ -221,6 +258,67 @@ pub fn collect(input: &Input) -> Result<Summary, String> {
     })
 }
 
+fn now() -> f64 {
+    std::time::SystemTime::now()
+        .duration_since(std::time::UNIX_EPOCH)
+        .map_or(0.0, |d| d.as_secs_f64())
+}
+
+/// The declared sex: the command line's, the run status's, the
+/// manifest's, then step 33's; the first that is not empty.
+fn declared(
+    input: &Input,
+    run: Option<&RunStatus>,
+    manifest: Option<&Manifest>,
+    qc: Option<&str>,
+) -> Option<String> {
+    [
+        input.declared_sex.as_deref(),
+        run.and_then(|r| r.declared_sex.as_deref()),
+        manifest.and_then(|m| m.value("run", "declared_sex")),
+        qc,
+    ]
+    .into_iter()
+    .flatten()
+    .find(|v| !v.is_empty())
+    .map(str::to_owned)
+}
+
+/// Reads the sections of group c, then runs the parts of the post pass
+/// that change them.
+fn read_group_c(
+    input: &Input,
+    index: &FolderIndex,
+    run: Option<&RunStatus>,
+    declared_sex: Option<&str>,
+) -> Vec<Section> {
+    let c = group_c::Ctx {
+        sample: &input.sample,
+        dir: &input.sample_dir,
+        assets: &input.assets,
+        tables: &input.tables,
+        index,
+    };
+    let reads = group_c::read_all(&c);
+    let mut sections: Vec<Section> = group_c::GROUP_C
+        .iter()
+        .zip(reads)
+        .map(|(spec, r)| {
+            let r = r.and_then(|found| match found {
+                None => Ok(None),
+                Some((f, v)) => with_mtime(&f.path, f.rel, Data::Json(v)),
+            });
+            Section::new(*spec, r, run)
+        })
+        .collect();
+    let mut view: Vec<(Spec, model::State, Option<&mut Data>)> = sections
+        .iter_mut()
+        .map(|s| (s.spec, s.state, s.data.as_mut()))
+        .collect();
+    group_c::post(&c, &mut view, &group_c::Others::default(), declared_sex);
+    sections
+}
+
 /// What one reader gave: the file (relative), its mtime and its data.
 type Read = Result<Option<(String, f64, Data)>, String>;