diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 8d4318ab..055335c6 100644 --- a/.github/workflows/public-evidence.yml +++ b/.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- 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: diff --git a/bin/annotsv_postfilter.py b/bin/annotsv_postfilter.py index e96a9466..f09905a7 100755 --- a/bin/annotsv_postfilter.py +++ b/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", ] 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") diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 514cce4c..225a40a9 100755 --- a/bin/collect_summary.py +++ b/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. diff --git a/docs/39-couple.md b/docs/39-couple.md index dbe6e0cf..019260bc 100644 --- a/docs/39-couple.md +++ b/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. diff --git a/docs/testing.md b/docs/testing.md index 30ab07e4..558dbc6b 100644 --- a/docs/testing.md +++ b/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 diff --git a/tests/test_annotsv_postfilter.py b/tests/test_annotsv_postfilter.py index 8592596c..0e84c9b2 100755 --- a/tests/test_annotsv_postfilter.py +++ b/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) diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index d381754e..94d90f13 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -163,6 +163,11 @@ reach both reports; bin/clinvar_hits.awk's sites mode words the calls the same way, and a reference call below the GQ floor reads "ref (low confidence)". + 29. step 05's checked table keeps AnnotSV's FORMAT and sample columns, so + the genotype and phase of a deletion reach the two-hit rule: a MUTYH + deletion phased 0|1 and a ClinVar hit phased 1|0 in the same phase set + are biallelic, in trans; the same deletion genotyped 0/0 is no + deletion and pairs with nothing. Run: python3 tests/test_collect_summary.py """ @@ -1177,6 +1182,40 @@ def annotsv_section(work, schema): and "Class 4 and 5 (AnnotSV)" not in txt + html_, summ["not_run"]) +# MUTYH spans chr1:45329163-45340893 (assets/carrier_genes.tsv); the calls are made up. +PHASE_HEAD = ("AnnotSV_ID\tSV_chrom\tSV_start\tSV_end\tSV_length\tSV_type\tID\tINFO\tFORMAT\tS\t" + "Annotation_mode\tGene_name\tB_loss_AFmax\tACMG_class\tAnnotSV_ranking_score\n") + + +def annotsv_phase_pairing(work): + """29. a deletion's GT and PS go through annotsv_postfilter.py into the two-hit rule.""" + import annotsv_postfilter + + def run(name, gt): + d = os.path.join(work, "sv-phase", name, "S") + put(f"{d}/annotsv/S_sv_class4_5.tsv", PHASE_HEAD + + f"1_45335001_45337000_DEL_1\t1\t45335000\t45337000\t-2000\tDEL\tSV1\t" + f"END=45337000;SVTYPE=DEL;SVCALLER=manta\tGT:PS\t{gt}:45330000\tfull\tMUTYH\t\t4\t0.95\n") + annotsv_postfilter.main(["--class45", f"{d}/annotsv/S_sv_class4_5.tsv", + "--out", f"{d}/annotsv/S_sv_class4_5_checked.tsv"]) + put(f"{d}/clinvar/S_clinvar_hits.vcf", VCF_HEAD + + "chr1\t45332000\t.\tC\tT\t50\tPASS\tGENEINFO=MUTYH:4595;CLNSIG=Pathogenic;" + "CLNREVSTAT=criteria_provided,_multiple_submitters,_no_conflicts\tGT:GQ:DP:PS\t1|0:50:30:45330000\n") + with open(f"{d}/annotsv/S_sv_class4_5_checked.tsv") as f: + head = f.readline().rstrip("\n").split("\t") + two = (collect_summary.sec_carrier(d, "S")[1].get("two_hits") or {}).get("genes", []) + return head, [x for x in two if x["gene"] == "MUTYH"] + + head, t = run("trans", "0|1") + check("sv phase: the checked table ends with AnnotSV's FORMAT and sample columns", head[-2:] == ["FORMAT", "S"], + head) + check("sv phase: a deletion phased 0|1 and a hit phased 1|0 in one phase set are biallelic, in trans", + len(t) == 1 and t[0]["verdict"] == "biallelic" and "in trans" in t[0]["why"] + and "phase set 45330000" in t[0]["why"], t) + _, t = run("homref", "0/0") + check("sv phase: a deletion genotyped 0/0 pairs with nothing", t == [], t) + + def report_followups_sections(work, schema): """26. step 42's genotype rules, step 9b's motif table of the other flagged motif loci and the loci the repeat catalog leaves without @@ -2719,6 +2758,7 @@ def crow(gene, result, risk, fa, fb): clinical_mane_quality(work, schema) comphet_phase_section(work, schema) annotsv_section(work, schema) + annotsv_phase_pairing(work) report_followups_sections(work, schema) str_extra_readings(work, schema) parascopy_section(work, schema)