fix(cn): one SMN1 exon-7 rule and raw MAPQ0 shares for the LPA null allele

.github/workflows/guard.yml

@@ -313,7 +313,7 @@ jobs:
           grep -q 'summary.json validates (og dev schema)' <<< "$out" \
             || { echo "ERROR: reports-stale-step did not run og dev schema"; exit 1; }
 
-      # The tests of the og filters steps 04, 05 and 07 run, og chip's (step 49), og family's (step 38), og run delta's, step 26b's og ancestry-fine,
+      # The tests of the og filters steps 04, 05 and 07 run, og chip's (step 49), og cn's (steps 35, 41 and 46), og family's (step 38), og run delta's, step 26b's og ancestry-fine,
       # the opt-in step reports' (45, 47, 48, 50 to 53), the
       # og tools setup.sh runs (og str catalog, og setup unzip and genetic-map), and the step 25
       # tests that compare og prs format with the folders they write by hand:
@@ -334,7 +334,9 @@ jobs:
                    tests/test_family_report.py tests/test_family_segments.py \
                    tests/test_vntyper_report.py tests/test_locityper_report.py tests/test_mocha_report.py \
                    tests/test_mei_report.py tests/test_viral_report.py tests/test_ystr_report.py \
-                   tests/test_exomiser_report.py; do
+                   tests/test_exomiser_report.py \
+                   tests/test_copy_number.py tests/test_hba_copy_number.py \
+                   tests/test_lpa_kiv2.py tests/test_traits_report.py tests/test_callability.py; do
             rc=0
             echo "::group::${t}"
             case "$t" in
@@ -346,6 +348,15 @@ jobs:
           done
           exit "$fail"
 
+      # The fake-docker job builds no og, so its cn-engine-calls case says it
+      # is skipped there; here bin/og is built and the case must run.
+      - name: The fake-docker case cn-engine-calls (steps 35, 41 and 46 through og cn) with og
+        working-directory: .
+        run: |
+          out=$(scripts/ci/run-fake-docker-suite.sh cn-engine-calls 2>&1) || { echo "$out"; exit 1; }
+          echo "$out"
+          if grep -q '^skip: bin/og is not built' <<<"$out"; then echo "ERROR: the case was skipped, with og built"; exit 1; fi
+
       # A Python file kept only for what other Python imports from it holds
       # those names and nothing else (tests/parity/remnants.tsv).
       - name: og dev remnants

.gitleaksignore

@@ -4,3 +4,8 @@
 # A test string in the Parascopy freshness self-test ("GRCh37_v9.9.tar.gz"),
 # not a key; the line was rewritten in 7d6dad4.
 63e12d6c2b29078a49e14de8ef4c4e57a9ebd8b9:scripts/ci/freshness.py:generic-api-key:940
+# A Rust constant name compared with a marker's name in og-cn, not a key; the
+# comparison was turned round in the same pull request. The second line is this
+# file's own first comment about it, which quoted the expression.
+0bc60efda244e5361c54fe411f46db9514e5e8c3:rust/crates/og-cn/src/cli.rs:generic-api-key:590
+63b5a9273c444970e469d04c8da52aa8890726b5:.gitleaksignore:generic-api-key:7

CHANGELOG.md

@@ -4,6 +4,17 @@ Each release's notes. From 1.0.0 on, a major or minor release has hand-written n
 
 ## Unreleased
 
+- The copy-number engine now runs in the Rust `og` binary, as `og cn bed`, `og cn call`, `og cn intervals`, `og cn check-bam`, `og cn marker-regions` and `og cn marker-read`. That covers step 35's depth counts (HBA, the gene-dosage rows, CYP21A2, the interval rows) and its SMN1 marker, step 46's KIV-2 count and 4925G>A reading, the RHD regions of step 41 and the ALT-contig check of steps 32 and 35. They replace `bin/hba_copy_number.py` and `bin/lpa_kiv2.py`, which are gone, the commands of `bin/copy_number.py` and the marker reader of `bin/traits_report.py`. Two parts stay in Python for now: step 41's report makes the RhD call with the catalog reader and `assess()` of `bin/copy_number.py`, and step 32's CYP2D6 depth check is its own tool. `bin/traits_report.py` keeps `MARKERS` for the summary. Every table, file name and option a step uses stays the same; the parity cases in `tests/parity` hold every output that changed on purpose:
+  - Step 46: the KIV-2 4925G>A reading is made in the same call as the count and judges the aligner on the count's own MAPQ 0 share. Before, it read the share the count table prints, rounded to three places, so at a share just under 0.20 the count said not assessed while the marker still read `present`. Both now say not assessed.
+  - Step 35: the SMN1 marker message takes the SMN1 psCN from the first row of the SMN table that covers exon 7, the rule of the Paraphase comparison. It took the last such row, so a table with two rows over exon 7 made the two read different counts.
+  - Steps 35 and 46: a status the HBA or LPA table layout did not name (a catalog edit could give one) now writes `not assessed: <reason>`. The HBA table wrote `called` with no count, and the LPA table stopped with an error.
+  - A catalog row with a `more=` reading and no numbered reading is refused when the catalog loads, naming the row. Before, a call that reached it failed without naming the row.
+  - `og cn call --sex` reads a gene row on chrX or chrY against the sample's count there. The Python call had no sex, so such a row could only read `sex_unknown`. No row of the shipped catalog is a gene row on chrX or chrY.
+  - The tables are written to a temporary file and renamed, so a failed run leaves no half-written table.
+  - `og cn bed` and `og cn marker-regions` exit 2 when they cannot write the region list, so a step never reads a truncated list.
+  - `copy_number.py call --json` is gone with no replacement: no step passed it.
+  - Error messages about a wrong regions table say to make it with `og cn bed`.
+
 - The opt-in steps 45, 47, 48 and 50 to 53 now read their tool's output with the Rust `og` binary: `og y str-bed` and `og y str-report` (step 45), `og viral report` (47), `og exomiser prepare` and `og exomiser report` (48), `og mei report` and `og mei annotsv-input` (50), `og mocha report` (51), `og vntyper report` (52) and `og locityper report` (53). They replace `bin/ystr_report.py`, `bin/viral_report.py`, `bin/exomiser_report.py`, `bin/mei_report.py`, `bin/mocha_report.py`, `bin/vntyper_report.py` and `bin/locityper_report.py`, which are gone. og is static, so each step still runs in the image it ran in. The tables are the same, byte for byte, on the 123 parity cases, except where a line below says otherwise. Every output is written whole or not at all, and an input a step cannot read stops it with one `ERROR:` line, not a Python traceback.
   - Step 48: the summary printed at the end of `og exomiser report` has one line per rank. The Python cut it by variant count, so a rank with several variants could push the closing note in or out. The text file is unchanged.
   - Step 50: the depth of each catalog window is its own. The Python looked it up by gene and chromosome, so two windows of one gene on one chromosome both got the last one's depth. The shipped catalog has no such pair, so its results do not change.

bin/bloodtype_report.py

@@ -26,8 +26,8 @@
 
 --depth is samtools bedcov -Q 20 -j over the table's sites and the RHD row of
 bin/cn_loci.tsv (scripts/41-bloodtype.sh and the BLOODTYPE_DEPTH process make
-it from the BAM; `copy_number.py bed --locus RHD` prints that row's regions,
-and bin/copy_number.py turns their depth into the copy number). It gives:
+it from the BAM; `og cn bed --locus RHD` prints that row's regions, and the
+catalog code of bin/copy_number.py turns their depth into the copy number). It gives:
   - RhD. The common RhD-negative allele (RHD*01N.01) deletes the whole RHD gene,
     which a VCF does not show. RHD depth over the depth of its two flanks gives
     the copy number (0, 1, 2); exons 3 to 7 without RHD reads in a gene that is

bin/collect_summary.py

@@ -128,7 +128,7 @@
     "Variants of uncertain significance (not interpreted)",
 ]
 
-# Added to NOT_ASSESSED unless step 35's depth count (bin/hba_copy_number.py,
+# Added to NOT_ASSESSED unless step 35's depth count (og cn call --locus HBA,
 # short reads) or step 35b (Paraphase, PacBio HiFi, both genes) made a call.
 HBA_NOT_ASSESSED = ("HBA1/HBA2 copy number (alpha-thalassaemia deletions): neither step 35's read-depth count "
                     "(short reads) nor step 35b (Paraphase, PacBio HiFi) made a call")
@@ -910,13 +910,13 @@ def sec_parascopy(d, s):
                "smn1_marker": read_marker_table(os.path.join(d, "paralogs", f"{s}_smn1_marker.tsv"), "35",
                                                 "SMN1_27134"),
                "rows": len(rows),
-               # HBA1/HBA2 from depth (bin/hba_copy_number.py); None when
+               # HBA1/HBA2 from depth (og cn call --locus HBA); None when
                # step 35 wrote no HBA table (a run from before the check)
                "hba": read_hba(hp) if os.path.isfile(hp) else None,
                # the interval rows of bin/cn_loci.tsv (recurrent CNV regions,
                # carrier deletions, AZF); None without step 35's table
                "cn_intervals": read_cn_intervals(ip) if os.path.isfile(ip) else None,
-               # CYP21A2 from depth and PSV sites (bin/copy_number.py, the
+               # CYP21A2 from depth and PSV sites (og cn call, the
                # CYP21A2 row): copies, module total, finding; None without it
                "cyp21a2": read_hba(cp) if os.path.isfile(cp) else None,
                # the gene-dosage rows (HP, LCE3BC, CCL3L1): total copies from
@@ -931,7 +931,7 @@ def sec_parascopy(d, s):
 
 
 def read_cn_intervals(path):
-    """Step 35's interval table (bin/copy_number.py intervals): the rows with
+    """Step 35's interval table (og cn intervals): the rows with
     a deletion or a duplication, and a count of the rows by status."""
     rows = read_tsv(path)
     cols = ("locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change")
@@ -2523,7 +2523,7 @@ def sec_archaic(d, s):
 
 
 def sec_lpa(d, s):
-    """Step 46's KIV-2 count (lpa/<id>_lpa_kiv2.tsv, bin/lpa_kiv2.py), with
+    """Step 46's KIV-2 count (lpa/<id>_lpa_kiv2.tsv, og cn call --locus LPA_KIV2), with
     step 42's two LPA risk variants when the traits table is there, and
     whether either points to a genetically high Lp(a)."""
     p = first_existing(d, [f"lpa/{s}_lpa_kiv2.tsv"])