fix(sv): keep AnnotSV genotype and phase so deletions pair in trans, and fix the archaic evidence summary

.github/workflows/public-evidence.yml

@@ -8,24 +8,36 @@ name: Public evidence
 #                      downloaded once (its md5 checked) and indexed, for
 #                      every chromosome job.
 #   archaic-hg002      setup.sh --archaic and step 43 on HG002's GIAB v4.2.1
-#                      genotypes, one job per autosome (about 150 GB of
-#                      archaic downloads in all; each job fetches its three
-#                      files side by side first, and setup.sh resumes them
-#                      and checks their md5), then archaic-share joins the
-#                      22 into one report per population and prints HG002's
-#                      Neanderthal and Denisovan shares, and how many of
-#                      HG002's positions carry two records, one an SNV.
+#                      genotypes, one job per autosome of the chromosomes
+#                      input (all 22 by default: about 150 GB of archaic
+#                      downloads in all; each job fetches its three files
+#                      side by side first, and setup.sh resumes them and
+#                      checks their md5), then archaic-share joins them into
+#                      one report per population and prints HG002's
+#                      Neanderthal and Denisovan shares, the Neanderthal
+#                      total against the published range when all 22 ran,
+#                      and how many of HG002's positions carry two records,
+#                      one an SNV.
 #                      The job "empty" is the red control: an empty region
 #                      list must make it fail.
 #   ancient-footprint  setup.sh --ancient-panel on the full AADR 1240K files,
 #                      then step 44 on HG002, each with the peak memory of
 #                      its containers (the cgroup's memory.peak), the peak
-#                      disk use and the time.
+#                      disk use and the time. The ancient input off skips it.
 #
 # It runs by hand only (workflow_dispatch). Public repo: the GitHub-hosted
 # runner is free here.
 on:
   workflow_dispatch:
+    inputs:
+      chromosomes:
+        description: "Autosomes for the archaic jobs, separated by spaces or commas (for example 22), or all"
+        type: string
+        default: all
+      ancient:
+        description: "Run the ancient-footprint job too"
+        type: boolean
+        default: true
 
 permissions:
   contents: read
@@ -42,10 +54,30 @@ jobs:
   archaic-sites:
     runs-on: ubuntu-latest
     timeout-minutes: 120
+    outputs:
+      parts: ${{ steps.plan.outputs.parts }}
+      count: ${{ steps.plan.outputs.count }}
     steps:
       - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
         with:
           persist-credentials: false
+      - name: The chromosomes to run
+        id: plan
+        env:
+          CHROMOSOMES: ${{ inputs.chromosomes }}
+        run: |
+          set -euo pipefail
+          read -ra want <<< "$(printf '%s' "$CHROMOSOMES" | tr ',\r\n' '   ')"
+          [ "${want[*]:-all}" != all ] || mapfile -t want < <(seq 1 22)
+          for c in "${want[@]}"; do
+            [[ "$c" =~ ^([1-9]|1[0-9]|2[0-2])$ ]] || { echo "chromosomes: '${c}' is not an autosome 1 to 22"; exit 1; }
+          done
+          list=$(printf '%s\n' "${want[@]}" | sort -nu)
+          # The empty part is the red control and always runs.
+          parts=$(printf '%s\n' empty "$list" | jq -Rnc '[inputs]')
+          echo "parts=${parts}" >> "$GITHUB_OUTPUT"
+          echo "count=$(grep -c . <<< "$list")" >> "$GITHUB_OUTPUT"
+          echo "archaic parts: ${parts}"
       - name: The 1000 Genomes phase 3 sites file and its index
         run: |
           set -euo pipefail
@@ -74,8 +106,8 @@ jobs:
       # Eight at a time, so the other pull requests keep runners.
       max-parallel: 8
       matrix:
-        part: [empty, '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17',
-               '18', '19', '20', '21', '22']
+        # empty, then the autosomes of the chromosomes input (archaic-sites' plan step).
+        part: ${{ fromJSON(needs.archaic-sites.outputs.parts) }}
     steps:
       - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
         with:
@@ -245,7 +277,7 @@ jobs:
           if-no-files-found: error
 
   archaic-share:
-    needs: archaic-hg002
+    needs: [archaic-sites, archaic-hg002]
     runs-on: ubuntu-latest
     timeout-minutes: 20
     steps:
@@ -256,11 +288,27 @@ jobs:
         with:
           pattern: archaic-part-*
           path: parts
-      - name: One report per population, over the 22 autosomes
+      - name: Give a lone part its own folder
+        env:
+          PARTS: ${{ needs.archaic-sites.outputs.parts }}
         run: |
           set -euo pipefail
+          # download-artifact puts a single matching artifact straight into
+          # parts/, without the archaic-part-<chromosome> folder the steps
+          # below read.
+          if ! compgen -G 'parts/archaic-part-*' > /dev/null; then
+            c=$(jq -r '.[1]' <<< "$PARTS")
+            mv parts "lone-${c}" && mkdir parts && mv "lone-${c}" "parts/archaic-part-${c}"
+          fi
+          ls parts
+      - name: One report per population, over the autosomes that ran
+        env:
+          WANT: ${{ needs.archaic-sites.outputs.count }}
+        run: |
+          set -euo pipefail
+          if [ "$WANT" -eq 22 ]; then what="whole genome"; else what="${WANT} of the 22 autosomes"; fi
           {
-            echo "### Step 43 on HG002, whole genome (GIAB v4.2.1 genotypes)"
+            echo "### Step 43 on HG002, ${what} (GIAB v4.2.1 genotypes)"
             echo ""
             echo "| chromosome | archaic download GB | download seconds | setup seconds | table GB | table rows |"
             echo "|---|---|---|---|---|---|"
@@ -281,11 +329,21 @@ jobs:
             python3 bin/archaic_report.py report --sample HG002 --ibdmix-dir "${d}/ibdmix" --spans "${d}/spans.tsv" \
               --pop "$pop" --out-dir "${d}/out"
             kv() { awk -F'\t' -v k="$1" '$1 == k {print $2}' "${d}/out/HG002_archaic.tsv"; }
-            [ "$(kv analysed_chromosomes)" -eq 22 ] || { echo "${pop}: $(kv analysed_chromosomes) chromosomes, want 22"; exit 1; }
+            [ "$(kv analysed_chromosomes)" -eq "$WANT" ] || { echo "${pop}: $(kv analysed_chromosomes) chromosomes, want ${WANT}"; exit 1; }
             [ "$(kv analysed_sites)" -gt 0 ] || { echo "${pop}: no site analysed"; exit 1; }
             nea=$(kv neanderthal_percent)
+            mb=$(awk -v b="$(kv neanderthal_bp)" 'BEGIN {printf "%.1f", b / 1e6}')
+            # Chen et al. 2020 (Cell 180:677), released per-person calls for
+            # Europeans: the union of the Altai and Vindija segments, LOD 4,
+            # 50 kb or longer, the cutoffs of step 43 (docs/43-archaic.md).
+            if [ "$WANT" -eq 22 ]; then
+              where=$(awk -v b="$(kv neanderthal_bp)" 'BEGIN {print (b < 49.9e6 ? "below" : (b > 82.9e6 ? "above" : "inside"))}')
+              cmp="HG002's ${mb} Mb of Neanderthal sequence is ${where} the published range: Chen et al. 2020 released per-person calls with the same method and cutoffs (Altai and Vindija union, LOD 4, 50 kb or longer), and Europeans carry 49.9 to 82.9 Mb in them (see docs/43-archaic.md)."
+            else
+              cmp="HG002 has ${mb} Mb of Neanderthal sequence on these ${WANT} autosomes. The published range of 49.9 to 82.9 Mb (Chen et al. 2020, Europeans, Altai and Vindija union, LOD 4, 50 kb or longer) is for the 22 autosomes, so it is not compared here."
+            fi
             {
-              echo "**1000 Genomes ${pop} frequencies:** Neanderthal ${nea}% ($(kv neanderthal_segments) segments, longest $(kv neanderthal_longest_bp) bp), Denisovan $(kv denisovan_percent)% ($(kv denisovan_segments) segments), over $(kv analysed_chromosomes) autosomes, $(kv analysed_bp) bp, $(kv analysed_sites) sites. Compare this share only with other shares from this step, not with published figures (see docs/43-archaic.md)."
+              echo "**1000 Genomes ${pop} frequencies:** Neanderthal ${nea}% ($(kv neanderthal_segments) segments, longest $(kv neanderthal_longest_bp) bp), Denisovan $(kv denisovan_percent)% ($(kv denisovan_segments) segments), over $(kv analysed_chromosomes) autosomes, $(kv analysed_bp) bp, $(kv analysed_sites) sites. ${cmp}"
               echo ""
             } >> "$GITHUB_STEP_SUMMARY"
           done
@@ -304,6 +362,7 @@ jobs:
           cat "$GITHUB_STEP_SUMMARY"
 
   ancient-footprint:
+    if: inputs.ancient
     runs-on: ubuntu-latest
     timeout-minutes: 350
     steps:

bin/annotsv_postfilter.py

@@ -46,6 +46,12 @@ class 4 or 5 without being a real rare deletion:
                      consensus VCF's sample columns, which step 22 names
                      after the callers), else the SV's own caller
 
+The table ends with AnnotSV's FORMAT column and the sample column after it,
+as the input has them (the caller's GT and, when the caller phased the call,
+PS): the reports pair a deletion with a small variant of the same gene by that
+phase (bin/collect_summary.py two_hits), and a 0/0 call is no deletion. An
+input without FORMAT writes neither.
+
 Some class 4 and 5 calls are common or unreliable for a known reason, and
 the table moves them to class 3. ACMG_class is the class after these checks;
 AnnotSV_class keeps AnnotSV's own. Rule names each check that fired and
@@ -191,13 +197,20 @@ def population(row):
 def sample_format(row):
     """{FORMAT key: value} of the first sample column, the one after FORMAT
     in AnnotSV's header; {} when the row has none."""
+    cols = sample_columns(row)
+    if not cols:
+        return {}
+    return dict(zip((row.get("FORMAT") or "").split(":"), (row.get(cols[1]) or "").split(":")))
+
+
+def sample_columns(row):
+    """["FORMAT", <the sample column>] of AnnotSV's header, the two columns
+    sample_format reads; [] when the row has no sample column."""
     keys = list(row)
     if "FORMAT" not in keys:
-        return {}
+        return []
     i = keys.index("FORMAT") + 1
-    if i >= len(keys) or keys[i] == "Annotation_mode":
-        return {}
-    return dict(zip((row.get("FORMAT") or "").split(":"), (row.get(keys[i]) or "").split(":")))
+    return [] if i >= len(keys) or keys[i] == "Annotation_mode" else keys[i - 1:i + 1]
 
 
 def number(v):
@@ -319,6 +332,7 @@ def check_rows(rows, transcripts=None, rules=None):
     out = []
     for r in rows:
         o = {k: r.get(k, "") for k in KEEP}
+        o.update({k: r.get(k) or "" for k in sample_columns(r)})
         o["Caller"] = info_value(r.get("INFO"), "SVCALLER") or "."
         if sv_type(r) != "DEL":
             o["Exon_edge"], o["Exon_edge_detail"] = ".", ""
@@ -439,6 +453,14 @@ def main(argv=None):
         for r, s in zip(checked, consensus_support(checked, a.consensus)):
             r.update({k: "" if v is None else v for k, v in s.items()})
         cols += SUPPORT
+    # Last. A sample named like a column above would lose its GT and PS, so
+    # that stops the step instead.
+    sample = sample_columns(rows[0]) if rows else []
+    clash = [k for k in sample if k in cols]
+    if clash:
+        sys.exit(f"ERROR: the sample column {clash[0]} has the name of a checked-table column, so its GT and PS "
+                 "would be lost; rename the sample")
+    cols += sample
     tmp = a.out + ".part"
     with open(tmp, "w", newline="") as f:
         w = csv.DictWriter(f, fieldnames=cols, delimiter="\t", lineterminator="\n", extrasaction="ignore")

bin/collect_summary.py

@@ -3794,6 +3794,9 @@ def over(chrom, a, b):
             a, _, b = span.partition("-")
             if x["type"] != "DEL" or x["common"] or not (a.isdigit() and b.isdigit()):
                 continue
+            # A call genotyped 0/0 (or 0|0) is no deletion in this sample: it pairs with nothing.
+            if x.get("gt") and set(x["gt"].replace("|", "/").split("/")) == {"0"}:
+                continue
             chrom = chrom if chrom.startswith("chr") else "chr" + chrom
             # SV_start is the base before the deleted bases (annotsv_postfilter.sv_span):
             # they are SV_start+1 to SV_end, 1-based. The text keeps AnnotSV's coordinates.

docs/39-couple.md

@@ -48,7 +48,7 @@ Only the ClinVar screen is required. Each of the others is used when its file is
 - marked as ACOG too, the genes of the conditions in ACOG Committee Opinions 690 and 691 (2017): cystic fibrosis, spinal muscular atrophy, the hemoglobinopathies, fragile X and the Ashkenazi Jewish panel. All 20 are also in the ACMG list;
 - every other gene of an Orphanet disorder with autosomal recessive or X-linked inheritance and a disease-causing germline association Orphanet assessed (Orphadata of 2026-06-23, CC BY 4.0). These rows cover the recessive genes of a ClinVar hit outside the ACMG list. A gene on chrX outside the pseudoautosomal regions is X-linked, every other gene autosomal recessive.
 
-The `severity` column names the alleles of a gene that are not a full loss of function: a hypomorph (BTD p.Asp444His, the GALT Duarte variants) or a pseudodeficiency allele (HEXA c.739C>T and c.745C>T). Such an allele is never a carrier finding on its own; the couple report lists it under what the screen does not see, and the sample report's own carrier list sets it aside.
+The `severity` column names the alleles of a gene that are not a full loss of function: a hypomorph (BTD p.Asp444His, the GALT Duarte variants) or a pseudodeficiency allele (HEXA c.739C>T and c.745C>T). Such an allele is never a carrier finding on its own; the couple report lists it under what the screen does not see, and the sample report's own carrier list sets it aside. The couple screen does count a hypomorph beside the other partner's full loss of function in the same gene: an autosomal recessive gene is then `shared`, because a child who inherits both can have the condition, for some genes in a milder form. An X-linked gene keeps the `x_linked` result. For example, an RBM8A low-expression allele in one partner and, in the other, a deletion that removes RBM8A (the 1q21.1 deletion) is the pair behind TAR syndrome, and the report names the allele as a hypomorph in that row. A hypomorph in both partners, or in one alone, stays set aside.
 
 A ClinVar hit in a gene outside the list (a dominant gene, for example) is not part of this screen; each partner's own ClinVar section lists it. An emptied or damaged list stops the step instead of reporting no shared gene.
 

docs/testing.md

@@ -155,7 +155,7 @@ This is every step script that no e2e case runs, plus the full-size index build
 | Step | Why no e2e case runs it | What checks it instead |
 |---|---|---|
 | 26b (fine-grained ancestry) on the full panel | the install reads about 19 GB from gnomAD's bucket | case `ancestry-fine-grained-1` runs it on a reduced panel of the fixture's regions plus chr22:20-35 Mb |
-| 43 (archaic ancestry) on the whole genome | `setup.sh --archaic` reads 150 GB of archaic genomes | case `archaic-ancestry-1` runs it on a slice; the `archaic-hg002` jobs of `public-evidence.yml` run setup and step 43 on HG002's 22 autosomes, by hand |
+| 43 (archaic ancestry) on the whole genome | `setup.sh --archaic` reads 150 GB of archaic genomes | case `archaic-ancestry-1` runs it on the whole of chr22; the `archaic-hg002` jobs of `public-evidence.yml` run setup and step 43 on HG002's 22 autosomes, by hand |
 | 44 (ancient populations) on the full AADR panel | the install reads the 7.1 GB AADR release and runs smartpca on it for an hour | case `ancient-similarity-1` runs it on a planted panel; the `ancient-footprint` job of `public-evidence.yml` runs setup and step 44 on HG002, by hand |
 | 13 (VEP, offline) | the cache is 26 GB; the fixture's VEP subset was annotated once with `--database` instead | the image's smoke test, which annotates 50 records with `--database` |
 | 17 (CPSR) | the PCGR bundle (about 8 GB) plus a second VEP cache. CI never runs real CPSR | the image's smoke test, and the fake-docker cases for the command line |
@@ -182,10 +182,14 @@ The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub`
 
 ## The public-evidence workflow
 
-`.github/workflows/public-evidence.yml` runs two opt-in steps whole, on public data only, and writes what they measured to the job summary. It is not a gate; it runs by hand only (`workflow_dispatch`).
+`.github/workflows/public-evidence.yml` runs two opt-in steps whole, on public data only, and writes what they measured to the job summary. It is not a gate; it runs by hand only (`workflow_dispatch`). Two inputs pick what runs: `chromosomes` (default `all`) names the autosomes for the archaic jobs, for example `22` for one chromosome job (plus the `empty` red control) and that chromosome's three archaic files instead of 22 chromosome jobs and 150 GB, and `ancient` (default on) runs the ancient-footprint job:
+
+```bash
+gh workflow run public-evidence.yml -f chromosomes=22 -f ancient=false
+```
 
 - `ancient-footprint`: `setup.sh --ancient-panel` on the full AADR release, then step 44 on a gVCF built from GIAB's HG002 v4.2.1 genotypes, with the peak memory of each step's containers (one maximum per step, not per container), the peak disk growth and the time. It fails if a command fails, if the sampler saw no container memory, or if step 44 used no site. What it measured is in [step 44, Runtime](44-ancient.md#runtime).
-- `archaic-hg002`: `setup.sh --archaic` and step 43 on HG002, one job per autosome; `archaic-share` joins them into one report per population. A job fails if a command fails or step 43 analysed no site, and the job `empty` is its red control: an empty region list must make setup fail. What it measured is in [step 43, Runtime](43-archaic.md#runtime).
+- `archaic-hg002`: `setup.sh --archaic` and step 43 on HG002, one job per autosome; `archaic-share` joins them into one report per population. When all 22 autosomes ran, the summary says whether HG002's Neanderthal total lies inside the 49.9 to 82.9 Mb that Chen et al.'s released per-person calls give Europeans with the same method and cutoffs; it prints the comparison and does not gate on it. A job fails if a command fails or step 43 analysed no site, and the job `empty` is its red control: an empty region list must make setup fail. What it measured is in [step 43, Runtime](43-archaic.md#runtime).
 
 ## The prs-memory job
 

tests/test_annotsv_postfilter.py

@@ -36,7 +36,9 @@
      "not checked" and keeps its class; a CNVpytor call has support; a
      Delly call with FORMAT GT only (as step 05 adds it) has support over
      1,000 bp and is "not checked" at 800 bp; the
-     summary line counts the rows moved to class 3;
+     summary line counts the rows moved to class 3; every written row keeps
+     AnnotSV's FORMAT and sample columns, and a sample named like a checked
+     column (Caller) stops the step instead of losing its GT;
   8. red controls: without --rules the DSPP and IGH deletions stay class 4;
      and the rule table itself: every row has a known rule, IGH lies on
      chr14 between 105.5 and 107 Mb, the MHC is GRCh38's 28,510,120 to
@@ -205,6 +207,16 @@ def col(i, c):
     want = sorted(["TR_IN", "IGH_IN", "IGH_MOST", "MHC_IN", "ALIAS", "D_LOW", "DUP_LOW", "DELLY2"])
     check("exactly the eight rows above are moved to class 3", moved == want, moved)
     check("the summary line counts them", "moved to class 3: 8; depth support no: 3, not checked: 2" in err, err)
+    check("every row keeps AnnotSV's FORMAT and sample column (GT and the depth fields)",
+          (col("D_OK", "FORMAT"), col("D_OK", "S1"), col("NODEPTH", "S1")) == ("GT:DHFFC", "0/1:0.5", "0/1"),
+          got.get("D_OK"))
+    clash = os.path.join(work, "clash.tsv")
+    with open(clash, "w") as f:
+        f.write("\t".join("Caller" if h == "S1" else h for h in HEAD7) + "\n" + ROWS7)
+    p = subprocess.run([sys.executable, os.path.join(REPO, "bin", "annotsv_postfilter.py"), "--class45", clash,
+                        "--out", os.path.join(work, "clash_checked.tsv")], capture_output=True, text=True)
+    check("a sample named Caller stops the step instead of losing its GT and PS",
+          p.returncode != 0 and "rename the sample" in p.stderr, (p.returncode, p.stderr[-300:]))
 
     # 8. red controls
     plain, _ = run7(work, "norules", with_rules=False)