ci(e2e): cache GENCODE, test Locityper, and fail a case that checks nothing

.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/cnvpytor_to_vcf.py

@@ -25,6 +25,11 @@
 on the same chromosome overlaps it reciprocally by at least --min-overlap
 (default 0.5) of the longer of the two. Delly's calls are added before
 CNVpytor's, so a deletion both of them call and Manta misses is added once.
+Only a base record with FILTER PASS (or .) and an allele other than 0 in its
+first sample covers a call: a filtered Manta record (MinGQ and the rest) or
+one genotyped 0/0 or ./. stays in the file but hides nothing. A sites-only
+base covers with every DEL and DUP. A Delly PASS call genotyped 0/0 or ./. is
+not added either, so it covers no CNVpytor call.
 
 Every record says where it came from in INFO SVCALLER (the base's records
 get --base-name, default manta). An added record keeps the caller's ID, END,
@@ -155,6 +160,24 @@ def first_gt(fields):
     return values[i] if i < len(values) and values[i] else "./."
 
 
+def call_holds(fields, n_samples):
+    """True when a record is a call this sample carries.
+
+    A record Manta filtered (MinGQ, MaxDepth and the rest) or genotyped as
+    reference (0/0) or no-call (./.) is not a variant this sample carries;
+    letting it cover a Delly or CNVpytor call would drop that call, and the
+    summary (collect_summary.py) skips a 0/0 record, so the deletion would be
+    lost. The same holds for a Delly PASS record genotyped 0/0 or ./.: it is
+    not added, and so does not cover a CNVpytor call. A sites-only record (no
+    sample column) has no genotype and counts as a call.
+    """
+    if fields[6] not in ("PASS", "."):
+        return False
+    if not n_samples:
+        return True
+    return any(a not in ("0", ".", "") for a in first_gt(fields).replace("|", "/").split("/"))
+
+
 def as_float(value):
     try:
         return float(value)
@@ -214,7 +237,7 @@ def main(argv=None):
     for f in base:
         info = info_dict(f[7])
         kind, sp = sv_type(f, info), span(f, info)
-        if kind in ("DEL", "DUP") and sp:
+        if kind in ("DEL", "DUP") and sp and call_holds(f, n_samples):
             known.add(f[0], kind, sp[0], sp[1])
         f[7] = ("" if f[7] in ("", ".") else f[7] + ";") + "SVCALLER=" + args.base_name
 
@@ -249,7 +272,7 @@ def take(caller, records, keep):
         no_depth = [0]
 
         def delly_keep(f, info):
-            if f[6] != "PASS":
+            if f[6] != "PASS" or not call_holds(f, 1 if len(f) > 9 else 0):
                 return False
             sp = span(f, info)
             fmt = f[8].split(":") if len(f) > 9 else []