fix(archaic): keep strict-mask sites marked so the HG002 Neanderthal share lands within the published range

.github/workflows/e2e.yml

@@ -65,6 +65,7 @@ on:
       - 'nextflow_schema.json'
       - 'versions.env'
       - '.github/workflows/e2e.yml'
+      - '.github/actions/**'
       - 'tests/fixtures/annotsv/class45.sh'
       - 'tests/fixtures/comphet/make_fixture.py'
       - 'tests/fixtures/consumer_export/sites.tsv.gz'
@@ -636,6 +637,7 @@ jobs:
           grep -qx archaic-ancestry-1 "$list" && echo "data_archaic=true" >> "$GITHUB_OUTPUT"
           grep -qx bash-step-08-hla "$list" && echo "data_hla=true" >> "$GITHUB_OUTPUT"
           grep -qx prs-scores-and-run-options-1-strand-sex "$list" && echo "data_pgs=true" >> "$GITHUB_OUTPUT"
+          grep -qxE 'callability-1-catalog|cn-intervals-1-catalog' "$list" && echo "data_catalogs=true" >> "$GITHUB_OUTPUT"
           echo "list=${list}" >> "$GITHUB_OUTPUT"
 
       # Before the base outputs: the base leg pulled the same images.
@@ -691,6 +693,10 @@ jobs:
       # of two IPD-IMGT/HLA releases and GENCODE's genes and exons
       # (bash-step-08-hla, which then runs with the network blocked), and the
       # four real PGS Catalog score files of prs-scores-and-run-options-1.
+      # A third keyed the same way: the pinned downloads of the two table
+      # builders (GENCODE's GTF from EBI among them) that
+      # callability-1-catalog and cn-intervals-1-catalog rebuild their
+      # tables from.
       - name: Data cache keys
         id: data-keys
         run: |
@@ -699,6 +705,7 @@ jobs:
           echo "archaic=e2e-data-1kg-sites-${ARCHAIC_1KG_SITES_MD5}" | tee -a "$GITHUB_OUTPUT"
           echo "hla=$(scripts/ci/e2e-run.sh --data-key hla)" | tee -a "$GITHUB_OUTPUT"
           echo "pgs=$(scripts/ci/e2e-run.sh --data-key pgs)" | tee -a "$GITHUB_OUTPUT"
+          echo "catalogs=$(scripts/ci/e2e-run.sh --data-key catalogs)" | tee -a "$GITHUB_OUTPUT"
 
       - name: Restore cached data, MoChA panel
         id: data-mocha
@@ -732,6 +739,14 @@ jobs:
           path: ${{ env.E2E_DATA_CACHE }}/pgs
           key: ${{ steps.data-keys.outputs.pgs }}
 
+      - name: Restore cached data, table builders' downloads
+        id: data-catalogs
+        if: steps.cases-list.outputs.data_catalogs == 'true'
+        uses: actions/cache/restore@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
+        with:
+          path: ${{ env.E2E_DATA_CACHE }}/catalogs
+          key: ${{ steps.data-keys.outputs.catalogs }}
+
       # The PGS Catalog's synthetic HAPNEST panel (prs-2-ancestry and
       # prs-scores-and-run-options-2, 268 MB), from EBI's FTP server, keyed
       # on the sha256 tests/e2e/lib.sh pins. The cases take it from
@@ -778,15 +793,16 @@ jobs:
           E2E_BASE_OUT: ${{ matrix.job == 'base' && format('{0}/base.tar.zst', runner.temp) || '' }}
           E2E_BASE_IN: ${{ matrix.job != 'base' && format('{0}/e2e-base/base.tar.zst', runner.temp) || '' }}
           # Keep the downloads to the end of the job only when a save step uses them.
-          E2E_DATA_KEEP: ${{ (github.ref == 'refs/heads/main' || inputs.save_images) && ((steps.cases-list.outputs.data_mocha == 'true' && steps.data-mocha.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_archaic == 'true' && steps.data-archaic.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_hla == 'true' && steps.data-hla.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_pgs == 'true' && steps.data-pgs.outputs.cache-hit != 'true')) && 'true' || 'false' }}
+          E2E_DATA_KEEP: ${{ (github.ref == 'refs/heads/main' || inputs.save_images) && ((steps.cases-list.outputs.data_mocha == 'true' && steps.data-mocha.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_archaic == 'true' && steps.data-archaic.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_hla == 'true' && steps.data-hla.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_pgs == 'true' && steps.data-pgs.outputs.cache-hit != 'true') || (steps.cases-list.outputs.data_catalogs == 'true' && steps.data-catalogs.outputs.cache-hit != 'true')) && 'true' || 'false' }}
           HIT_MOCHA: ${{ steps.data-mocha.outputs.cache-hit || 'false' }}
           HIT_ARCHAIC: ${{ steps.data-archaic.outputs.cache-hit || 'false' }}
           HIT_HLA: ${{ steps.data-hla.outputs.cache-hit || 'false' }}
           HIT_PGS: ${{ steps.data-pgs.outputs.cache-hit || 'false' }}
-          DATA_ANY: ${{ steps.cases-list.outputs.data_mocha || steps.cases-list.outputs.data_archaic || steps.cases-list.outputs.data_hla || steps.cases-list.outputs.data_pgs || 'false' }}
+          HIT_CATALOGS: ${{ steps.data-catalogs.outputs.cache-hit || 'false' }}
+          DATA_ANY: ${{ steps.cases-list.outputs.data_mocha || steps.cases-list.outputs.data_archaic || steps.cases-list.outputs.data_hla || steps.cases-list.outputs.data_pgs || steps.cases-list.outputs.data_catalogs || 'false' }}
         run: |
           if [ "$DATA_ANY" = true ]; then
-            echo "Data cache (${JOB}): exact key hit, MoChA panel ${HIT_MOCHA}, 1000 Genomes sites ${HIT_ARCHAIC}, IPD-IMGT/HLA and GENCODE ${HIT_HLA}, PGS scores ${HIT_PGS}" | tee -a "$GITHUB_STEP_SUMMARY"
+            echo "Data cache (${JOB}): exact key hit, MoChA panel ${HIT_MOCHA}, 1000 Genomes sites ${HIT_ARCHAIC}, IPD-IMGT/HLA and GENCODE ${HIT_HLA}, PGS scores ${HIT_PGS}, table builders ${HIT_CATALOGS}" | tee -a "$GITHUB_STEP_SUMMARY"
           fi
           scripts/ci/e2e-run.sh
 
@@ -802,6 +818,7 @@ jobs:
           HIT_ARCHAIC: ${{ steps.data-archaic.outputs.cache-hit }}
           HIT_HLA: ${{ steps.data-hla.outputs.cache-hit }}
           HIT_PGS: ${{ steps.data-pgs.outputs.cache-hit }}
+          HIT_CATALOGS: ${{ steps.data-catalogs.outputs.cache-hit }}
         run: |
           d=$E2E_DATA_CACHE
           if [ "$HIT_MOCHA" != true ] && [ -s "${d}/mocha/manifest.txt" ] \
@@ -812,7 +829,7 @@ jobs:
             echo "archaic=true" >> "$GITHUB_OUTPUT"
           fi
           # Whole and checked (data_ok in e2e-run.sh), or not saved.
-          for n in hla pgs; do
+          for n in hla pgs catalogs; do
             hit=HIT_${n^^}
             if [ "${!hit}" != true ] && [ -d "${d}/${n}" ] && scripts/ci/e2e-run.sh --data-check "$n"; then
               echo "${n}=true" >> "$GITHUB_OUTPUT"
@@ -844,6 +861,12 @@ jobs:
           path: ${{ env.E2E_DATA_CACHE }}/pgs
           key: ${{ steps.data-keys.outputs.pgs }}
 
+      - uses: actions/cache/save@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
+        if: always() && steps.data-save.outputs.catalogs == 'true'
+        with:
+          path: ${{ env.E2E_DATA_CACHE }}/catalogs
+          key: ${{ steps.data-keys.outputs.catalogs }}
+
       # The case leaves the panel there only once its sha256 matched.
       - name: Data cache, HAPNEST panel, what to save
         id: hapnest-save
@@ -936,9 +959,10 @@ jobs:
       # Each group keeps its smallest images up to GROUP_BUDGET_MB of zstd
       # output; the rest are pulled on every run. Four groups of 1500 MB are
       # 5.9 GB, and while main replaces one group the old and the new entry
-      # both count: 7.3 GB. With the five data entries above (about 2.3 GB:
-      # MoChA 0.64, 1000 Genomes sites 1.36, HAPNEST 0.27, hla and GENCODE
-      # 0.03, PGS a few KB) that is 9.6 GB, under the 10 GB limit. An image pinned by digest
+      # both count: 7.3 GB. With the six data entries above (about 2.4 GB:
+      # MoChA 0.64, 1000 Genomes sites 1.36, HAPNEST 0.27, the table
+      # builders' downloads 0.08, hla and GENCODE 0.03, PGS a few KB) that is
+      # 9.7 GB, under the 10 GB limit. An image pinned by digest
       # is never saved: docker load gives it back without its digest, so
       # e2e-run.sh would pull it again anyway.
       - name: Image sizes, and save images for the next run

bin/annotsv_postfilter.py

@@ -61,9 +61,13 @@ class 4 or 5 without being a real rare deletion:
                      record's own read depth, duphold's DHFFC (a deletion
                      needs < 0.7) or DHBFC (a duplication needs > 1.3), step
                      15's thresholds; else, for a Delly call over 1,000 bp,
-                     its RDCN (a deletion needs 0 or 1 copies, a duplication
-                     3 or more, or 2 or more on chrX and chrY, where the sex
-                     is not known here); a CNVpytor call is a read-depth call,
+                     its RDCN, read as bin/cnvpytor_to_vcf.py reads it when
+                     step 05 adds a Delly call: a deletion needs fewer copies
+                     than normal and a duplication more, where normal is 1 on
+                     a male's chrX and chrY outside the PARs (--sex, --par-bed)
+                     and 2 elsewhere; with no --sex either 1 or 2 is normal
+                     there, so a deletion needs 0 or 1 copies and a
+                     duplication 2 or more; a CNVpytor call is a read-depth call,
                      so yes. A Delly call step 05 added (FORMAT GT only)
                      over 1,000 bp is yes: cnvpytor_to_vcf.py set aside the
                      ones whose RDCN does not show the change. "not checked"
@@ -94,14 +98,19 @@ class 4 or 5 without being a real rare deletion:
                      other gene is left.
 
 Usage: annotsv_postfilter.py --class45 CLASS45.tsv [--genes GENES.bed]
-       [--rules RULES.tsv] [--consensus CONSENSUS.vcf.gz] --out CHECKED.tsv
+       [--rules RULES.tsv] [--consensus CONSENSUS.vcf.gz]
+       [--sex male|female] [--par-bed PAR.bed] --out CHECKED.tsv
 """
 import argparse
 import csv
 import gzip
 import os
 import sys
 
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
+import cnvpytor_to_vcf  # noqa: E402  the Delly RDCN rule step 05 applies when it adds a call
+import family_report  # noqa: E402  the PAR intervals
+
 EDGE_BP = 5
 COMMON_AF = 0.01
 CONSENSUS_BP = 1000
@@ -229,8 +238,10 @@ def sv_span(row):
     return start, max(end, start + 1)
 
 
-def depth_support(row, caller):
-    """(yes|no|not checked|., detail) for a deletion or duplication."""
+def depth_support(row, caller, sex="", par=None):
+    """(yes|no|not checked|., detail) for a deletion or duplication. sex:
+    male, female or "" (not given); par: family_report.read_par() of the
+    pseudoautosomal regions, or None."""
     t = sv_type(row)
     if t not in ("DEL", "DUP"):
         return ".", ""
@@ -247,9 +258,18 @@ def depth_support(row, caller):
     except (KeyError, ValueError):
         length = 0
     if rdcn is not None and rdcn >= 0 and length > DELLY_DEPTH_BP:
-        sex_chrom = chrom_key(row.get("SV_chrom", "")).upper() in ("X", "Y")
-        yes = rdcn <= 1 if t == "DEL" else rdcn >= (2 if sex_chrom else 3)
-        return ("yes" if yes else "no"), f"Delly RDCN {rdcn:g}"
+        # cnvpytor_to_vcf.py's rule, so a call reads the same in both places:
+        # the normal copies are 1 on a male's chrX and chrY outside the PARs,
+        # 2 elsewhere, and either one there when the sex is not given.
+        chrom = "chr" + chrom_key(row.get("SV_chrom", "")).upper()
+        try:
+            start = int(row["SV_start"])
+        except (KeyError, ValueError):
+            start = 0
+        cn = str(int(rdcn)) if rdcn.is_integer() else f"{rdcn:g}"
+        no = cnvpytor_to_vcf.delly_depth_contradicts(t, chrom, start, start + length, ["RDCN"], [cn], sex,
+                                                     family_report.in_par(par or {}, chrom, start))
+        return ("no" if no else "yes"), f"Delly RDCN {rdcn:g}"
     if caller == "delly" and rdcn is None and length > DELLY_DEPTH_BP:
         # Step 05 adds a Delly call through cnvpytor_to_vcf.py, which keeps
         # FORMAT GT only, after it has set aside every call over 1,000 bp
@@ -325,10 +345,11 @@ def region_rules(row, rules, transcripts):
     return hits, new_genes, downgrade
 
 
-def check_rows(rows, transcripts=None, rules=None):
+def check_rows(rows, transcripts=None, rules=None, sex="", par=None):
     """The class 4 and 5 rows with the added columns. transcripts None: the
     genes file was not given, so the exon edges and gene aliases are not
-    checked. rules: read_rules() of the rule table, or None."""
+    checked. rules: read_rules() of the rule table, or None. sex and par:
+    as depth_support takes them."""
     out = []
     for r in rows:
         o = {k: r.get(k, "") for k in KEEP}
@@ -349,7 +370,7 @@ def check_rows(rows, transcripts=None, rules=None):
         o["Population_source"] = src
         o["Common"] = "yes" if af is not None and af >= COMMON_AF else "no"
         o["AnnotSV_class"] = r.get("ACMG_class", "")
-        o["Depth_support"], depth_detail = depth_support(r, o["Caller"])
+        o["Depth_support"], depth_detail = depth_support(r, o["Caller"], sex, par)
         hits, new_genes, downgrade = region_rules(r, rules or {}, transcripts)
         if o["Depth_support"] == "no":
             hits.insert(0, ("no_depth_support", f"read depth does not show the change ({depth_detail})"))
@@ -432,8 +453,13 @@ def main(argv=None):
     ap.add_argument("--genes", help="AnnotSV's genes.RefSeq.sorted.bed")
     ap.add_argument("--rules", help="assets/sv_rule_regions.tsv")
     ap.add_argument("--consensus", help="step 22's consensus VCF")
+    ap.add_argument("--sex", choices=("male", "female", ""), default="",
+                    help="the sample's sex, for the normal copies on chrX and chrY (default: not given)")
+    ap.add_argument("--par-bed", help="the pseudoautosomal regions (assets/par_grch38.bed)")
     ap.add_argument("--out", required=True)
     a = ap.parse_args(argv)
+    if a.par_bed and not os.path.isfile(a.par_bed):
+        sys.exit(f"ERROR: no PAR file {a.par_bed}")
     rows = read_tsv(a.class45) if os.path.getsize(a.class45) else []
     if rows and not {"SV_chrom", "SV_start", "SV_end", "SV_type", "ACMG_class"} <= set(rows[0]):
         sys.exit(f"ERROR: {a.class45} lacks one of SV_chrom, SV_start, SV_end, SV_type, ACMG_class")
@@ -447,7 +473,7 @@ def main(argv=None):
         if not os.path.isfile(a.rules):
             sys.exit(f"ERROR: no rule table {a.rules}")
         rules = read_rules(a.rules, {chrom_key(r.get("SV_chrom", "")) for r in rows})
-    checked = check_rows(rows, transcripts, rules)
+    checked = check_rows(rows, transcripts, rules, a.sex, family_report.read_par(a.par_bed))
     cols = KEEP + ADDED
     if a.consensus:
         for r, s in zip(checked, consensus_support(checked, a.consensus)):

bin/archaic_report.py

@@ -32,6 +32,10 @@
              the site stays in IBDmix's input and IBDmix's --mask stops it
              scoring, except that a homozygous mismatch there still counts
              against a segment.
+  strict     (setup.sh --archaic) the lifted GRCh38 rows of one chromosome with
+             every row outside the 1000 Genomes strict mask marked masked in
+             all three archaic columns (the same code 10 to 12), instead of
+             dropped, so IBDmix still scores a mismatch there.
   genotypes  (step 43) IBDmix's genotype file of each chromosome: the three
              archaic columns, SAMPLE (yours, from the VCF) and the stand-ins;
              the mask of each archaic genome (chrN.<genome>.mask.bed, the
@@ -55,6 +59,7 @@
 library only.
 """
 import argparse
+import bisect
 import gzip
 import heapq
 import os
@@ -292,6 +297,41 @@ def flush():
     print(f"{chrom}: {rows} rows", file=sys.stderr)
 
 
+def cmd_strict(a):
+    """The lifted GRCh38 rows of one chromosome (no header, the 13 table
+    columns), once per position, with the calls of every row outside the 1000
+    Genomes strict mask marked masked (MASKED plus the ALT count, 9 stays 9).
+    The strict mask is the same for the three genomes, so all three are
+    marked; step 43 then hands IBDmix those positions as its mask."""
+    iv = read_mask(a.strict, a.chrom)
+    if not iv:
+        die(f"{a.strict} holds no interval on chromosome {a.chrom}: a mask that matches nothing "
+            "would mark every site masked")
+    starts = [s for s, _ in iv]
+    seen = set()
+    rows = marked = 0
+    with open(a.rows) as f, open(a.out, "w") as out:
+        for line in f:
+            c = line.rstrip("\n").split("\t")
+            if len(c) != len(TABLE_HEADER):
+                die(f"{a.rows}: a row has {len(c)} columns, want {len(TABLE_HEADER)}: {line[:80]!r}")
+            pos = int(c[1])
+            if pos in seen:
+                continue
+            seen.add(pos)
+            i = bisect.bisect_right(starts, pos - 1) - 1
+            if not (i >= 0 and pos - 1 < iv[i][1]):
+                hit = False
+                for k in (4, 5, 6):
+                    if c[k] in ("0", "1", "2"):
+                        c[k] = str(int(c[k]) + MASKED)
+                        hit = True
+                marked += hit
+            out.write("\t".join(c) + "\n")
+            rows += 1
+    print(f"{a.chrom}: {rows} rows, {marked} outside the strict mask and marked masked", file=sys.stderr)
+
+
 # setup.sh --archaic stamps its table with the sha256 of everything above this
 # section, so an edit below it does not rebuild the table (about 165 GB of downloads).
 
@@ -677,6 +717,11 @@ def main(argv=None):
     t.add_argument("--exclude", action="append", default=[],
                    help="BED: leave out every site inside these regions (segmental duplications, indel neighbourhoods)")
     t.add_argument("--out", required=True)
+    m = sub.add_parser("strict", help="mark the rows outside the GRCh38 strict mask masked (setup.sh --archaic)")
+    m.add_argument("--chrom", required=True)
+    m.add_argument("--strict", required=True, help="BED: the 1000 Genomes GRCh38 strict mask")
+    m.add_argument("--rows", required=True, help="the lifted rows, tab-separated, no header")
+    m.add_argument("--out", required=True)
     g = sub.add_parser("genotypes", help="IBDmix's genotype files for one sample (step 43)")
     g.add_argument("--table-dir", required=True)
     g.add_argument("--vcf", required=True)
@@ -693,7 +738,7 @@ def main(argv=None):
     a = p.parse_args(argv)
     if a.cmd == "genotypes" and a.pseudo < 1:
         die("--pseudo must be 1 or more", 2)
-    {"table": cmd_table, "genotypes": cmd_genotypes, "report": cmd_report}[a.cmd](a)
+    {"table": cmd_table, "strict": cmd_strict, "genotypes": cmd_genotypes, "report": cmd_report}[a.cmd](a)
 
 
 if __name__ == "__main__":