diff --git a/.github/workflows/e2e.yml b/.github/workflows/e2e.yml index d7835267..760b239d 100644 --- a/.github/workflows/e2e.yml +++ b/.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 diff --git a/bin/annotsv_postfilter.py b/bin/annotsv_postfilter.py index f09905a7..004a021a 100755 --- a/bin/annotsv_postfilter.py +++ b/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,7 +98,8 @@ 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 @@ -102,6 +107,10 @@ class 4 or 5 without being a real rare deletion: 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)): diff --git a/bin/archaic_report.py b/bin/archaic_report.py index 2ad55061..6ada76bd 100755 --- a/bin/archaic_report.py +++ b/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..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__": diff --git a/bin/cnvpytor_to_vcf.py b/bin/cnvpytor_to_vcf.py index 6747c720..873ca582 100755 --- a/bin/cnvpytor_to_vcf.py +++ b/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 [] diff --git a/bin/owngenome_recheck.py b/bin/owngenome_recheck.py index 6dd4c69d..5e940be7 100755 --- a/bin/owngenome_recheck.py +++ b/bin/owngenome_recheck.py @@ -22,6 +22,8 @@ in-pipeline-release a newer pin is in a pipeline release: update the pipeline upstream-only upstream has a newer release that no pipeline release pins no-newer-release the lookup answered and nothing is newer + unknown a pin with no upstream lookup that the run's manifest + does not name: nothing was compared error the lookup failed (never read as up to date); exit 1 The steps an item feeds come from the reads column of assets/steps.tsv, then every step downstream of them through its upstream column, then the reports. @@ -37,7 +39,9 @@ check named, and lists the named script steps (04b, 14, 29, 03a-03d), which skip when their output exists, as not redone. --refresh-clinvar runs setup.sh --refresh clinvar first, and refuses while a run of any sample of -GENOME_DIR is going on, since they share the ClinVar files. Standard library only. +GENOME_DIR is going on, since they share the ClinVar files. It holds +GENOME_DIR/clinvar/.refresh.lock while it refreshes: run-all.sh and step 06 +wait on it. Standard library only. """ import argparse import fcntl @@ -66,8 +70,11 @@ REPORT_STEP = "24" # the HTML report; the text report is written after it # ClinVar's ##fileDate is the release; NCBI's Last-Modified comes a day or two later. CLINVAR_SLACK_DAYS = 2 -CLASSES = ("rerun-ready", "refresh-available", "in-pipeline-release", "upstream-only", "no-newer-release", "error") +CLASSES = ("rerun-ready", "refresh-available", "in-pipeline-release", "upstream-only", "no-newer-release", "unknown", + "error") LOG_HEADER = "utc\tevent\titem\tclass\tdetail\n" +# Held exclusive by --refresh-clinvar, shared by scripts/run-all.sh and scripts/06-clinvar-screen.sh. +CLINVAR_LOCK = "clinvar/.refresh.lock" CLINVAR_FILES = ("clinvar/clinvar_pathogenic_chr.vcf.gz", "clinvar/clinvar_chr.vcf.gz", "clinvar/clinvar.vcf.gz") # The pinned items: (manifest key, label, freshness lookup, how to read its upstream). PINNED = ( @@ -290,7 +297,10 @@ def pgs(self, pgs_id): except self.fr.LOOKUP_ERRORS as e: raise Failed(str(e)) import email.utils - return email.utils.parsedate_to_datetime(lm).strftime("%Y%m%d") + try: # a malformed date is a failed lookup of this score, never a crash of the whole check + return email.utils.parsedate_to_datetime(lm).strftime("%Y%m%d") + except (TypeError, ValueError, AttributeError, OverflowError): + raise Failed(f"cannot read the Last-Modified of {url}: {lm!r}") def release(self): try: @@ -334,7 +344,10 @@ def pgs(self, pgs_id): def release(self): e = self._get("release") - return e["tag"], dict(e.get("pins") or {}) + tag, pins = e.get("tag"), e.get("pins") or {} + if not isinstance(tag, str) or not tag or not isinstance(pins, dict): + raise Failed("the offline answers' release needs a tag and a pins object") + return tag, dict(pins) # --- the check ----------------------------------------------------------------------------- @@ -362,6 +375,10 @@ def classify(item, release_tag, release_pins, release_error, clinvar=False): ran, inst, up = item.get("ran_with"), item.get("installed"), item.get("upstream") if inst and ran and ver(inst) > ver(ran): return "rerun-ready" + # An image re-pinned at the same tag (a new digest): the tags compare equal, the references do not. + ran_ref, inst_ref = item.get("ran_with_ref"), item.get("installed_ref") + if inst and ran and ver(inst) == ver(ran) and ran_ref and inst_ref and ran_ref != inst_ref: + return "rerun-ready" # ClinVar's date is in the manifest whenever step 06 had a file, so a run without one cannot be compared # with the file here. A pin with no run value (HLA data not installed then) compares this checkout's. if clinvar and not ran: @@ -391,6 +408,21 @@ def classify(item, release_tag, release_pins, release_error, clinvar=False): return "upstream-only" +def steps_not_run(sample, sample_dir): + """The steps this run did not run: the ones run-all.sh recorded as skipped (no BAM, missing + data), and the ones it wrote no row for, which is how it leaves out an opt-in step not in TOOLS + (run_state's "not planned"). Empty when logs/run_status.tsv records no step at all (a run of an + older version), so nothing is left out on a guess.""" + try: + status = run_state.read_run_status(os.path.join(sample_dir, "logs", "run_status.tsv")) + if not (status or {}).get("steps"): + return set() + return {s["step"] for s in run_state.collect(sample, sample_dir)["steps"] + if s["state"] in ("skipped", "not planned")} + except (OSError, ValueError, KeyError): + return set() + + def action(item, sample): c, key = item["class"], item["item"] if c == "rerun-ready": @@ -415,6 +447,9 @@ def action(item, sample): return "no pipeline release pins it yet; nothing to do until one does" if c == "error": return "the lookup failed; check again later" + if c == "unknown": + return ("the run's manifest does not name this pin (a run of an older version): nothing to compare; " + "the next run records it") return "" @@ -446,18 +481,31 @@ def add(item, fn, *args, clinvar=False): ran = pin_value(key, man.get(("data", key)) if key == "hla_database" else version.get(key)) inst = this_hla if key == "hla_database" else pin_value(key, this_pins.get(key)) item = {"item": key, "label": label, "kind": "pin", "ran_with": ran, "installed": inst} + if key.endswith("_IMAGE"): + item["ran_with_ref"], item["installed_ref"] = version.get(key) or None, this_pins.get(key) or None if key == "hla_database" and not this_hla: item["error"] = "could not read this checkout's HLA_DB_RELEASE (scripts/lib/common.sh)" add(item, answers.pinned, key, lookup, pattern) # The other versions.env pins a step reads (steps 50 to 53): no freshness source exists for them, so - # nothing is looked up; a pin here that differs from the run's makes its step rerun-ready. + # nothing is looked up; a pin here that differs from the run's (an image by its full reference, so a + # new digest at the same tag counts) makes its step rerun-ready. A pin the run's manifest lacks is + # unknown, never no-newer-release. A pin read only by steps this run did not run (an opt-in step not + # in TOOLS, or one skipped for missing data) is left out: a rerun would not run that step either. known = {k for k, *_ in PINNED} + skipped = steps_not_run(sample, sample_dir) for key in dict.fromkeys(r[:-5] if r.endswith("@self") else r for s in steps for r in reads_of(s)): if key in known or key not in this_pins: continue - ran, inst = pin_value(key, version.get(key)), pin_value(key, this_pins.get(key)) + readers = [s["step"] for s in steps if any(key_matches(key, r) for r in reads_of(s))] + if readers and all(n in skipped for n in readers): + continue + ran_ref, inst_ref = version.get(key) or None, this_pins.get(key) or None + ran, inst = pin_value(key, ran_ref), pin_value(key, inst_ref) + cls = "unknown" if not (ran_ref and inst_ref) else "rerun-ready" if ran_ref != inst_ref else "no-newer-release" item = {"item": key, "label": f"{key} (no upstream lookup)", "kind": "pin", "ran_with": ran, - "installed": inst, "class": "rerun-ready" if ran and inst and ran != inst else "no-newer-release"} + "installed": inst, "class": cls} + if key.endswith("_IMAGE"): + item["ran_with_ref"], item["installed_ref"] = ran_ref, inst_ref item["steps"] = affected_steps(steps, [key]) item["action"] = action(item, sample) items.append(item) @@ -520,6 +568,9 @@ def print_check(state, out=None): w(f"Database re-check of {state['sample']} ({state['checked_utc']}):\n") for i in state["items"]: detail = i.get("error") or f"run used {i.get('ran_with') or 'unknown'}, here {i.get('installed') or 'unknown'}, upstream {i.get('upstream') or 'unknown'}" + if not i.get("error") and i.get("ran_with_ref") and i.get("installed_ref") and i["ran_with_ref"] != i["installed_ref"] \ + and i.get("ran_with") == i.get("installed"): + detail += f"; same tag, new image: run used {i['ran_with_ref']}, here {i['installed_ref']}" w(f" {i['label']}: {i['class']} ({detail})\n") if i["class"] not in ("no-newer-release",) and i.get("action"): w(f" {i['action']}\n") @@ -624,6 +675,31 @@ def running_samples(genome_dir): and run_state.collect(n, os.path.join(genome_dir, n))["run"]["state"] == "running"] +def refresh_clinvar(genome_dir): + """setup.sh --refresh clinvar under the ClinVar lock of GENOME_DIR (CLINVAR_LOCK). The ClinVar files + are shared by every sample, so a run of any of them blocks the refresh: run-all.sh and step 06 hold + the lock shared while they read the files, and wait while this holds it. A run with no lock (a bare + nextflow run, an older checkout) is caught by its logs/run_status.tsv. The lock is let go before + --rerun starts run-all.sh, which takes it shared.""" + os.makedirs(os.path.join(genome_dir, "clinvar"), exist_ok=True) + with open(os.path.join(genome_dir, CLINVAR_LOCK), "a") as lock: + try: + fcntl.flock(lock, fcntl.LOCK_EX | fcntl.LOCK_NB) + except OSError: + print(f"ERROR: a run or step 06 holds {os.path.join(genome_dir, CLINVAR_LOCK)}: not refreshing the " + "ClinVar files it reads.", file=sys.stderr) + return 3 + busy = running_samples(genome_dir) + if busy: + print(f"ERROR: a run of {', '.join(busy)} is going on: not refreshing the ClinVar files it reads.", + file=sys.stderr) + return 3 + r = subprocess.run([os.path.join(ROOT, "scripts", "setup.sh"), "--refresh", "clinvar", genome_dir]).returncode + if r != 0: + print("ERROR: setup.sh --refresh clinvar failed: see above.", file=sys.stderr) + return r + + def main(argv=None): ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) ap.add_argument("--sample") @@ -659,16 +735,8 @@ def main(argv=None): print(f"ERROR: no {sample_dir}/run_manifest.tsv: run the pipeline on {a.sample} first.", file=sys.stderr) return 2 if a.refresh_clinvar: - # The ClinVar files are shared by every sample of GENOME_DIR, so a run of any of them blocks the - # refresh. This narrows the window; a run that starts during the refresh is not stopped. - busy = running_samples(a.genome_dir) - if busy: - print(f"ERROR: a run of {', '.join(busy)} is going on: not refreshing the ClinVar files it reads.", - file=sys.stderr) - return 3 - r = subprocess.run([os.path.join(ROOT, "scripts", "setup.sh"), "--refresh", "clinvar", a.genome_dir]).returncode + r = refresh_clinvar(a.genome_dir) if r != 0: - print("ERROR: setup.sh --refresh clinvar failed: see above.", file=sys.stderr) return r fr = load_freshness() this_pins, _ = fr.read_pins(os.path.join(ROOT, "versions.env")) diff --git a/bin/sv_svlen_fix.py b/bin/sv_svlen_fix.py index 30d8ac56..7c119463 100755 --- a/bin/sv_svlen_fix.py +++ b/bin/sv_svlen_fix.py @@ -18,6 +18,12 @@ breakend records pass unchanged. A rewritten record keeps its old values in ORIG_SVLEN and ORIG_END. +The input holds one sample. Step 05 and the reports read the genotype and +phase of the first sample column (bin/annotsv_postfilter.py sample_format), +so a VCF with more sample columns would let another person's genotype, or a +0/0, decide whether a deletion counts: such an input stops the script (exit +2) with the sample names. + Usage: sv_svlen_fix.py IN.vcf[.gz]|- > OUT.vcf """ import gzip @@ -82,6 +88,12 @@ def main(argv): for line in fh: if line.startswith("#"): if line.startswith("#CHROM"): + samples = line.rstrip("\n").split("\t")[9:] + if len(samples) > 1: + sys.stderr.write("ERROR: the SV VCF holds {0} samples ({1}); AnnotSV's input must hold one. " + "Pick it first: bcftools view -s \n".format( + len(samples), ", ".join(samples[:5]) + (", ..." if len(samples) > 5 else ""))) + return 2 for h in NEW_INFO: if h not in seen: out.write(h + "\n") diff --git a/docs/05-annotsv.md b/docs/05-annotsv.md index 2bd0da75..7c343a39 100644 --- a/docs/05-annotsv.md +++ b/docs/05-annotsv.md @@ -1,7 +1,7 @@ # Step 5: Structural Variant Annotation (AnnotSV) ## What This Does -Classifies every structural variant from Manta that passes duphold's depth filter (step 15; when that file is missing or older than Manta's calls, the script annotates Manta's calls and says so), plus the Delly and CNVpytor deletions and duplications that no Manta call covers, using ACMG guidelines (class 1-5), adding gene overlap, population frequency, and clinical significance. +Classifies every structural variant from Manta that passes duphold's depth filter (step 15; when that file is missing or older than Manta's calls, the script annotates Manta's calls and says so), plus the Delly and CNVpytor deletions and duplications that no Manta call covers (a Manta record that is filtered or genotyped 0/0 covers nothing), using ACMG guidelines (class 1-5), adding gene overlap, population frequency, and clinical significance. ## Why Raw Manta output contains thousands of SVs with no clinical interpretation. AnnotSV tells you which ones matter by cross-referencing known pathogenic SVs, gene databases, and population data. @@ -45,7 +45,7 @@ Step 05 and the Nextflow `ANNOTSV` stop with a pointer to `--annotsv-data` when AnnotSV classifies only the records of the one VCF it is given. The step builds that VCF, `annotsv/_annotsv_input.vcf`, in two passes inside the AnnotSV image: - **The net length change** (`bin/sv_svlen_fix.py`). Manta writes a deletion that also inserts a few bases with its REF and ALT sequences, `END` at the last deleted base and `SVLEN` as the deleted length. AnnotSV measures the overlapped coding sequence over `POS-1..END`, so REF 137 bp with ALT 2 bp read as 137 deleted coding bases: `Frameshift` "yes" and criterion 2E-1 (+0.90, class 4), although the change is -135 bp and keeps the frame. For every record whose REF and ALT are both sequence, the script sets `SVLEN` to `len(ALT) - len(REF)` and, for a net deletion, `END` to `POS - 1 + len(REF) - len(ALT)`: one base less than AnnotSV's own default `END` (`POS + len(REF) - len(ALT)`), so that its `POS-1..END` interval, which counts the anchor base, is as long as the net change. The old values stay in `ORIG_SVLEN` and `ORIG_END`. Symbolic records (``) and breakends pass unchanged. -- **The other callers' calls** (`bin/cnvpytor_to_vcf.py`). A large deletion that CNVpytor (step 18) or Delly (step 19) calls, but Manta calls as something else or misses, got no class before. The script adds Delly's PASS deletions and duplications that read depth does not contradict (below), then CNVpytor's with e-val (`pytorP1`) below 0.01 and a fraction of mapping-quality-0 reads (`pytorQ0`) below 0.5, each only when no record already in the file has the same type on the same chromosome with a reciprocal overlap of at least half the longer call. A deletion both of them call is added once, as Delly's. CNVpytor's copy-neutral LOH calls are left out. Each file is used when it exists, and the step says when one is missing. +- **The other callers' calls** (`bin/cnvpytor_to_vcf.py`). A large deletion that CNVpytor (step 18) or Delly (step 19) calls, but Manta calls as something else or misses, got no class before. The script adds Delly's PASS deletions and duplications with a non-reference genotype (not 0/0 or ./.) that read depth does not contradict (below), then CNVpytor's with e-val (`pytorP1`) below 0.01 and a fraction of mapping-quality-0 reads (`pytorQ0`) below 0.5, each only when no record already in the file has the same type on the same chromosome with a reciprocal overlap of at least half the longer call. A deletion both of them call is added once, as Delly's. CNVpytor's copy-neutral LOH calls are left out. Each file is used when it exists, and the step says when one is missing. - **Read depth for long Delly calls.** Short-read Delly can mark PASS a long deletion or duplication, even one of many megabases, whose read depth shows no change, and AnnotSV then classes it as the loss or gain of every gene it spans. So a Delly call longer than 1,000 bp joins the input only when Delly's own read-depth copy number (`FORMAT/RDCN`) agrees: below a normal region's copies for a deletion, above them for a duplication. The normal is 1 on a male's chrX or chrY outside the pseudoautosomal regions (`assets/par_grch38.bed`) and 2 everywhere else. With no sex given, a stretch of chrX or chrY outside them counts 1 and 2 as normal, so there a deletion needs RDCN 0 or 1 and a duplication 2 or more: give the sex (the second argument, or the samplesheet's `sex` column) for the strict check. A long call with no estimate (`.` or a negative RDCN, no reads in the flanks) is set aside too. A call of 1,000 bp or less rests on its split reads and read pairs alone, and a Delly VCF without the RDCN field is read as before. The couple screen of [step 39](39-couple.md) uses the same rule, from the same code. The step's log line `Callers beside manta: ...`, repeated in the input VCF's `##svcallers_added` header line, counts the Delly calls set aside. Every record names its caller in `INFO/SVCALLER` (`manta`, `delly` or `cnvpytor`), which AnnotSV copies into its `INFO` column. With `SV_VCF` set, that file alone is annotated, after the length pass, with `SVCALLER=sv_vcf`. @@ -91,7 +91,7 @@ docker run --rm --user "$(id -u):$(id -g)" \ -bracketedAltMode DEL,DUP,INS ``` -To annotate another SV VCF (for example Sniffles2 output), set `SV_VCF` to its host path; the file must be inside `${GENOME_DIR}`. +To annotate another SV VCF (for example Sniffles2 output), set `SV_VCF` to its host path; the file must be inside `${GENOME_DIR}` and hold one sample. The reports read the genotype and phase of the first sample column, so a VCF with more sample columns stops the step (`bin/sv_svlen_fix.py`) instead of letting another person's genotype decide; pick the sample first with `bcftools view -s `. ## Output - `${GENOME_DIR}/${SAMPLE}/annotsv/${SAMPLE}_sv_annotated.tsv` — main annotated output (one row per SV, with ACMG class) @@ -118,7 +118,7 @@ AnnotSV gives class 4 or 5 to some calls that are not a rare deletion: a call on Some class 4 and 5 calls are common or unreliable for a reason AnnotSV does not weigh. The table keeps them, with `ACMG_class` 3, AnnotSV's class in `AnnotSV_class` and the reason in `Rule_detail`. The reports list such a call apart from the class 4 and 5 section, with AnnotSV's class, the rule and `Rule_detail`, and leave it out of the section's count. -- **No read-depth support** (`no_depth_support`). A real deletion lowers the read depth and a real duplication raises it. `Depth_support` reads the record's own value: duphold's `DHFFC` below 0.7 for a deletion and `DHBFC` above 1.3 for a duplication, the thresholds of [step 15](15-duphold.md) (Pedersen and Quinlan 2019, PMID 31222198); for a Delly call over 1,000 bp without them, Delly's `RDCN`, 0 or 1 copies for a deletion and 3 or more for a duplication (2 or more on chrX and chrY, since the table does not know the sex). A CNVpytor call is a read-depth call, so it has support. A Delly call step 05 added to AnnotSV's input keeps only `FORMAT/GT`, so its `RDCN` does not reach the table; over 1,000 bp it says `yes`, because `bin/cnvpytor_to_vcf.py` already set aside every such call whose `RDCN` does not show the change. `RDCN` itself is read here only for `SV_VCF` input, where the caller's FORMAT fields survive. `no` moves the call to class 3. A record with no depth value says `not checked` and keeps its class: that happens for a Delly call of 1,000 bp or less, when step 15 did not run (the step then annotates Manta's calls and says so), or for `SV_VCF` input without these fields. Step 15 already drops the Manta calls its depth contradicts, so in a default run this check mostly confirms; it catches an `SV_VCF` file duphold annotated but did not filter. +- **No read-depth support** (`no_depth_support`). A real deletion lowers the read depth and a real duplication raises it. `Depth_support` reads the record's own value: duphold's `DHFFC` below 0.7 for a deletion and `DHBFC` above 1.3 for a duplication, the thresholds of [step 15](15-duphold.md) (Pedersen and Quinlan 2019, PMID 31222198); for a Delly call over 1,000 bp without them, Delly's `RDCN`, judged by the rule `bin/cnvpytor_to_vcf.py` applies when it adds a Delly call: fewer copies than normal for a deletion and more for a duplication, where normal is 1 on a male's chrX and chrY outside the pseudoautosomal regions and 2 everywhere else, PARs included. Step 05 and `ANNOTSV` pass the sample's sex and `assets/par_grch38.bed` for this. With no sex given, 1 and 2 both count as normal on chrX and chrY outside the PARs, so there a deletion needs 0 or 1 copies and a duplication 2 or more. A CNVpytor call is a read-depth call, so it has support. A Delly call step 05 added to AnnotSV's input keeps only `FORMAT/GT`, so its `RDCN` does not reach the table; over 1,000 bp it says `yes`, because `bin/cnvpytor_to_vcf.py` already set aside every such call whose `RDCN` does not show the change. `RDCN` itself is read here only for `SV_VCF` input, where the caller's FORMAT fields survive. `no` moves the call to class 3. A record with no depth value says `not checked` and keeps its class: that happens for a Delly call of 1,000 bp or less, when step 15 did not run (the step then annotates Manta's calls and says so), or for `SV_VCF` input without these fields. Step 15 already drops the Manta calls its depth contradicts, so in a default run this check mostly confirms; it catches an `SV_VCF` file duphold annotated but did not filter. - **Inside a coding tandem repeat** (`tandem_repeat`). An SV whose two breakpoints lie within 10 bp of a tandem repeat that overlaps a MANE Select coding exon changes the repeat's length. Short-read SV callers place such breakpoints loosely and do not size repeats, and AnnotSV can read the change as a frameshift. A length change of a repeat is a question for a repeat genotyper (the repeats of [step 9](09-str-expansions.md)'s catalog are read there), not an SV call. ACMG's PM4 and BP3 (Richards 2015, PMID 25741868) already weigh in-frame changes inside a repeat differently from other length changes. - **The IG and TR loci** (`ig_tr_locus`). An SV more than half inside IGH, IGK, IGL, TRA (with TRD inside it), TRB or TRG. B and T cells cut DNA out of these loci when they rearrange their receptor genes (V(D)J recombination), so DNA from blood carries such deletions in a share of its cells without any inherited change. - **The MHC** (`mhc`). An SV more than half inside the MHC, GRCh38 chr6:28,510,120-33,480,577. Its haplotypes differ by tens of kilobases between people, and gnomAD-SV marks every SV over half inside IGH or the MHC as low confidence (FILTER `IGH_MHC_OVERLAP`, Collins 2020, PMID 32461652). CYP21A2 lies in the MHC; short reads cannot call it here, and [step 35b](35-paralogs.md) reads it from PacBio HiFi reads. diff --git a/docs/43-archaic.md b/docs/43-archaic.md index c4f45847..7dbd38a6 100644 --- a/docs/43-archaic.md +++ b/docs/43-archaic.md @@ -57,13 +57,13 @@ Chen et al. do not score every site. Each archaic genome comes with its own mask - a genome's call outside that genome's FilterBed is kept and marked masked in the table (10, 11 or 12: 10 plus its ALT count); - every call at a site inside a segmental duplication (UCSC's GRCh37 track), or inside the window of a 1000 Genomes indel, is marked masked the same way. The window is 5 bp on each side of the indel, reaching past the longer of its REF and ALT, as the published script writes it; - step 43 writes each masked call to IBDmix as its ALT count, and the masked positions to a mask file per chromosome and archaic genome, which it passes to IBDmix with `-r`; -- a table row outside the 1000 Genomes GRCh38 strict mask is dropped. +- every call at a site outside the 1000 Genomes GRCh38 strict mask is marked masked the same way, in all three archaic columns. `setup.sh --archaic` does that after the liftover (`archaic_report.py strict`), because the mask is on GRCh38 and the table is built on GRCh37. -The first version of the masks dropped the masked calls instead (code 9, or no row). IBDmix then had nothing to score at those sites, the mismatches included, so a segment could run through masked sequence that contradicts it. That is one reason HG002's share stayed above the published figure after the masks; see Runtime below. The strict mask still drops rows rather than marking them: that filter runs in `setup.sh` after the liftover, and moving it to the same marking is the remaining difference from the published workflow in how masked sites are handled. +The first version of the masks dropped the masked calls instead (code 9, or no row), and a second one still dropped every row outside the strict mask. IBDmix then had nothing to score at those sites, the mismatches included, so a segment could run through masked sequence that contradicts it. That is why HG002's share stayed above the published figure after the masks; see Runtime below. One of their exclusions is not applied here: the CpG mask, which needs a GRCh37 reference and the chimpanzee, orangutan and macaque SNP sets, none of which this pipeline installs. What that leaves out is written down in the Notes below. -The filters cut hard: on chr22 the strict mask covers 25.2 Mb of 50.8, the Vindija FilterBed 19.6 Mb, and segmental duplications take out 4.3 Mb. Without them the step called sequence from positions the archaic genomes cannot be read at, and HG002's Neanderthal total came out about twice the published figure (see Runtime below). The filters are part of the data version: a change to a pin or to a mask means `snpad2016-1kgp3-v3` and a rebuild. `setup.sh --archaic` also rebuilds a table whose build code changed (its stamp holds the sha256 of the table half of `bin/archaic_report.py`), so a table built before the masked codes is rebuilt on the next `setup.sh --archaic`. Until then step 43 reads it as before: a table with no masked code gives empty masks. +The filters cut hard: on chr22 the strict mask covers 25.2 Mb of 50.8, the Vindija FilterBed 19.6 Mb, and segmental duplications take out 4.3 Mb. Without them the step called sequence from positions the archaic genomes cannot be read at, and HG002's Neanderthal total came out about twice the published figure (see Runtime below). The filters are part of the data version: a change to a pin or to a mask means `snpad2016-1kgp3-v3` and a rebuild. `setup.sh --archaic` also rebuilds a table whose build code changed (its stamp holds the sha256 of the table half of `bin/archaic_report.py`), so a table built before the masked codes, or before the strict-mask rows were kept, is rebuilt on the next `setup.sh --archaic`. Until then step 43 reads it as before: a table with no masked code gives empty masks. The table is built on GRCh37 and lifted once to GRCh38 with Picard `LiftoverVcf`. A site that does not lift, or lifts to another chromosome, is dropped. The archaic files are downloaded one chromosome at a time (150 GB in all, at most 13.8 GB at once, for chromosome 2) and deleted once the chromosome is built, so an interrupted install resumes at the chromosome it stopped in. `ARCHAIC_REGIONS` (a GRCh37 BED, chromosomes named 1 to 22) builds only those regions. @@ -101,7 +101,7 @@ Measured on public HG002 by the `archaic-hg002` jobs of `public-evidence.yml`, o That run was made before the masks, from a table that held every lifted site. It gave HG002 115.2 Mb of Neanderthal sequence over the 22 autosomes, about twice what Chen et al. report for Europeans with the same method and cutoffs (about 51 Mb per person, and 49.9 to 82.9 Mb across their released per-person calls). The cause was the missing masks: the table carried the positions each archaic genome's own filters exclude, and IBDmix scored them. On chr22, 43% of that run's Neanderthal sequence lay outside the 1000 Genomes GRCh38 strict mask alone. -With the masks dropping sites, a run on chr22 (run 42089578470, before the masked codes) gave HG002 1.3 Mb of Neanderthal sequence with the `ALL` frequencies and 1.7 Mb with `EUR`, still at the top of the published spread below. Keeping the masked sites as IBDmix's mask, as the published workflow does, lets the mismatches there end a segment; the chr22 figure on this table is printed by the e2e case on every run. The e2e case holds it inside the published spread: HG002's chr22 Neanderthal total must land inside the spread of Chen et al.'s own released calls for the 497 Europeans of CEU, TSI, GBR, FIN and IBS (51 kb at the lowest, 541 kb at the median, 1,132 kb at the 95th percentile, 1,716 kb at the highest). The pre-mask run gave 1,748 kb on chr22, more than all 497 of them. On this table the e2e fork run gave 1,177 kb: closer to the published figures, but still above the 95th percentile and about twice the median, so the gap is narrowed, not closed. The strict mask still drops rows in `scripts/setup.sh` (marking them instead is the remaining step), and the whole-genome total must land in 49.9 to 82.9 Mb before the share counts as in line with the published one. The case prints where HG002 falls inside that spread; it does not gate on the percentile, for two reasons: one person's chromosome is one draw from that distribution, and the published calls come from imputed low-coverage genotypes, which lose rare alleles and so lose archaic ones, while HG002 is a high-coverage benchmark. The whole-genome figure is measured by a `public-evidence.yml` run, which installs the 150 GB of archaic data once per release. +With the masks dropping sites, a run on chr22 (run 42089578470, before the masked codes) gave HG002 1.3 Mb of Neanderthal sequence with the `ALL` frequencies and 1.7 Mb with `EUR`, still at the top of the published spread below. Keeping the masked sites as IBDmix's mask, as the published workflow does, lets the mismatches there end a segment; the chr22 figure on this table is printed by the e2e case on every run. The e2e case compares it with the spread of Chen et al.'s own released calls for the 497 Europeans of CEU, TSI, GBR, FIN and IBS (51 kb at the lowest, 541 kb at the median, 1,132 kb at the 95th percentile, 1,716 kb at the highest). The pre-mask run gave 1,748 kb on chr22, more than all 497 of them. With the masked codes but the strict-mask rows still dropped, the fork runs gave 1,177 to 1,190 kb. Keeping those rows marked, the last difference from the published workflow in how masked sites are handled, gave 1,138 kb on chr22 (run 42335592918). One chromosome of one person varies widely, so the case gates on a regression guard of 1,200 kb (and on the lowest published value, 51 kb), and prints the published figures beside the measured one. It does not claim the figure is inside the published range. The whole-genome total (49.9 to 82.9 Mb across the published per-person calls) is checked by the release run of `public-evidence.yml`, which installs the 150 GB of archaic data once per release, not by the e2e case; that run decides whether the share is in range. The case prints where HG002 falls inside the spread; it does not gate on the median or the 95th percentile, for two reasons: one person's chromosome is one draw from that distribution, and the published calls come from imputed low-coverage genotypes, which lose rare alleles and so lose archaic ones, while HG002 is a high-coverage benchmark. ## Notes @@ -110,9 +110,9 @@ With the masks dropping sites, a run on chr22 (run 42089578470, before the maske - A segment counts when either of your two copies carries it. The share is of the analysed span, not of the whole genome. - The table has autosomes only: chrX and chrY are not assessed. - Sites the VCF has no record for are read as the reference allele. A VCF that holds only variants (no reference calls) cannot tell a reference call from a gap in coverage. -- The share is of the analysed span, and the span still runs from the first to the last table site of each chromosome whatever the masks remove. A segment runs from its first site to the site after its last, so it can cross masked sequence and count those bases (a homozygous mismatch at a masked site still ends it, as in the published workflow); a wrong genotype at a site that was kept can raise the share as well. The number is an estimate in both directions, not a floor. Compare it with the published per-person figures, which come from masked data too; today it is closer to them than before, not yet in line. +- The share is of the analysed span, and the span still runs from the first to the last table site of each chromosome whatever the masks remove. A segment runs from its first site to the site after its last, so it can cross masked sequence and count those bases (a homozygous mismatch at a masked site still ends it, as in the published workflow); a wrong genotype at a site that was kept can raise the share as well. The number is an estimate in both directions, not a floor. Compare it with the published per-person figures, which come from masked data too. - The CpG sites Chen et al. exclude are kept here. A C to T change at a CpG site is the commonest ancient-DNA error, so a few of the archaic genotypes at those sites are wrong in both directions, and the sites are over-represented among the matches a method like this finds. Their own filter needs data this pipeline does not install, so the share carries that error; it is one of the reasons a number from this step is an estimate. -- On the e2e fixture the case builds the table for GRCh37 chr22 and checks HG002's GIAB v4.2.1 genotypes against the published chr22 spread, a copy of HG002 with one chromosome over chr22:28-29 Mb taken from Vindija (a Neanderthal segment must cover it), a person built from the table who carries only the majority allele (a quarter of HG002's Neanderthal sequence at most), that no table row falls outside the strict mask, and that the FilterBeds and the exclusions marked calls masked. +- On the e2e fixture the case builds the table for GRCh37 chr22 and checks HG002's GIAB v4.2.1 genotypes against the published chr22 spread, a copy of HG002 with one chromosome over chr22:28-29 Mb taken from Vindija (a Neanderthal segment must cover it), a person built from the table who carries only the majority allele (a quarter of HG002's Neanderthal sequence at most), that every table row outside the strict mask carries only masked calls, and that the FilterBeds and the exclusions marked calls masked. ## References diff --git a/docs/recheck.md b/docs/recheck.md index e3ece402..221c6fc0 100644 --- a/docs/recheck.md +++ b/docs/recheck.md @@ -23,6 +23,8 @@ Standard library only, so the plain `python3` on the host is enough. It needs th | IPD-IMGT/HLA | `hla_database` | The `release_version.txt` of the IPD-IMGT/HLA repository | 08 | | Each PGS score | `pgs:`, its `HmPOS_date` | The PGS Catalog's REST API, then the `Last-Modified` date of the harmonised GRCh38 file | 25, then 26 | +The pins of steps 50 to 53 (`MEGANE_IMAGE`, `MOCHA_PANEL_VERSION`, `MOCHA_IMAGE`, `VNTYPER_IMAGE`, `LOCITYPER_PANEL_VERSION`, `LOCITYPER_IMAGE`) have no upstream lookup: each is compared only with this checkout's `versions.env`. An image is compared by its full reference, so a new `@sha256` digest at the same tag counts as a change, for these and for `PHARMCAT_IMAGE`. A pin read only by steps this run did not run (an opt-in step not in `TOOLS`, or one skipped for missing data) is not listed: a rerun would not run that step either. + The HTML report (24) and the text report are written again after any of them. The ClinVar release inside the PCGR bundle (`cpsr_clinvar_release`) moves with the PCGR bundle, never with step 06's ClinVar file. The lookups are the ones the monthly CI check uses (`scripts/ci/freshness.py`). It also reads the newest release of this pipeline (`--owngenome-repo`, default `GeiserX/OwnGenome`) and that release's `versions.env`, to say whether a newer pin is already in a release. @@ -37,11 +39,12 @@ Each item gets one class: | Class | Meaning | What to do | |---|---|---| -| `rerun-ready` | What is installed here is newer than what the run used: a newer ClinVar file or PGS file on disk, or a pin in this checkout | Rerun (for a pin, run `./scripts/setup.sh` first so its data is there) | -| `refresh-available` | ClinVar only: NCBI's file is more than 2 days newer than the one installed | `./scripts/setup.sh --refresh clinvar $GENOME_DIR`, then rerun; `--rerun --refresh-clinvar` does both, and refuses while any sample of `$GENOME_DIR` is running, since they share the ClinVar files | +| `rerun-ready` | What is installed here is newer than what the run used: a newer ClinVar file or PGS file on disk, or a pin in this checkout (an image re-pinned to a new digest at the same tag too) | Rerun (for a pin, run `./scripts/setup.sh` first so its data is there) | +| `refresh-available` | ClinVar only: NCBI's file is more than 2 days newer than the one installed | `./scripts/setup.sh --refresh clinvar $GENOME_DIR`, then rerun; `--rerun --refresh-clinvar` does both, and refuses while any sample of `$GENOME_DIR` is running, since they share the ClinVar files ([the ClinVar lock](#the-clinvar-lock)) | | `in-pipeline-release` | A pipeline release pins a newer version | Update the pipeline to that release, run `setup.sh`, then rerun | | `upstream-only` | Upstream has a newer release that no pipeline release pins | Nothing yet for a pinned item; for a PGS score, download it again ([step 25](25-prs.md)) | | `no-newer-release` | The lookup answered and nothing is newer | Nothing | +| `unknown` | A pin of steps 50 to 53 that the run's manifest does not name (a run of an older version): nothing was compared | Nothing; the next run records it. It does not make the command exit 1 | | `error` | The lookup failed | Check again later. A failed lookup is never read as "up to date", and the command exits 1 | The 2-day margin for ClinVar is there because NCBI dates the file a day or two after the release in its `##fileDate`. A newer PGS file date means the catalog wrote the harmonised file again; it does not prove the weights changed. @@ -67,6 +70,10 @@ When the run ends, it compares the steps that ran a task again with the steps th A sample last run by an older `run-all.sh` has no `meta command` row: run the same command once by hand. +### The ClinVar lock + +The ClinVar files are shared by every sample of `$GENOME_DIR`. `--refresh-clinvar` holds `$GENOME_DIR/clinvar/.refresh.lock` exclusive while `setup.sh --refresh clinvar` replaces them, and lets go before the rerun starts. `run-all.sh` holds the same lock shared from its start to its end, and `scripts/06-clinvar-screen.sh` while it runs: a run that starts during a refresh waits for it ("Waiting for the ClinVar refresh ..."), and a refresh refuses (exit 3) while a run holds the lock. It also refuses while any sample's `logs/run_status.tsv` says a run is going on, which covers a run of an older checkout. The shell side needs `flock` (part of util-linux on Linux; on macOS `brew install flock`); without it nothing waits. A bare `nextflow run`, and `setup.sh --refresh clinvar` run by hand, take no lock. + ### A bare Nextflow run A bare `nextflow run main.nf` writes no run manifest, so there is nothing to check. To get the command back, `--launch-dir` prints the last command of that launch directory (its `.nextflow/history`) with `-resume` added. It never runs it. @@ -77,9 +84,10 @@ python3 bin/owngenome_recheck.py --launch-dir /path/to/launch/dir ## The files it writes -- `/recheck/state.json`, replaced whole on each check: `checked_utc`, the hosts asked, the pipeline release it read, one entry per item (`item`, `label`, `ran_with`, `installed`, `upstream`, `class`, `steps`, `action`, and `error` or `in_owngenome_release` when they apply), `affected_steps` (the steps a rerun would redo for the `rerun-ready` and `refresh-available` items), and after a rerun a `rerun` block (`started_utc`, `finished_utc`, `rc`, `command`, `predicted_vs_redone`). The last rerun stays until another one replaces it. +- `/recheck/state.json`, replaced whole on each check: `checked_utc`, the hosts asked, the pipeline release it read, one entry per item (`item`, `label`, `ran_with`, `installed`, `upstream`, `class`, `steps`, `action`, and `error` or `in_owngenome_release` when they apply; an image item also has `ran_with_ref` and `installed_ref`, the full references with any digest), `affected_steps` (the steps a rerun would redo for the `rerun-ready` and `refresh-available` items), and after a rerun a `rerun` block (`started_utc`, `finished_utc`, `rc`, `command`, `predicted_vs_redone`). The last rerun stays until another one replaces it. - `/recheck/log.tsv`, only ever appended: one row per item and check, and a row as each rerun starts and ends. - `/recheck/rerun.lock`, an empty file each `--rerun` locks while it runs. +- `$GENOME_DIR/clinvar/.refresh.lock`, an empty file: [the ClinVar lock](#the-clinvar-lock). `run-all.sh` and step 06 make it too when it is missing. ## In the reports and on the live site diff --git a/docs/testing.md b/docs/testing.md index fdabfd79..e93a6159 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -92,17 +92,17 @@ Every pull request starts the workflow, so main can require the e2e check. Its ` A pull request's e2e run used to run every case: 2 h 17 min from start to end and 339 job-minutes with free runners (run 41774860458, 2026-11-04), while most cases cannot be affected by a given change. [`scripts/ci/e2e-scope.sh`](https://github.com/GeiserX/OwnGenome/blob/main/scripts/ci/e2e-scope.sh) now picks the cases from the files the pull request changed since it left main (`git diff BASE...HEAD`, three dots): -- **None** when no changed path is in the e2e scope: `scripts/`, `modules/`, `workflows/`, `bin/`, `conf/`, `assets/`, `tests/e2e/`, `tests/fixtures/VERSION`, `main.nf`, `nextflow.config`, `nextflow_schema.json`, `versions.env`, `.github/workflows/e2e.yml`, and the committed files under `tests/fixtures/` that a case reads (each named in `SCOPE_FILES` and in that case's row; the other files there feed only the unit tests). -- **Every case** when a changed path is shared wiring (`FULL` in the script): `main.nf`, `nextflow.config`, `nextflow_schema.json`, `versions.env`, `assets/steps.tsv`, `workflows/`, `conf/`, `scripts/lib/`, `scripts/run-all.sh`, steps 02 and 03 (the BAM and VCF every case starts from), the fixture (`tests/fixtures/VERSION`, `scripts/ci/build-fixture.sh`), and the e2e runner and planner (`e2e-run.sh`, `e2e-shards.sh`, `e2e-scope.sh`, `tests/e2e/lib.sh`, `tests/e2e/bin/`, `e2e.yml`). Also when the pull request carries the label `e2e-full`: the scope job reads the labels when it runs, so add the label and re-run the workflow (or push). +- **None** when no changed path is in the e2e scope: `scripts/`, `modules/`, `workflows/`, `bin/`, `conf/`, `assets/`, `tests/e2e/`, `tests/fixtures/VERSION`, `main.nf`, `nextflow.config`, `nextflow_schema.json`, `versions.env`, `.github/workflows/e2e.yml`, `.github/actions/`, and the committed files under `tests/fixtures/` that a case reads (each named in `SCOPE_FILES` and in that case's row; the other files there feed only the unit tests). +- **Every case** when a changed path is shared wiring (`FULL` in the script): `main.nf`, `nextflow.config`, `nextflow_schema.json`, `versions.env`, `assets/steps.tsv`, `workflows/`, `conf/`, `scripts/lib/`, `scripts/run-all.sh`, steps 02 and 03 (the BAM and VCF every case starts from), the fixture (`tests/fixtures/VERSION`, `scripts/ci/build-fixture.sh`), and the e2e runner and planner (`e2e-run.sh`, `e2e-shards.sh`, `e2e-scope.sh`, `tests/e2e/lib.sh`, `tests/e2e/bin/`, `e2e.yml`) with the composite actions every leg runs (`.github/actions/`). Also when the pull request carries the label `e2e-full`: the scope job reads the labels when it runs, so add the label and re-run the workflow (or push). - **Otherwise** the smoke cases `nextflow-from-fastq-1-bash` and `nextflow-from-fastq-2-nextflow` (the bash steps and the pipeline from FASTQ, and their parity check), every case whose own file changed, every case whose row in `e2e-case-paths.tsv` or `e2e-case-times.tsv` changed, every case with a path glob that matches a changed file, and the cases those need (`needs`, followed through). The base leg runs every base case whenever any case runs, because every shard starts from what they wrote. A change to a base step other than 02 and 03 therefore runs the base cases, the smoke cases and the cases that list the step, but not a shard case that only reads that step's output: the full run on main after the merge is what tests those. The job summary of `e2e-scope` lists every picked case with the reason it was picked. -[`scripts/ci/e2e-case-paths.tsv`](https://github.com/GeiserX/OwnGenome/blob/main/scripts/ci/e2e-case-paths.tsv) holds, for each case, the files it runs: the step scripts, bin helpers, modules (`modules/local//*`) and assets it calls, and what those call in turn. The globs are matched as in bash `[[ PATH == GLOB ]]`, so `*` matches slashes too. A base case's row is `-`. The table was first written from the text of each case, followed through the scripts, helpers and modules it names (comments left out). `e2e-scope.sh check` runs first in `e2e-scope` and fails when a case has no row or two, a row names no case, a base row lists paths, or a glob matches no tracked file; then `e2e-scope` is red and so is `e2e`. A glob too many costs runner time; a glob too few lets a break through to main, where the full run after the merge finds it. +[`scripts/ci/e2e-case-paths.tsv`](https://github.com/GeiserX/OwnGenome/blob/main/scripts/ci/e2e-case-paths.tsv) holds, for each case, the files it runs: the step scripts, bin helpers, modules (`modules/local//*`) and assets it calls, and what those call in turn. The globs are matched as in bash `[[ PATH == GLOB ]]`, so `*` matches slashes too. A base case's row is `-`. The table was first written from the text of each case, followed through the scripts, helpers and modules it names (comments left out). `e2e-scope.sh check` runs first in `e2e-scope` and fails when a case has no row or two, a row names no case, a base row lists paths, a glob matches no tracked file, or a glob matches a file outside the e2e scope (a change to that file would start no run, so the case would never see it); then `e2e-scope` is red and so is `e2e`. A glob too many costs runner time; a glob too few lets a break through to main, where the full run after the merge finds it. -The images the cases use are cached in four GitHub Actions cache entries, one per group. `scripts/ci/e2e-run.sh --image-groups` names each image's group, taken from its repository name without the tag, so a version bump changes the key of one group only; a group whose key missed restores the newest older entry of that group and pulls only what changed. A group keeps its smallest images up to 1500 MB of zstd output, and the larger ones are pulled on every run. Four groups of 1500 MB and the five data entries below stay under the repository's 10 GB cache limit even while main replaces one group: 5.9 GB of images, 1.5 GB more for the group being replaced, and about 2.3 GB of data. Images pinned by digest are not cached, because `docker load` gives them back without the digest, and they are left out of the group's key, so a digest bump misses no entry. Only the base leg saves, and only on main (or in a run by hand with `save_images`, on any branch), after it pulled every image the cases use and got through its cases, and only the groups whose key missed: a cache saved by a pull request can be read by that pull request alone and would only push main's entries out. The job summary lists every image with its size and what the cache did with it: restored on an exact key hit, restored from an older entry, saved, or pulled. Every merge in the e2e scope starts a run on main, so the cache follows main's images. GitHub drops an entry no run has read for 7 days, so after a quiet week the next run on main fills the cache again. Until this layout the cache was one entry, and by October 2026 it had grown to about 11 GB, over the limit, so every run missed it. +The images the cases use are cached in four GitHub Actions cache entries, one per group. `scripts/ci/e2e-run.sh --image-groups` names each image's group, taken from its repository name without the tag, so a version bump changes the key of one group only; a group whose key missed restores the newest older entry of that group and pulls only what changed. A group keeps its smallest images up to 1500 MB of zstd output, and the larger ones are pulled on every run. Four groups of 1500 MB and the six data entries below stay under the repository's 10 GB cache limit even while main replaces one group: 5.9 GB of images, 1.5 GB more for the group being replaced, and about 2.4 GB of data. Images pinned by digest are not cached, because `docker load` gives them back without the digest, and they are left out of the group's key, so a digest bump misses no entry. Only the base leg saves, and only on main (or in a run by hand with `save_images`, on any branch), after it pulled every image the cases use and got through its cases, and only the groups whose key missed: a cache saved by a pull request can be read by that pull request alone and would only push main's entries out. The job summary lists every image with its size and what the cache did with it: restored on an exact key hit, restored from an older entry, saved, or pulled. Every merge in the e2e scope starts a run on main, so the cache follows main's images. GitHub drops an entry no run has read for 7 days, so after a quiet week the next run on main fills the cache again. Until this layout the cache was one entry, and by October 2026 it had grown to about 11 GB, over the limit, so every run missed it. -Five downloads are cached the same way. Two come from EBI's 1000 Genomes server, which is slow and at times stalls: MoChA's chr20 panel with its manifest (0.64 GB, case `mocha-1`), keyed on `MOCHA_PANEL_MANIFEST_SHA256` and the chromosome, and the 1000 Genomes phase 3 sites VCF (1.36 GB, case `archaic-ancestry-1`), keyed on `ARCHAIC_1KG_SITES_MD5`. On a miss, each of those two is tried first at `E2E_DATA_MIRROR`, the GitHub release `e2e-data-1` that holds copies of the pinned files, and then at EBI; the pin checks both. The release also holds the four PGS Catalog score files and the HAPNEST panel below, so an EBI outage costs none of them on a cold cache. Two are small and failed whole runs on 2026-11-04, because their case is in the base leg or a shard that cannot start without it: `hla.dat` of the two IPD-IMGT/HLA releases `bash-step-08-hla` installs, from raw.githubusercontent.com (which served it under 1 KB/s), with GENCODE's gene and exon files from EBI's FTP server (curl exit 7 and 28), and the four real PGS Catalog score files of `prs-scores-and-run-options-1-strand-sex`, from the same server. `scripts/ci/e2e-run.sh --data-key` names their keys: the releases and a hash of the code in `scripts/lib/common.sh` that installs them, and a hash of the score ids and their server. `scripts/ci/e2e-run.sh` hands each case its files just before it runs, at the paths `setup.sh` and `install_data_file` read them from, so they download nothing; the score files the case copies from the cache itself. On a miss it downloads them first, checked against their pins or the checksum files the servers publish beside them (`MD5SUMS`, `md5checksum.txt`, the `.md5` of each score), trying up to 6 times 30 s apart. The four score files are pinned by md5 in `PGS_MD5` of `e2e-run.sh` (the md5 the PGS Catalog served beside each on 2026-11-06), tried on the mirror first and then at EBI; a restored or downloaded score must match its pin as well as its `.md5`. A file gets at most `E2E_MIRROR_BUDGET` seconds on the mirror (900, or half the fill's budget when that is less), so a mirror that trickles leaves EBI its turn. One fill may take at most `E2E_DATA_BUDGET` seconds (2400 in `e2e.yml`): a source that trickles above `fetch`'s stall limit would otherwise run on to the leg's 240-minute limit, as EBI's 1000 Genomes server did twice on 2026-11-05. A fill still running then is stopped and its case fails as `FAIL (data fill over 2400s)` without running, since its own download would be as slow. A restored entry is checked again before use (`--data-check`): the IPD-IMGT/HLA and GENCODE entry against the sha256 list written when it was filled, the score files against their `.md5`; one that fails is downloaded again. `bash-step-08-hla` then runs with every HTTP and HTTPS request sent to a closed port, after a control shows that a request under that block fails at the closed port, so a warm cache cannot reach either host. Only the leg whose cases include the case restores and saves its entry, and it saves only on main or in a run by hand with `save_images`, as with the images. The job summary says whether each key hit, and the case's log says where its files came from. The fifth is the PGS Catalog's synthetic HAPNEST panel (0.27 GB, cases `prs-2-ancestry` and `prs-scores-and-run-options-2-panel-chrx`) from EBI's FTP server, which failed cases with curl exit 7 on 2026-11-04. It is keyed on the sha256 `tests/e2e/lib.sh` pins (`HAPNEST_SHA256`). The cases handle it themselves (`setup_ancestry_panel` in `lib.sh`): with the panel in the cache, or on a miss from `E2E_DATA_MIRROR`, they check the sha256 and run `setup.sh --ancestry-panel` with every HTTP and HTTPS request sent to a closed port, so an EBI outage cannot fail them; when neither has it, `setup.sh` downloads the panel from EBI as a user's run does. The first case puts a panel that did not come from the cache in the cache for the save step. A red control in `prs-2-ancestry` runs `setup.sh` with the network blocked and no panel, which must fail. Every download, there and in `setup.sh`, gives up on a transfer that stalls (`FETCH_STALL`, see `fetch` in `scripts/lib/common.sh`), and the next try resumes it. +Six downloads are cached the same way. Two come from EBI's 1000 Genomes server, which is slow and at times stalls: MoChA's chr20 panel with its manifest (0.64 GB, case `mocha-1`), keyed on `MOCHA_PANEL_MANIFEST_SHA256` and the chromosome, and the 1000 Genomes phase 3 sites VCF (1.36 GB, case `archaic-ancestry-1`), keyed on `ARCHAIC_1KG_SITES_MD5`. On a miss, each of those two is tried first at `E2E_DATA_MIRROR`, the GitHub release `e2e-data-1` that holds copies of the pinned files, and then at EBI; the pin checks both. The release also holds the four PGS Catalog score files, the HAPNEST panel and the table builders' downloads below, so an EBI outage costs none of them on a cold cache. Two are small and failed whole runs on 2026-11-04, because their case is in the base leg or a shard that cannot start without it: `hla.dat` of the two IPD-IMGT/HLA releases `bash-step-08-hla` installs, from raw.githubusercontent.com (which served it under 1 KB/s), with GENCODE's gene and exon files from EBI's FTP server (curl exit 7 and 28), and the four real PGS Catalog score files of `prs-scores-and-run-options-1-strand-sex`, from the same server. `scripts/ci/e2e-run.sh --data-key` names their keys: the releases and a hash of the code in `scripts/lib/common.sh` that installs them, and a hash of the score ids and their server. `scripts/ci/e2e-run.sh` hands each case its files just before it runs, at the paths `setup.sh` and `install_data_file` read them from, so they download nothing; the score files the case copies from the cache itself. On a miss it downloads them first, checked against their pins or the checksum files the servers publish beside them (`MD5SUMS`, `md5checksum.txt`, the `.md5` of each score), trying up to 6 times 30 s apart. The four score files are pinned by md5 in `PGS_MD5` of `e2e-run.sh` (the md5 the PGS Catalog served beside each on 2026-11-06), tried on the mirror first and then at EBI; a restored or downloaded score must match its pin as well as its `.md5`. A file gets at most `E2E_MIRROR_BUDGET` seconds on the mirror (900, or half the fill's budget when that is less), so a mirror that trickles leaves EBI its turn. One fill may take at most `E2E_DATA_BUDGET` seconds (2400 in `e2e.yml`): a source that trickles above `fetch`'s stall limit would otherwise run on to the leg's 240-minute limit, as EBI's 1000 Genomes server did twice on 2026-11-05. A fill still running then is stopped and its case fails as `FAIL (data fill over 2400s)` without running, since its own download would be as slow. A restored entry is checked again before use (`--data-check`): the IPD-IMGT/HLA and GENCODE entry against the sha256 list written when it was filled, the score files against their `.md5`; one that fails is downloaded again. `bash-step-08-hla` then runs with every HTTP and HTTPS request sent to a closed port, after a control shows that a request under that block fails at the closed port, so a warm cache cannot reach either host. Only the leg whose cases include the case restores and saves its entry, and it saves only on main or in a run by hand with `save_images`, as with the images. The job summary says whether each key hit, and the case's log says where its files came from. The fifth is the PGS Catalog's synthetic HAPNEST panel (0.27 GB, cases `prs-2-ancestry` and `prs-scores-and-run-options-2-panel-chrx`) from EBI's FTP server, which failed cases with curl exit 7 on 2026-11-04. It is keyed on the sha256 `tests/e2e/lib.sh` pins (`HAPNEST_SHA256`). The cases handle it themselves (`setup_ancestry_panel` in `lib.sh`): with the panel in the cache, or on a miss from `E2E_DATA_MIRROR`, they check the sha256 and run `setup.sh --ancestry-panel` with every HTTP and HTTPS request sent to a closed port, so an EBI outage cannot fail them; when neither has it, `setup.sh` downloads the panel from EBI as a user's run does. The first case puts a panel that did not come from the cache in the cache for the save step. A red control in `prs-2-ancestry` runs `setup.sh` with the network blocked and no panel, which must fail. The sixth, entry `catalogs` (0.08 GB), holds the pinned downloads of `scripts/ci/build-callability.sh` and `scripts/ci/build-cn-intervals.py`, whose cases `callability-1-catalog` and `cn-intervals-1-catalog` rebuild `assets/callability_genes.tsv` and the interval rows of `bin/cn_loci.tsv` and compare them with the committed files. Each builder lists its downloads with `--list`; the key is a hash of those lists, so a new pin is a new key. GENCODE's GTF among them comes from EBI's FTP server, which refused every connection for minutes on 2026-11-06 and failed `callability-1-catalog` on main (run 42259432681). The runner fills the entry as it does the others, mirror first and every file checked against its sha256, and each case points its builder's cache at its folder there; a builder that still has to download tries for about five minutes, 30 s apart. Every download, there and in `setup.sh`, gives up on a transfer that stalls (`FETCH_STALL`, see `fetch` in `scripts/lib/common.sh`), and the next try resumes it. How long it took in one job, before the split: a passing e2e job took about 105 to 205 minutes with 84 cases, and 230 minutes of cases with 130 (run 41390434697, 2026-11-02). The runner decides much of it, not the cache: the same CPU model took 106 minutes in one run and 166 in another. In October 2026, with 84 cases, the same code took 106 and 166 minutes on an AMD EPYC 9V45 (runs 39197502832 and 39075176392, the second of which also saved the cache), 112 minutes on an Intel Xeon 6973P (run 39088086725), 156 and 191 minutes on an AMD EPYC 7763 (runs 39075764278 and 39075201465) and 205 minutes on an AMD EPYC 9V74 (run 39197512594, which saved three groups). The image cache saves about 7 minutes of that: pulling every image took 9 minutes, restoring and loading the four groups took 2.5 (run 39075201465 restored 39 images in 118 seconds), and saving them all took 12. Each leg has a limit of 240 minutes (355 for a traced shard leg); the base leg waits up to 75 minutes for a fixture a pull request bumps. Each leg's first steps print the runner's CPU model and CPU count to the log and the job summary, so a later timing can name the hardware. @@ -121,10 +121,14 @@ What it covers: - **Index times** (`run-all-sample-qc-and-stale-index-1-index-times.sh`): every index the direct `nextflow run` and the launcher published is not older than its data file, and samtools and bcftools print no stale-index warning on the BAM and the VCF. Then a short VCF-only run with three planted pairs checks `main.nf`'s completion handler: it sets the index times of the run's own files after the last copy, and leaves an index older than the run as it was. A copied pair with its index set an hour older is the negative control. - **The Nextflow report's tables** (`paralogs-report-section-1-nf-report-tables.sh`): `HTML_REPORT` runs for real on the CYP2D6 depth check, the cDNA read check and the phased compound-het pairs that earlier cases wrote, and its reports must show the depth-check line, the cDNA column with its flagged count and the pairs by phase; the same run without the three tables must show none of them. A stub run of every tool then checks that `main.nf` passes the cDNA table to each `HTML_REPORT` task that got ClinVar hits, against a copy of the pipeline with that join removed. +- **Locityper, step 53** (`ci-e2e-2-locityper.sh`): `setup.sh --locityper-data` builds the database for the HLA-DQB1 row of `assets/locityper_loci.bed` from a chr6 and chr20 reference cut from the fixture's, and `scripts/53-locityper.sh` genotypes the fixture's reads aligned to it, with chr20 as the background region. The fixture has no reads over FCGR or KIR, and HLA-DRB1 is left out for time: with it the case took 59 minutes, 35 of them jellyfish's k-mer count and 23 `locityper augment`. With HLA-DQB1 alone `augment` takes about a minute, but jellyfish, which does not depend on the loci, still takes 36 (38 minutes in all on an AMD EPYC 9V74), so the case raises its own limit to 90 minutes (`# case-timeout: 5400`; every other case keeps 60). The locus must come out called, with two panel haplotypes, a quality, reads and the authors' filter, in the columns `bin/locityper_report.py` writes. Red controls: the step stops without the database, without the BAM's index, and with a database whose `DONE` names another reference. + - **The 1.0 output contract** (`output-contract-1-hg002.sh`): every stable file of [`assets/output_contract.tsv`](https://github.com/GeiserX/OwnGenome/blob/main/assets/output_contract.tsv) whose `e2e` column names a leg is in what the base cases' bash steps and case `60-nextflow` wrote, under its contracted name, with its required columns, VCF fields and JSON keys. Its red control drops a column and renames a file on a copy and must fail on both. `tests/test_output_contract.py` checks the same contract against the code with no run ([What 1.0 keeps stable](output-contract.md)). - **Report assets:** a last case lists every external `http(s)` address in a `src=` or `href=` attribute of the generated HTML reports, so the remote files a report loads when opened are known. The list goes to the job summary. +A case that exits 0 must have run at least one check: `finish` in `lib.sh` fails a case that ran none, and `e2e-run.sh` fails one whose log has no `[PASS]` line (`FAIL (no check ran)`), so a case that loses its checks to an early `finish` cannot stay green. A case with nothing to do on a run sets `NO_CHECKS_REASON` before `finish`, which prints the reason and adds it to the job summary; `prs-3-panel-measure` does that on every run that does not measure the panel. `tests/test_e2e_finish.sh` checks `finish` on planted cases. + Each check is a count or a column: a VCF with records, the sample name in the VCF header, a gene in the ClinVar hit, at least one called gene in PharmCAT's `report.json`, a CYP2D6 row from pypgx, the planted ClinVar row in the HTML report with its gene, its significance and no empty cell. An exit code alone never passes a case. Many steps ask Docker for `--cpus 8`, which Docker refuses on a 4-CPU machine. The job puts a small shim first on `PATH` ([`tests/e2e/bin/docker`](https://github.com/GeiserX/OwnGenome/blob/main/tests/e2e/bin/docker)) that lowers `--cpus` to the CPUs the machine has and passes everything else through. diff --git a/modules/local/annotsv/main.nf b/modules/local/annotsv/main.nf index e2fc74c4..b4de6897 100644 --- a/modules/local/annotsv/main.nf +++ b/modules/local/annotsv/main.nf @@ -25,7 +25,8 @@ 5 rows are written to their own table, found by column name, and bin/annotsv_postfilter.py writes each of them with its caller, whether a deletion sits on the exon edges of one transcript of AnnotSV's genes - file, its benign-set frequency, its read-depth support and the class + file, its benign-set frequency, its read-depth support (a Delly RDCN + against the same sex- and PAR-aware normal copies) and the class after the rule table assets/sv_rule_regions.tsv, which moves to class 3 a call inside a coding tandem repeat, the IG and TR loci or the MHC, one the depth contradicts, and one whose gene symbol is a pseudogene's @@ -101,7 +102,7 @@ process ANNOTSV { genes_arg=() if [ -s "\$genes" ]; then genes_arg=(--genes "\$genes"); fi annotsv_postfilter.py --class45 ${meta.id}_sv_class4_5.tsv \${genes_arg[@]+"\${genes_arg[@]}"} \\ - --rules ${sv_rules} --out ${meta.id}_sv_class4_5_checked.tsv + --rules ${sv_rules} --sex '${meta.sex ?: ''}' --par-bed ${par_bed} --out ${meta.id}_sv_class4_5_checked.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/renovate.json b/renovate.json index 4cf6b687..76858499 100644 --- a/renovate.json +++ b/renovate.json @@ -63,6 +63,14 @@ "extractVersionTemplate": "^v(?.+)$", "versioningTemplate": "regex:^(?\\d+)\\.(?\\d+)\\.(?\\d+)$" }, + { + "customType": "regex", + "description": "nf-weblog, the Nextflow plugin of the live site: NF_WEBLOG_PLUGIN in scripts/run-all.sh and the plugins block docs/live-site.md shows, in one PR. Its releases are tagged without a v.", + "managerFilePatterns": ["/^scripts/run-all\\.sh$/", "/^docs/live-site\\.md$/"], + "matchStrings": ["NF_WEBLOG_PLUGIN=nf-weblog@(?\\d+\\.\\d+\\.\\d+)", "id 'nf-weblog@(?\\d+\\.\\d+\\.\\d+)'"], + "depNameTemplate": "nextflow-io/nf-weblog", + "datasourceTemplate": "github-releases" + }, { "customType": "regex", "description": "actionlint, downloaded by version in the lint workflow.", diff --git a/scripts/05-annotsv.sh b/scripts/05-annotsv.sh index bfbcec2e..ec4638a6 100755 --- a/scripts/05-annotsv.sh +++ b/scripts/05-annotsv.sh @@ -174,7 +174,8 @@ fi # caller, whether a deletion sits on the exon edges of one transcript of # AnnotSV's genes file (a processed transcript aligned across an intron), the # benign-set frequency of its type, and whether its read depth shows the -# change. The rule table moves to class 3 a call inside a coding tandem +# change (a Delly RDCN against the same sex- and PAR-aware normal copies as +# above). The rule table moves to class 3 a call inside a coding tandem # repeat, the IG and TR loci or the MHC, one the depth contradicts, and one # whose gene symbol is a pseudogene's alias. How many callers saw it comes # from step 22's consensus, which the reports join when they are written. @@ -187,9 +188,10 @@ else echo "Exon edges: no ${GENES_BED#"$GENOME_DIR"/} in the annotation data, so no deletion is checked against them." fi run_in --cpus 1 --memory 2g -v "${PGP_ROOT}/bin:/pgp-bin:ro" \ - -v "${RULES_TSV}:/pgp/sv_rule_regions.tsv:ro" "${ANNOTSV_IMAGE}" \ + -v "${RULES_TSV}:/pgp/sv_rule_regions.tsv:ro" \ + -v "${PGP_ROOT}/assets/par_grch38.bed:/pgp/par_grch38.bed:ro" "${ANNOTSV_IMAGE}" \ python3 /pgp-bin/annotsv_postfilter.py --class45 "$(cpath "$CLASS45")" ${GENES_ARGS[@]+"${GENES_ARGS[@]}"} \ - --rules /pgp/sv_rule_regions.tsv --out "$(cpath "$CHECKED")" + --rules /pgp/sv_rule_regions.tsv --sex "$SEX" --par-bed /pgp/par_grch38.bed --out "$(cpath "$CHECKED")" echo "=== AnnotSV complete ===" echo "Input: ${INPUT_VCF}" diff --git a/scripts/06-clinvar-screen.sh b/scripts/06-clinvar-screen.sh index 7023c861..4d938d89 100755 --- a/scripts/06-clinvar-screen.sh +++ b/scripts/06-clinvar-screen.sh @@ -61,6 +61,18 @@ done echo "=== ClinVar Pathogenic Screen: ${SAMPLE} ===" +# bin/owngenome_recheck.py --refresh-clinvar holds this lock while it replaces the ClinVar files every sample +# shares; the screen holds it shared on fd 7 until it ends, so it waits for a refresh going on and a refresh +# refuses while it reads. Without flock (macOS: brew install flock) nothing waits here. +CVLOCK="${GENOME_DIR}/clinvar/.refresh.lock" +if [ -d "${GENOME_DIR}/clinvar" ] && command -v flock >/dev/null 2>&1; then + [ -e "$CVLOCK" ] || : 2>/dev/null >> "$CVLOCK" || true + if [ -r "$CVLOCK" ]; then + exec 7<"$CVLOCK" + flock -s -n 7 || { echo "Waiting for the ClinVar refresh that holds ${CVLOCK} to end ..."; flock -s 7; } + fi +fi + for f in "$VCF" "${VCF}.tbi" "$CLINVAR" "${CLINVAR}.tbi" "$REF" "${REF}.fai" "$SITES"; do if [ ! -f "$f" ]; then echo "ERROR: File not found: ${f}" >&2 diff --git a/scripts/ci/build-callability.sh b/scripts/ci/build-callability.sh index 05927b0d..e4705848 100755 --- a/scripts/ci/build-callability.sh +++ b/scripts/ci/build-callability.sh @@ -5,6 +5,7 @@ # not read). # # Usage: scripts/ci/build-callability.sh [OUT] (default: assets/callability_genes.tsv) +# scripts/ci/build-callability.sh --list the downloads, one "NAME URL SHA256" per line # # Downloads seven pinned files (about 85 MB, each checked against its # sha256) into a temporary folder, or reads them from CALLABILITY_CACHE when @@ -72,6 +73,14 @@ FILES=( "report.txt|https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/GCA_000001405.15_GRCh38_assembly_report.txt|aa733ce92719f6339b3a4540137d4322894434363314c4b7704ba6691e9eab66" ) +if [ "${1:-}" = --list ]; then + for f in "${FILES[@]}"; do + IFS='|' read -r name url sum <<<"$f" + echo "${name} ${url} ${sum}" + done + exit 0 +fi + sha256() { if command -v sha256sum >/dev/null 2>&1; then sha256sum "$1" | cut -d' ' -f1; else shasum -a 256 "$1" | cut -d' ' -f1; fi } @@ -88,7 +97,9 @@ get() { return fi echo "Downloading $2" >&2 - curl -fsSL --retry 3 --retry-all-errors -o "${WORK}/$1" "$2" + # Tries spread over about five minutes: EBI's FTP server refused every + # connection for minutes at a time (curl exit 7, 2026-11-06). + curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "${WORK}/$1" "$2" got=$(sha256 "${WORK}/$1") if [ "$got" != "$3" ]; then echo "ERROR: $1 has sha256 ${got}, want $3" >&2 diff --git a/scripts/ci/build-cn-intervals.py b/scripts/ci/build-cn-intervals.py index c676deab..6c186bf7 100755 --- a/scripts/ci/build-cn-intervals.py +++ b/scripts/ci/build-cn-intervals.py @@ -49,6 +49,7 @@ Usage: scripts/ci/build-cn-intervals.py [--catalog bin/cn_loci.tsv] [--cache DIR] [--check] + scripts/ci/build-cn-intervals.py --list the downloads, one "NAME URL SHA256" per line Rewrites the catalog's rows from the marker line on (MARKER), or with --check exits 1 when the catalog is not what the downloads give. The downloads (about @@ -62,6 +63,7 @@ import re import sys import tempfile +import time import urllib.request REPO = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) @@ -140,8 +142,10 @@ def fetch(name, url, want, folder): os.replace(path + ".part", path) break except OSError as e: - if i == 2: + # A missing or forbidden file stays so: no point waiting. + if i == 2 or getattr(e, "code", None) in (403, 404): raise BuildError(f"cannot download {url}: {e}") + time.sleep(30) got = sha256(path) if got != want: raise BuildError(f"{url}: sha256 {got}, pinned {want}") @@ -369,7 +373,20 @@ def main(argv=None): ap.add_argument("--catalog", default=os.path.join(REPO, "bin", "cn_loci.tsv")) ap.add_argument("--cache") ap.add_argument("--check", action="store_true") + ap.add_argument("--list", action="store_true") a = ap.parse_args(argv) + if a.list and a.check: + ap.error("--list and --check do not go together") + if a.list: + try: + files = [("clingen_recurrent_cnv.bed", versions_env("CLINGEN_RECURRENT_CNV_URL"), + versions_env("CLINGEN_RECURRENT_CNV_SHA256"))] + MASKS + except BuildError as e: + print(f"ERROR: {e}", file=sys.stderr) + return 2 + for n, u, h in files: + print(n, u, h) + return 0 try: with tempfile.TemporaryDirectory() as tmp: folder = a.cache or tmp diff --git a/scripts/ci/build-sv-rule-regions.sh b/scripts/ci/build-sv-rule-regions.sh index 6edc6bf2..4c92fc27 100755 --- a/scripts/ci/build-sv-rule-regions.sh +++ b/scripts/ci/build-sv-rule-regions.sh @@ -3,9 +3,16 @@ # where bin/annotsv_postfilter.py overrides AnnotSV's class 4 or 5. # # Usage: scripts/ci/build-sv-rule-regions.sh [OUT] (default: assets/sv_rule_regions.tsv) +# scripts/ci/build-sv-rule-regions.sh --check [FILE] exit 1 when FILE +# (default: assets/sv_rule_regions.tsv) is not what the downloads +# give, byte for byte, and show the first lines that differ +# scripts/ci/build-sv-rule-regions.sh --list the downloads, one +# "NAME URL SHA256" per line # # Downloads four pinned files (about 145 MB, each checked against its sha256) -# into a temporary folder, writes OUT and deletes the downloads: +# into a temporary folder, or reads them from SV_RULE_CACHE when that folder +# holds them with the right sha256 (and keeps what it downloads there), writes +# OUT and deletes the downloads: # # TRExplorer v2.1's GRCh38 tandem repeat catalog (broadinstitute/trexplorer-catalog) # GENCODE 50's basic annotation (the release scripts/lib/common.sh pins) @@ -33,7 +40,12 @@ set -euo pipefail REPO=$(cd "$(dirname "$0")/../.." && pwd) -OUT=${1:-${REPO}/assets/sv_rule_regions.tsv} +CHECK="" +case "${1:-}" in + --list) ;; + --check) CHECK=${2:-${REPO}/assets/sv_rule_regions.tsv} ;; + *) OUT=${1:-${REPO}/assets/sv_rule_regions.tsv} ;; +esac TRX_URL=https://github.com/broadinstitute/trexplorer-catalog/releases/download/v2.1/TRExplorer.repeat_catalog_v2.1.hg38.1_to_1000bp_motifs.bed.gz TRX_SHA256=4eb0fbbd36d40a801c2d764cdda247fea8edfae1604087526459064624ec1417 @@ -44,23 +56,42 @@ HGNC_SHA256=8bf6f686e640dbc864a346484626953621b0500d900d26cd085a35c2782d725d REGIONS_URL=https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/GCA_000001405.15_GRCh38_assembly_regions.txt REGIONS_SHA256=6b91d2c7961a322936db8fd09b7fb4a9143a590882f28c0a628a4ef26897f6c7 +if [ "${1:-}" = --list ]; then + echo "trx.bed.gz ${TRX_URL} ${TRX_SHA256}" + echo "gencode.gtf.gz ${GENCODE_URL} ${GENCODE_SHA256}" + echo "hgnc.txt ${HGNC_URL} ${HGNC_SHA256}" + echo "regions.txt ${REGIONS_URL} ${REGIONS_SHA256}" + exit 0 +fi + sha256() { if command -v sha256sum >/dev/null 2>&1; then sha256sum "$1" | cut -d' ' -f1; else shasum -a 256 "$1" | cut -d' ' -f1; fi } WORK=$(mktemp -d) trap 'rm -rf "$WORK"' EXIT +[ -z "$CHECK" ] || OUT="${WORK}/rebuilt.tsv" -# get NAME URL SHA256: download URL to WORK/NAME and stop unless it matches. +# get NAME URL SHA256: WORK/NAME from SV_RULE_CACHE when it matches, else +# downloaded; stop unless it matches. get() { + local cached="${SV_RULE_CACHE:-}/$1" got + if [ -n "${SV_RULE_CACHE:-}" ] && [ -f "$cached" ] && [ "$(sha256 "$cached")" = "$3" ]; then + cp "$cached" "${WORK}/$1" + return + fi echo "Downloading $2" >&2 - curl -fsSL --retry 3 -o "${WORK}/$1" "$2" - local got + # Tries spread over about five minutes: EBI's FTP server refused every + # connection for minutes at a time (curl exit 7, 2026-11-06). + curl -fsSL --retry 6 --retry-all-errors --retry-delay 30 --retry-max-time 300 -o "${WORK}/$1" "$2" got=$(sha256 "${WORK}/$1") if [ "$got" != "$3" ]; then echo "ERROR: $1 has sha256 ${got}, want $3" >&2 exit 1 fi + if [ -n "${SV_RULE_CACHE:-}" ]; then + mkdir -p "$SV_RULE_CACHE" && cp "${WORK}/$1" "$cached" + fi } get trx.bed.gz "$TRX_URL" "$TRX_SHA256" get gencode.gtf.gz "$GENCODE_URL" "$GENCODE_SHA256" @@ -188,4 +219,13 @@ with open(out, "w") as f: print("Rows: " + ", ".join(f"{k} {counts[k]}" for k in sorted(counts)), file=sys.stderr) PY mv "${OUT}.part" "$OUT" +if [ -n "$CHECK" ]; then + if cmp -s "$OUT" "$CHECK"; then + echo "${CHECK} is what the pinned downloads give, byte for byte" >&2 + exit 0 + fi + echo "ERROR: ${CHECK} is not what the pinned downloads give; rebuild it with scripts/ci/build-sv-rule-regions.sh. The first differences (< the file, > the rebuild):" >&2 + diff "$CHECK" "$OUT" | head -n 20 >&2 || true + exit 1 +fi echo "Wrote ${OUT}" >&2 diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 1254dfcb..c302ec67 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -65,6 +65,7 @@ bumps-drop-in-baseline scripts/ci/fake-docker/lib.sh callability-1-catalog assets/callability_genes.tsv assets/paralog_loci.tsv assets/pgx_genes.tsv bin/cn_loci.tsv scripts/ci/build-callability.sh scripts/ci/fake-docker/lib.sh callability-2-genes assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/stub/sample.vcf.gz bin/collect_summary.py bin/render_report.py bin/results_index.py modules/local/html_report/* modules/local/mosdepth/* modules/local/vcf_precheck/* scripts/16b-mosdepth.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh chip-clonal-1 assets/chip_genes.tsv bin/chip_report.py bin/render_report.py modules/local/chip/* modules/local/vcf_precheck/* scripts/49-chip.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh +ci-e2e-2-locityper assets/locityper_loci.bed bin/locityper_report.py scripts/53-locityper.sh scripts/setup.sh clinical-filter-mane-and-quality-1-tiers assets/construct_variants.tsv bin/clip_check.py bin/comphet_phase.py bin/constraint_join.awk bin/csq_rank.awk scripts/23-clinical-filter.sh scripts/31-slivar.sh scripts/ci/fake-docker/lib.sh clinical-filter-mane-and-quality-2-vep-norm scripts/13-vep-annotation.sh scripts/ci/fake-docker/lib.sh clinical-filter-mane-and-quality-3-modules assets/construct_variants.tsv assets/stub/sample.bam assets/stub/sample.bam.bai bin/clip_check.py bin/comphet_phase.py bin/constraint_join.awk bin/csq_rank.awk modules/local/clinical_filter/* modules/local/slivar/* modules/local/vcf_precheck/* modules/local/vep/* scripts/13-vep-annotation.sh scripts/23-clinical-filter.sh scripts/31-slivar.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh @@ -119,7 +120,7 @@ paralog-checks-and-doc-fixes-1-hba assets/par_grch38.bed assets/paralog_loci.tsv paralog-rows-1 assets/paralog_loci.tsv assets/parascopy_models_v1.7.txt bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/paraphase_report.py bin/render_report.py scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh paralogs-paraphase-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/stub/sample.vcf.gz bin/collect_summary.py bin/paraphase_report.py bin/render_report.py bin/results_index.py modules/local/align_minimap2/* modules/local/deepvariant/* modules/local/fastp/* modules/local/html_report/* modules/local/paraphase/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/24-html-report.sh scripts/35b-paraphase.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh paralogs-report-section-1-nf-report-tables assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/collect_summary.py bin/comphet_phase.py bin/constraint_join.awk bin/copy_number.py bin/csq_rank.awk bin/cyp2d6_depth_check.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py bin/results_index.py modules/local/clinvar_screen/* modules/local/export_consumer/* modules/local/html_report/* modules/local/pypgx/* modules/local/slivar/* modules/local/vcf_precheck/* scripts/06-clinvar-screen.sh scripts/24-html-report.sh scripts/31-slivar.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz -pgx-1-depth-check bin/cyp2d6_depth_check.py scripts/ci/alt-depth-ab.sh scripts/ci/fake-docker/lib.sh +pgx-1-depth-check bin/cyp2d6_depth_check.py scripts/21-cyrius.sh scripts/32-pypgx.sh scripts/ci/alt-depth-ab.sh scripts/ci/fake-docker/lib.sh pgx-2-outside-calls assets/pgx_contraindications.tsv assets/pgx_genes.tsv assets/pgx_inhibitors.tsv bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_parse.py bin/pgx_refblock_fill.py scripts/07-pharmacogenomics.sh scripts/27-cpic-lookup.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh pgx-4-parascopy assets/paralog_loci.tsv bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh pgx-5-kir assets/pgx_genes.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat bin/cn_loci.tsv bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py scripts/08-hla-typing.sh scripts/24-html-report.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index 6bf12c02..82e0a424 100644 --- a/scripts/ci/e2e-case-times.tsv +++ b/scripts/ci/e2e-case-times.tsv @@ -64,6 +64,7 @@ bumps-drop-in-baseline shard 4 bash-step-11-roh,bash-step-14-imputation,bash-ste callability-1-catalog shard 74 - callability-2-genes shard 32 - chip-clonal-1 shard 277 - +ci-e2e-2-locityper shard 2289 - clinical-filter-mane-and-quality-1-tiers shard 14 - clinical-filter-mane-and-quality-2-vep-norm shard 2 - clinical-filter-mane-and-quality-3-modules shard 8 clinical-filter-mane-and-quality-1-tiers diff --git a/scripts/ci/e2e-run.sh b/scripts/ci/e2e-run.sh index c0496df1..2546b6aa 100755 --- a/scripts/ci/e2e-run.sh +++ b/scripts/ci/e2e-run.sh @@ -14,7 +14,8 @@ # --pull-images pulls those images, as step 3 below does, and exits. # --pull-image REF pulls one image with the same retry (the docker shim's # pulls of the lazy images at their first use), and exits. -# --data-key NAME prints the Actions cache key of data entry NAME (hla, pgs). +# --data-key NAME prints the Actions cache key of data entry NAME (hla, pgs, +# catalogs). # --data-check NAME exits 0 when entry NAME in E2E_DATA_CACHE is whole and # matches its checksums (the e2e workflow saves only such an entry). # --data-fill NAME fills entry NAME in E2E_DATA_CACHE as a case's run does, @@ -173,8 +174,8 @@ if [ "${1:-}" = --image-groups ]; then fi # --- Data cache keys and checks ------------------------------------------------- -# The e2e workflow asks here for the Actions cache key of the entries hla and -# pgs (`--data-key NAME`) and whether an entry is whole before it saves it +# The e2e workflow asks here for the Actions cache key of the entries hla, +# pgs and catalogs (`--data-key NAME`) and whether an entry is whole before it saves it # (`--data-check NAME`), so the key and the check are written once. The other # entries (mocha, archaic) are keyed in e2e.yml on their versions.env pins. DATA_CACHE=${E2E_DATA_CACHE:-} @@ -187,6 +188,23 @@ PGS_MD5=${E2E_PGS_MD5-"PGS002257 bb1b2b11a499c5930e8c835647e48296 PGS000065 90ac7aabbbde9bdaa681dbdb66c0ebd8 PGS000064 43b39c4450f423f0772bbd8779dd640f PGS000066 b34f31b8b3a7dbbbb58882a70e14649d"} +# CATALOGS: the builders of the committed tables two cases rebuild and +# compare, as "FOLDER COMMAND...": each lists its pinned downloads with +# --list ("NAME URL SHA256"), and entry catalogs holds them in FOLDER under +# NAME, where the builder's cache reads them (callability-1-catalog, +# cn-intervals-1-catalog). GENCODE's GTFs come from +# EBI's FTP server, which refused every connection on 2026-11-06 (run +# 42259432681). +CATALOGS="callability bash ${REPO}/scripts/ci/build-callability.sh +cn-intervals python3 ${REPO}/scripts/ci/build-cn-intervals.py" +# catalog_files: "FOLDER NAME URL SHA256" of every download of CATALOGS. +catalog_files() { + local dir interp builder out + while read -r dir interp builder; do + out=$("$interp" "$builder" --list) && [ -n "$out" ] || { echo "ERROR: ${builder} --list gave nothing" >&2; return 1; } + sed "s|^|${dir} |" <<< "$out" + done <<< "$CATALOGS" +} # hla_releases: the IPD-IMGT/HLA releases bash-step-08-hla installs (the # release before the pin, then the pin), from the case and scripts/lib/common.sh. hla_releases() { @@ -211,6 +229,9 @@ data_key() { pgs) h=$(grep -E '^(NEW_IDS|PGS_BASE)=' "$PGS_CASE" | sha256sum | cut -c1-16) echo "e2e-data-pgs-${h}" ;; + catalogs) + h=$(catalog_files | sha256sum | cut -c1-16) && [ "${#h}" -eq 16 ] || return 1 + echo "e2e-data-catalogs-${h}" ;; *) echo "ERROR: no data cache key for ${1}" >&2; return 1 ;; esac } @@ -234,6 +255,13 @@ data_ok() { && { [ -z "$pin" ] || [ "$sum" = "$pin" ]; } \ || return 1 done ;; + catalogs) + # Every download of every builder, matching its pin. + local list dir name url + list=$(catalog_files) || return 1 + while read -r dir name url pin; do + [ -s "${d}/${dir}/${name}" ] && [ "$(sha256sum "${d}/${dir}/${name}" | cut -d' ' -f1)" = "$pin" ] || return 1 + done <<< "$list" ;; *) return 1 ;; esac } @@ -258,10 +286,14 @@ data_ok() { # same way, so the mirror can only save time, never change what a case reads. # kill_tree PID: PID and every process under it (data_fetch's subshells and # their curl). +# The children are listed first and the parent is stopped before them: a +# parent still running when its child dies would start the next command +# (mfetch's source download after the mirror try), which the list missed. kill_tree() { - local c - for c in $(pgrep -P "$1" 2> /dev/null); do kill_tree "$c"; done + local c kids + kids=$(pgrep -P "$1" 2> /dev/null || true) kill -TERM "$1" 2> /dev/null || true + for c in $kids; do kill_tree "$c"; done } # within SECS CMD...: runs CMD (a function too) and stops it with everything # under it after SECS seconds; returns CMD's status, or 124 when it was stopped. @@ -380,6 +412,15 @@ data_fetch() ( fi done data_ok pgs || exit 1 ;; + catalogs) + list=$(catalog_files) || exit 1 + while read -r dir name url sum; do + f="${d}/${dir}/${name}" + [ -s "$f" ] && [ "$(_digest sha256 "$f")" = "$sum" ] && continue + rm -f "$f" + mfetch "$url" "$f" sha256 "$sum" || exit 1 + done <<< "$list" + data_ok catalogs || exit 1 ;; *) echo "ERROR: no data cache entry ${1}" >&2; exit 1 ;; esac ) @@ -677,6 +718,8 @@ trap 'kill "$LOW_PID" 2> /dev/null || true' EXIT # finds them and downloads nothing. # prs-scores-and-run-options-1: the four real PGS Catalog score files # (ftp.ebi.ac.uk); the case copies them from E2E_DATA_CACHE/pgs. +# callability-1-catalog, cn-intervals-1-catalog: the pinned downloads of their builders (CATALOGS above); each case +# points its builder's cache at E2E_DATA_CACHE/catalogs/. data_name() { # shellcheck source=../../versions.env case "$1" in @@ -684,6 +727,7 @@ data_name() { archaic-ancestry-1.sh) echo "archaic ${GENOME_DIR}/reference/archaic/$(. "${REPO}/versions.env" && echo "$ARCHAIC_DATA_VERSION")/download" ;; bash-step-08-hla.sh) echo "hla ${GENOME_DIR} offline" ;; prs-scores-and-run-options-1-strand-sex.sh) echo "pgs -" ;; + callability-1-catalog.sh|cn-intervals-1-catalog.sh) echo "catalogs -" ;; esac } # Every HTTP(S) request of a process run under it goes to a closed port. @@ -820,6 +864,10 @@ for c in "${CASES[@]}"; do [[ "$c" == $PATTERN ]] || continue fi log="${LOG_DIR}/${name}.log" + # A case may raise its own limit with a line "# case-timeout: SECONDS" + # (never below CASE_TIMEOUT). + ctimeout=$(sed -n 's/^# case-timeout: \([0-9][0-9]*\)$/\1/p' "${REPO}/tests/e2e/${c}" | head -n 1) + [ "${ctimeout:-0}" -gt "$CASE_TIMEOUT" ] || ctimeout=$CASE_TIMEOUT start=$(date +%s) tstart=$(date +%s.%N) echo "::group::${name}" @@ -829,7 +877,7 @@ for c in "${CASES[@]}"; do rc=-2 else set +e - timeout -k 30 "$CASE_TIMEOUT" "${CASE_ENV[@]}" bash "${REPO}/tests/e2e/${c}" 2>&1 | tee "$log" + timeout -k 30 "$ctimeout" "${CASE_ENV[@]}" bash "${REPO}/tests/e2e/${c}" 2>&1 | tee "$log" rc=${PIPESTATUS[0]} set -e fi @@ -843,11 +891,15 @@ for c in "${CASES[@]}"; do echo "::endgroup::" # A case whose network block could not be shown to hold proves nothing offline. [ "$rc" -eq 0 ] && $CASE_CONTROL_FAILED && rc=-1 + # Nor does a case that exited 0 without one [PASS] line and without saying + # why it had nothing to check (NO_CHECKS_REASON, finish in lib.sh). + [ "$rc" -eq 0 ] && ! grep -q '^\[PASS\] ' "$log" && ! grep -q "^${name}: no checks run: " "$log" && rc=-3 case "$rc" in 0) result=pass ;; -1) result="FAIL (network-block control)"; failed=$((failed + 1)) ;; + -3) result="FAIL (no check ran)"; failed=$((failed + 1)) ;; -2) result="FAIL (data fill over ${E2E_DATA_BUDGET:-2400}s)"; failed=$((failed + 1)) ;; - 124) result="FAIL (timeout ${CASE_TIMEOUT}s)"; failed=$((failed + 1)) ;; + 124) result="FAIL (timeout ${ctimeout}s)"; failed=$((failed + 1)) ;; *) result="FAIL (exit ${rc})"; failed=$((failed + 1)) ;; esac secs=$(( $(date +%s) - start )) diff --git a/scripts/ci/e2e-scope.sh b/scripts/ci/e2e-scope.sh index 0793b294..4e254dc5 100755 --- a/scripts/ci/e2e-scope.sh +++ b/scripts/ci/e2e-scope.sh @@ -8,7 +8,8 @@ # # check exits 1 unless every case of tests/e2e has exactly one row in # e2e-case-paths.tsv, every row names a case, a base case's row is -# "-" and every glob matches a tracked file. It names each fault. +# "-" and every glob matches a tracked file, and only files inside +# the e2e scope. It names each fault. # scope prints the e2e scope as the push paths of e2e.yml are written # (a directory as DIR/**, then the single files), sorted. # rows reads `git diff -U0` output of the two case tables @@ -42,7 +43,7 @@ export LC_ALL=C # The paths the e2e cases exercise. A pull request that changes none of them # runs no case. The push paths of e2e.yml are the same list (e2e-plan checks). -SCOPE_DIRS=(scripts modules workflows bin conf assets tests/e2e) +SCOPE_DIRS=(scripts modules workflows bin conf assets tests/e2e .github/actions) SCOPE_FILES=(tests/fixtures/VERSION main.nf nextflow.config nextflow_schema.json versions.env .github/workflows/e2e.yml # The committed fixtures the cases read (their rows name them too). tests/fixtures/annotsv/class45.sh tests/fixtures/comphet/make_fixture.py tests/fixtures/consumer_export/sites.tsv.gz @@ -59,7 +60,7 @@ FULL=(main.nf nextflow.config nextflow_schema.json versions.env assets/steps.tsv scripts/02-alignment.sh scripts/03-deepvariant.sh tests/fixtures/VERSION scripts/ci/build-fixture.sh scripts/ci/e2e-run.sh scripts/ci/e2e-shards.sh scripts/ci/e2e-scope.sh - tests/e2e/lib.sh 'tests/e2e/bin/*' .github/workflows/e2e.yml) + tests/e2e/lib.sh 'tests/e2e/bin/*' .github/workflows/e2e.yml '.github/actions/*') # Run on every pull request that runs any case: the bash steps and the # Nextflow pipeline from the fixture's FASTQ, and the parity check of the two. @@ -87,7 +88,7 @@ tracked() { } check() { - local bad="" all base c globs g f hit gl + local bad="" all base c globs g f hit gl out [ -f "$PATHS" ] || die "no table ${PATHS}" all=$(cases); base=$(base_cases) for c in "${SMOKE[@]}"; do @@ -104,6 +105,9 @@ check() { } END { for (k in known) if (!(k in seen)) print "case " k " has no row: add one with the globs of the files it runs" }' | sort) local -A have=() + local scope_globs=() p + for p in "${SCOPE_DIRS[@]}"; do scope_globs+=("${p}/*"); done + scope_globs+=("${SCOPE_FILES[@]}") mapfile -t files < <(tracked) for f in "${files[@]}"; do have[$f]=1; done while IFS=$'\t' read -r c globs; do @@ -120,11 +124,23 @@ check() { done fi [ "$hit" = 1 ] || bad+=$'\n'"case ${c}: ${g} matches no tracked file" + # A file outside the e2e scope never starts a run, so a glob that + # reaches one would let a change to it merge untested. + out="" + if [[ $g != *[*?[]* ]]; then + in_list "$g" "${scope_globs[@]}" || out=$g + else + for f in "${files[@]}"; do + # shellcheck disable=SC2053 # g is a glob on purpose + if [[ $f == $g ]] && ! in_list "$f" "${scope_globs[@]}"; then out=$f; break; fi + done + fi + [ -z "$out" ] || bad+=$'\n'"case ${c}: ${g} matches ${out}, which is outside the e2e scope (SCOPE_DIRS, SCOPE_FILES): a change to it would run no case" done done < <(rows_of "$PATHS" | awk -F'\t' 'NF == 2') bad=$(sed '/^$/d' <<< "$bad") [ -z "$bad" ] || die "$(printf '%s\n' "in ${PATHS}:" "$bad")" - echo "Every one of the $(wc -l <<< "$all" | tr -d ' ') cases has one paths row, and every glob matches a file." + echo "Every one of the $(wc -l <<< "$all" | tr -d ' ') cases has one paths row, and every glob matches a file, inside the e2e scope." } scope() { { printf '%s/**\n' "${SCOPE_DIRS[@]}"; printf '%s\n' "${SCOPE_FILES[@]}"; } | sort; } diff --git a/scripts/run-all.sh b/scripts/run-all.sh index ec353a40..de0149f4 100755 --- a/scripts/run-all.sh +++ b/scripts/run-all.sh @@ -63,6 +63,17 @@ export GENOME_DIR=${GENOME_DIR:?Set GENOME_DIR to your data directory} validate_sample "$SAMPLE" G=$GENOME_DIR S="${GENOME_DIR}/${SAMPLE}" LOG_DIR="${GENOME_DIR}/${SAMPLE}/logs" T=${TOOLS:-} C=${EXTRA_CALLERS:-} on() { [[ "${!1:-}" =~ ^(true|1)$ ]]; } +# The ClinVar files are shared by every sample of GENOME_DIR. bin/owngenome_recheck.py --refresh-clinvar holds this +# lock while it replaces them; the run holds it shared on fd 7 until it ends, so a run waits for a refresh going on +# and a refresh refuses while a run goes on. Without flock (macOS: brew install flock) nothing waits here. +CVLOCK="${G}/clinvar/.refresh.lock" +if [ -d "${G}/clinvar" ] && command -v flock >/dev/null 2>&1; then + [ -e "$CVLOCK" ] || : 2>/dev/null >> "$CVLOCK" || true + if [ -r "$CVLOCK" ]; then + exec 7<"$CVLOCK" + flock -s -n 7 || { echo "Waiting for the ClinVar refresh that holds ${CVLOCK} to end ..."; flock -s 7; } + fi +fi java_major=$(java -version 2>&1 | awk -F'"' '/version "/ {split($2, v, "."); print (v[1] == 1 ? v[2] : v[1]); exit}') || true if ! command -v nextflow >/dev/null 2>&1 || [[ ! "${java_major:-}" =~ ^[0-9]+$ ]] || [ "$java_major" -lt 17 ]; then printf '%s\n' "ERROR: run-all.sh starts the Nextflow pipeline and needs Java 17 or later (found: ${java_major:-none}) and Nextflow ($(command -v nextflow || echo 'not on PATH'))." \ @@ -296,7 +307,7 @@ else echo " benchmark-variants skipped (only one caller VCF: set EXTRA_CALLERS, export NXF_VER=${NXF_VER:-$NEXTFLOW_VERSION} # The live site (docs/live-site.md). Nextflow's weblog is the nf-weblog plugin, which Nextflow 26.10 loads only when # a config names it, so the block lives in a config of this run's own, passed only when OWNGENOME_LIVE_URL is set: a run -# without it neither downloads the plugin nor warns about an unknown weblog scope. Bump the pin with its release. +# without it neither downloads the plugin nor warns about an unknown weblog scope. Renovate bumps the pin (renovate.json). NF_WEBLOG_PLUGIN=nf-weblog@1.2.0 if [ -n "${OWNGENOME_LIVE_URL:-}" ]; then # The ingest listener speaks plain HTTP and serves /ingest/ at its root, so https or a path would lose every post. diff --git a/scripts/setup.sh b/scripts/setup.sh index f12bcd94..5b65b475 100755 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -1985,19 +1985,28 @@ install_archaic() { tail -n 20 "${wk}/liftover.log" echo "[WARN] Picard could not lift chr${c}"; return 1 fi - # The rows that land on the same chromosome and inside the 1000 Genomes - # GRCh38 strict mask, once per position. The mask is sliced to this - # chromosome first: bcftools -T reads the whole targets file into memory, - # and the mask has 4.1 million lines for the genome. + # The rows that land on the same chromosome, once per position. A row + # outside the 1000 Genomes GRCh38 strict mask is kept with its archaic + # calls marked masked (as the FilterBed ones are), not dropped: the + # published workflow hands IBDmix that mask, and IBDmix still scores a + # mismatch at a masked site. + # shellcheck disable=SC2016 # $1 and $2 belong to the inner bash + if ! run_in --rw "$dir" "$BCFTOOLS_IMAGE" bash -euo pipefail -c \ + 'bcftools query -i "CHROM==\"chr$2\"" -f "%CHROM\t%POS\t%REF\t%ALT\t%INFO/Altai\t%INFO/Vindija\t%INFO/Denisova\t%INFO/AF_ALL\t%INFO/AF_AFR\t%INFO/AF_AMR\t%INFO/AF_EAS\t%INFO/AF_EUR\t%INFO/AF_SAS\n" "$1" > "$3/lifted.tsv"' \ + _ "$(cpath "${wk}/hg38.vcf")" "$c" "$(cpath "$wk")"; then + echo "[WARN] Could not read the lifted chr${c}"; return 1 + fi awk -v c="chr${c}" '$1 == c' "${dir}/strict_GRCh38.bed" > "${wk}/strict.bed" if [ ! -s "${wk}/strict.bed" ]; then echo "[WARN] the GRCh38 strict mask has no interval on chr${c}: ${dir}/strict_GRCh38.bed"; return 1 fi - # shellcheck disable=SC2016 # $1 to $3 belong to the inner bash + if ! run_in --rw "$dir" -v "${PGP_ROOT}/bin:/pgp-bin:ro" "$PYTHON_IMAGE" \ + python3 /pgp-bin/archaic_report.py strict --chrom "$c" --strict "$(cpath "${wk}/strict.bed")" \ + --rows "$(cpath "${wk}/lifted.tsv")" --out "$(cpath "${wk}/marked.tsv")"; then + echo "[WARN] Could not mark the strict mask on chr${c}"; return 1 + fi if ! { printf 'chrom\tpos\tref\talt\tAltai\tVindija\tDenisova\tAF_ALL\tAF_AFR\tAF_AMR\tAF_EAS\tAF_EUR\tAF_SAS\n' - run_in --rw "$dir" "$BCFTOOLS_IMAGE" bash -euo pipefail -c \ - 'bcftools query -T "$3" -i "CHROM==\"chr$2\"" -f "%CHROM\t%POS\t%REF\t%ALT\t%INFO/Altai\t%INFO/Vindija\t%INFO/Denisova\t%INFO/AF_ALL\t%INFO/AF_AFR\t%INFO/AF_AMR\t%INFO/AF_EAS\t%INFO/AF_EUR\t%INFO/AF_SAS\n" "$1"' \ - _ "$(cpath "${wk}/hg38.vcf")" "$c" "$(cpath "${wk}/strict.bed")" | awk -F'\t' '!seen[$2]++' + cat "${wk}/marked.tsv" } | gzip -c > "${dir}/chr${c}.tsv.gz.part"; then echo "[WARN] Could not write chr${c} of the table"; return 1 fi diff --git a/tests/e2e/README.md b/tests/e2e/README.md index 774f73a2..b7e09bdb 100644 --- a/tests/e2e/README.md +++ b/tests/e2e/README.md @@ -32,6 +32,8 @@ parallel jobs by `scripts/ci/e2e-shards.sh`, and each case runs at most once. - `lib.sh` holds the shared checks (`check`, `check_eq`, `check_ge`, `has`, `lacks`, `vcf_count`, ...). Every check is a count, a column or a content match, never an exit code alone. +- A case gets 3600 seconds (`CASE_TIMEOUT` of `scripts/ci/e2e-run.sh`). One + that needs longer says so in a line `# case-timeout: SECONDS` of its file. - `bin/docker` lowers `--cpus` to what the runner has; everything else passes through to the real Docker. diff --git a/tests/e2e/archaic-ancestry-1.sh b/tests/e2e/archaic-ancestry-1.sh index 2b616f0c..4baa9348 100755 --- a/tests/e2e/archaic-ancestry-1.sh +++ b/tests/e2e/archaic-ancestry-1.sh @@ -4,15 +4,17 @@ # - setup.sh --archaic builds the table for GRCh37 chr22 only # (ARCHAIC_REGIONS): the three archaic chr22 files from MPI-EVA, checked # against their pinned md5 lists, and the 1000 Genomes sites, lifted to -# GRCh38, every row inside the 1000 Genomes GRCh38 strict mask, and each -# archaic call outside its genome's FilterBed, in a segmental duplication -# or within 5 bp of a 1000 Genomes indel marked masked (code 10 to 12), -# which step 43 hands IBDmix as its mask; +# GRCh38; every archaic call outside its genome's FilterBed, outside the +# 1000 Genomes GRCh38 strict mask, in a segmental duplication or within +# 5 bp of a 1000 Genomes indel is marked masked (code 10 to 12), which +# step 43 hands IBDmix as its mask; # - HG002 (GIAB's v4.2.1 genotypes on chr22): the step runs, the shares are # numbers, the analysed span is most of the chromosome, and the -# Neanderthal total sits inside the range Chen et al.'s own released -# calls give Europeans on chr22 (the check that found the unmasked table -# calling about three times too much); +# Neanderthal total lies between the lowest value (51 kb) of the range +# Chen et al.'s own released calls give Europeans on chr22 and a 1,200 kb +# regression guard (the guard that catches the unmasked table calling +# about three times too much; the whole-genome check decides whether the +# share is in range); # - PLANT, HG002 with one chromosome over chr22:28-29 Mb (GRCh38) copied # from Vindija 33.19 (0/1 wherever Vindija carries the ALT, no record # where it does not): a Neanderthal segment covers 70% of that window or @@ -75,20 +77,26 @@ check "every row is on GRCh38 chr22" \ bash -c 'gzip -dc "$1" | awk -F"\t" "NR > 1 && \$1 != \"chr22\" {bad = 1} END {exit bad}"' _ "$TABLE" check_ge "rows where Vindija carries the ALT inside the plant window" \ "$(gzip -dc "$TABLE" | awk -F'\t' -v f="$W_FROM" -v t="$W_TO" 'NR > 1 && $2 >= f && $2 < t && ($6 == 1 || $6 == 2) {n++} END {print n + 0}')" 100 -# The masks: every FilterBed of the chromosome is installed, and no table row -# sits outside the 1000 Genomes GRCh38 strict mask. Without the masks the -# table carries the sites the method cannot read, and this check goes red. +# The masks: every FilterBed of the chromosome is installed, and every table +# row outside the 1000 Genomes GRCh38 strict mask is kept with its archaic +# calls marked masked (10 to 12), none left as 0, 1 or 2. A table that drops +# those rows has no row outside the mask; one that keeps them unmarked has +# calls 0 to 2 there. Either way this goes red. for f in Altai Vindija33.19 Denisova; do check "the FilterBed of ${f} chr22 is installed" test -s "${DATA}/mask/${f}.chr22.bed.gz" done check "the GRCh38 strict mask is installed" test -s "${DATA}/strict_GRCh38.bed" -OUTSIDE=$(gzip -dc "$TABLE" | awk -F'\t' -v m="${DATA}/strict_GRCh38.bed" ' +gzip -dc "$TABLE" | awk -F'\t' -v m="${DATA}/strict_GRCh38.bed" ' BEGIN { while ((getline l < m) > 0) { split(l, c, "\t"); if (c[1] == "chr22") { s[++n] = c[2]; e[n] = c[3] } } } NR > 1 { p = $2 - 1 while (i < n && e[i + 1] <= p) i++ - if (!(i < n && s[i + 1] <= p && p < e[i + 1])) bad++ } - END { print bad + 0 }') -check_eq "no table row is outside the GRCh38 strict mask" "$OUTSIDE" 0 + if (!(i < n && s[i + 1] <= p && p < e[i + 1])) { + out++ + for (k = 5; k <= 7; k++) if ($k ~ /^(0|1|2)$/) open++ } } + END { print out + 0, open + 0 }' > "${T}/strict-counts.txt" +read -r OUT_ROWS OUT_OPEN < "${T}/strict-counts.txt" +check_ge "table rows outside the GRCh38 strict mask are kept" "$OUT_ROWS" 7000 +check_eq "none of them has an unmarked archaic call" "$OUT_OPEN" 0 # The excluded regions: segmental duplications and the 5 bp around a 1000 # Genomes indel. GRCh37 coordinates, so the test lifts nothing: it checks the # GRCh37 positions of the rejected sites are gone from the GRCh37 input, by @@ -115,15 +123,15 @@ check_ge "analysed span (bp) covers most of chr22" "$(kv "$H" analysed_bp)" 3000 # on the 1000 Genomes phase 3 genomes (LOD 4, union of the Altai and Vindija # segments of 50 kb or more). Over the 497 people of CEU, TSI, GBR, FIN and # IBS, chr22 gives 51 kb at the lowest, 153 kb at the 5th percentile, 541 kb -# at the median, 1,132 kb at the 95th and 1,716 kb at the highest. The gate is -# that spread: a run of this step before the masks landed gave HG002 1,748 kb -# on chr22, more than every one of those 497 people. Where HG002 sits inside -# the spread is printed, not gated: one person's chromosome is one draw, and -# the published calls come from imputed low-coverage genotypes while HG002 is -# a high-coverage benchmark. +# at the median, 1,132 kb at the 95th and 1,716 kb at the highest. One +# chromosome of one person varies widely, so the gate is a regression guard, +# not a claim about the range: 51 kb to 1,200 kb. Runs before the strict-mask +# rows were kept marked gave 1,748 kb (no masks) and 1,177 to 1,190 kb (masked +# codes, rows dropped); the marked table gave 1,138 kb. The published figures +# are printed beside the measured value, and the whole-genome check decides. NEA=$(kv "$H" neanderthal_bp) check_ge "HG002's chr22 Neanderthal total is not under the published lowest (51 kb)" "${NEA:-0}" 51000 -check "and not over the published highest (1,716 kb)" test "${NEA:-0}" -le 1716000 +check "and not over the chr22 regression guard (1,200 kb)" test "${NEA:-0}" -le 1200000 echo " chr22 Neanderthal bp: ${NEA} (published EUR: lowest 51000, p5 153000, median 541000, p95 1132000, highest 1716000)" printf '#### Step 43 HG002 chr22 Neanderthal\n\n%s bp. Published EUR calls of Chen et al. on chr22: lowest 51000, median 541000, 95th 1132000, highest 1716000 bp.\n\n' \ "${NEA:-0}" >> "$E2E_NOTES" diff --git a/tests/e2e/callability-1-catalog.sh b/tests/e2e/callability-1-catalog.sh index f28fd9d3..1f486180 100755 --- a/tests/e2e/callability-1-catalog.sh +++ b/tests/e2e/callability-1-catalog.sh @@ -3,15 +3,20 @@ # from its pinned downloads, byte for byte: a hand edit, or a change to # assets/paralog_loci.tsv, assets/pgx_genes.tsv or bin/cn_loci.tsv that moves # an owner without a rebuild, fails here. Control: a copy of the committed -# file with one owner cell changed must differ from the rebuild. +# file with one owner cell changed must differ from the rebuild. The +# downloads come from the e2e data cache when the runner filled it (entry +# catalogs, mirror first; scripts/ci/e2e-run.sh), else from their sources. . "$(dirname "$0")/lib.sh" OUT="${CASE_TMP}/callability_genes.tsv" COMMITTED="${REPO}/assets/callability_genes.tsv" rm -f "$OUT" echo "+ scripts/ci/build-callability.sh" -CALLABILITY_CACHE="${CASE_TMP}/downloads" bash "${REPO}/scripts/ci/build-callability.sh" "$OUT" 2>&1 | tee "$STEP_LOG" +CACHE="${CASE_TMP}/downloads" +[ -z "${E2E_DATA_CACHE:-}" ] || CACHE="${E2E_DATA_CACHE}/catalogs/callability" +CALLABILITY_CACHE="$CACHE" bash "${REPO}/scripts/ci/build-callability.sh" "$OUT" 2>&1 | tee "$STEP_LOG" check_eq "build-callability.sh exits 0" "${PIPESTATUS[0]}" 0 +echo "${CASE_NAME}: $(grep -c '^Downloading ' "$STEP_LOG") of 7 downloads not found in ${CACHE}" >> "$E2E_NOTES" # same A B: the two files hold the same bytes same() { cmp -s "$1" "$2"; } diff --git a/tests/e2e/ci-e2e-2-locityper.sh b/tests/e2e/ci-e2e-2-locityper.sh new file mode 100755 index 00000000..41722bca --- /dev/null +++ b/tests/e2e/ci-e2e-2-locityper.sh @@ -0,0 +1,113 @@ +#!/usr/bin/env bash +# Step 53 (Locityper) and the database it reads, the shipped scripts as a +# user runs them: setup.sh --locityper-data builds the database, then +# scripts/53-locityper.sh genotypes HG002 and writes the table +# bin/locityper_report.py makes. The image's smoke test repeats their +# commands; this case runs the scripts. +# +# The fixture holds reads over the MHC but none over FCGR or KIR, so the +# database is built for the HLA-DQB1 row of assets/locityper_loci.bed +# (LOCITYPER_LOCI). HLA-DRB1, the other MHC row, is left out for time: with +# both, jellyfish took 35 minutes and the pairwise alignments of +# `locityper augment` 23 more (fork run 42301973505). With HLA-DQB1 alone +# augment takes about a minute, but jellyfish, whose time does not depend on +# the loci, still takes 36 (38 minutes in all, fork run 42325306077), so the +# case keeps a limit of 90 minutes instead of the runner's 60 (the +# case-timeout line below). It runs in a GENOME_DIR of its own with a +# reference of chr6 and chr20 cut from the fixture's, so jellyfish +# counts 235 Mb instead of the fixture's 1.8 Gb; the fixture's reads over the +# MHC and chr20:10.0-10.5 Mb are aligned to it, since Locityper needs the +# BAM's contigs to be the reference's. Step 53 measures depth on chr20 +# (LOCITYPER_BG_REGION): the fixture has no chr17, where Locityper's own +# background region is. +# +# Red controls: the step stops without a database, without the BAM's index, +# and with a database whose DONE names another reference. +# case-timeout: 5400 +. "$(dirname "$0")/lib.sh" + +G="${CASE_TMP}/genome" +S="${G}/${SAMPLE}" +rm -rf "$G" +mkdir -p "${G}/reference" "${S}/aligned" +U="$(id -u):$(id -g)" + +# --- the reference and the reads ------------------------------------------------- +echo "+ a chr6 and chr20 reference cut from the fixture's, and the reads over the MHC and chr20 aligned to it" +in_genome "$SAMTOOLS_IMAGE" samtools faidx reference/GRCh38_no_alt_analysis_set.fasta chr6 chr20 \ + > "${G}/reference/GRCh38_no_alt_analysis_set.fasta" +docker run --rm --user "$U" -v "${G}:/g" "$SAMTOOLS_IMAGE" samtools faidx /g/reference/GRCh38_no_alt_analysis_set.fasta +check_eq "the cut reference has chr6 and chr20" "$(cut -f1 "${G}/reference/GRCh38_no_alt_analysis_set.fasta.fai" 2>/dev/null | paste -sd' ' -)" "chr6 chr20" +docker run --rm --user "$U" -v "${FIXTURE_DIR}:/in:ro" -v "${G}:/g" "$SAMTOOLS_IMAGE" bash -c ' + samtools view -u -F 0x900 /in/HG002_slice.bam chr6:31500001-33100000 chr20:10000001-10500000 \ + | samtools collate -O -u - /g/collate | samtools fastq -n -1 /g/r1.fq.gz -2 /g/r2.fq.gz -0 /dev/null -s /dev/null -' +docker run --rm --user "$U" -v "${G}:/g" "$MINIMAP2_IMAGE" \ + minimap2 -t "${THREADS:-4}" -a -x sr -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}" /g/reference/GRCh38_no_alt_analysis_set.fasta /g/r1.fq.gz /g/r2.fq.gz \ + | docker run --rm -i --user "$U" -v "${G}:/g" "$SAMTOOLS_IMAGE" bash -c \ + "samtools sort -o /g/${SAMPLE}/aligned/${SAMPLE}_sorted.bam - && samtools index /g/${SAMPLE}/aligned/${SAMPLE}_sorted.bam" +rm -f "${G}/r1.fq.gz" "${G}/r2.fq.gz" +check_ge "reads aligned to the cut reference" \ + "$(docker run --rm -v "${G}:/g:ro" "$SAMTOOLS_IMAGE" samtools view -c -F 0x904 "/g/${SAMPLE}/aligned/${SAMPLE}_sorted.bam" 2>/dev/null)" 10000 + +LOCI="${CASE_TMP}/loci.bed" +awk -F'\t' '!/^#/ && $4 == "HLA-DQB1"' "${REPO}/assets/locityper_loci.bed" > "$LOCI" +check_eq "the HLA-DQB1 row of assets/locityper_loci.bed" "$(cut -f1,4 "$LOCI" | tr '\t' ' ')" "chr6 HLA-DQB1" + +# run53 LOG: step 53 on HG002 in G, its exit code in STEP_RC. +run53() { + echo "+ scripts/53-locityper.sh ${SAMPLE} (GENOME_DIR=${G})" + GENOME_DIR="$G" LOCITYPER_BG_REGION=chr20:10000001-10500000 THREADS="${THREADS:-4}" \ + "${REPO}/scripts/53-locityper.sh" "$SAMPLE" > "$1" 2>&1 + STEP_RC=$? + cat "$1" +} + +# --- red control: no database yet ------------------------------------------------ +run53 "${CASE_TMP}/nodb.log" +check_eq "red control: step 53 stops (exit 1) before the database is built" "$STEP_RC" 1 +check "red control: and says how to build it" grep -q "Locityper's database not found: .*setup.sh --locityper-data" "${CASE_TMP}/nodb.log" + +# --- setup.sh --locityper-data --------------------------------------------------- +DB="${G}/reference/locityper-${LOCITYPER_PANEL_VERSION}" +echo "+ LOCITYPER_LOCI=${LOCI} scripts/setup.sh --locityper-data ${G}" +GENOME_DIR="$G" LOCITYPER_LOCI="$LOCI" THREADS="${THREADS:-4}" "${REPO}/scripts/setup.sh" --locityper-data "$G" 2>&1 | tee "$STEP_LOG" +check_eq "setup.sh --locityper-data exits 0" "${PIPESTATUS[0]}" 0 +check "and says the database is built" grep -q "\[OK\] Locityper's database (step 53, 1 loci" "$STEP_LOG" +check "the database holds HLA-DQB1 (Locityper's success file)" test -e "${DB}/db/loci/HLA-DQB1/success" +check_eq "its loci.bed is the HLA-DQB1 row" "$(cut -f4 "${DB}/loci.bed" 2>/dev/null | paste -sd' ' -)" "HLA-DQB1" +check "its DONE names the reference's .fai" grep -q "fai=$(sha256sum "${G}/reference/GRCh38_no_alt_analysis_set.fasta.fai" | cut -d' ' -f1) " "${DB}/DONE" +GENOME_DIR="$G" LOCITYPER_LOCI="$LOCI" THREADS="${THREADS:-4}" "${REPO}/scripts/setup.sh" --locityper-data "$G" > "${CASE_TMP}/again.log" 2>&1 +check_eq "a second setup.sh --locityper-data exits 0" "$?" 0 +check "and builds nothing again" grep -q "already built" "${CASE_TMP}/again.log" + +# --- step 53 --------------------------------------------------------------------- +run53 "${CASE_TMP}/step.log" +check_eq "scripts/53-locityper.sh exits 0" "$STEP_RC" 0 +TABLE="${S}/locityper/${SAMPLE}_locityper.tsv" +check_eq "the table's columns are those of bin/locityper_report.py" "$(head -n 1 "$TABLE" 2>/dev/null | tr '\t' ' ')" \ + "locus region status genotype haplotype_1 haplotype_2 quality total_reads unexplained_reads weighted_distance warnings filter reason" +check_eq "one row per locus of the database" "$(awk -F'\t' 'NR > 1 {print $1}' "$TABLE" 2>/dev/null | paste -sd' ' -)" "HLA-DQB1" +l=HLA-DQB1 +row=$(awk -F'\t' -v l="$l" '$1 == l' "$TABLE" 2>/dev/null) +check_eq "${l}: called" "$(cut -f3 <<< "$row")" called +check "${l}: two panel haplotypes" has '^[^,.][^,]*,[^,.][^,]*$' "$(cut -f4 <<< "$row")" +check "${l}: a numeric quality" has '^[0-9]+(\.[0-9]+)?$' "$(cut -f7 <<< "$row")" +check_ge "${l}: reads used" "$(cut -f8 <<< "$row")" 1 +check "${l}: the authors' filter is pass or fail" has '^(pass|fail)$' "$(cut -f12 <<< "$row")" +check "${l}: Locityper's result file" test -s "${S}/locityper/work/gt/loci/${l}/res.json.gz" +check "Locityper's log beside the table" test -s "${S}/locityper/${SAMPLE}_locityper.log" +echo "${CASE_NAME}: $(awk -F'\t' 'NR > 1 {printf "%s %s q%s %s; ", $1, $4, $7, $12}' "$TABLE" 2>/dev/null)" >> "$E2E_NOTES" + +# --- red controls: no BAM index, a database of another reference ------------------- +mv "${S}/aligned/${SAMPLE}_sorted.bam.bai" "${CASE_TMP}/bai" +run53 "${CASE_TMP}/nobai.log" +check_eq "red control: step 53 stops (exit 1) without the BAM's index" "$STEP_RC" 1 +check "red control: and names it" grep -q "File not found: .*${SAMPLE}_sorted.bam.bai" "${CASE_TMP}/nobai.log" +mv "${CASE_TMP}/bai" "${S}/aligned/${SAMPLE}_sorted.bam.bai" +sed -i 's/fai=[0-9a-f]* /fai=0000 /' "${DB}/DONE" +run53 "${CASE_TMP}/stale.log" +check_eq "red control: step 53 stops (exit 1) on a database built from another reference" "$STEP_RC" 1 +check "red control: and says so" grep -q "was built from another reference" "${CASE_TMP}/stale.log" + +rm -rf "${G:?}" +finish diff --git a/tests/e2e/cn-intervals-1-catalog.sh b/tests/e2e/cn-intervals-1-catalog.sh index 21b741f7..5b0db18c 100755 --- a/tests/e2e/cn-intervals-1-catalog.sh +++ b/tests/e2e/cn-intervals-1-catalog.sh @@ -3,10 +3,13 @@ # writes from its pinned downloads (ClinGen's recurrent CNV BED V2.1 and three # GIAB v3.6 stratifications, each checked against its sha256): a hand edit of # an interval row fails here. Control: a copy of the catalog with one target -# moved by 1 kb must fail the same check. +# moved by 1 kb must fail the same check. The downloads come from the e2e +# data cache when the runner filled it (entry catalogs, mirror first; +# scripts/ci/e2e-run.sh), else from their sources. . "$(dirname "$0")/lib.sh" CACHE="${CASE_TMP}/downloads" +[ -z "${E2E_DATA_CACHE:-}" ] || CACHE="${E2E_DATA_CACHE}/catalogs/cn-intervals" mkdir -p "$CACHE" echo "+ scripts/ci/build-cn-intervals.py --check" python3 "${REPO}/scripts/ci/build-cn-intervals.py" --cache "$CACHE" --check 2>&1 | tee "$STEP_LOG" diff --git a/tests/e2e/lib.sh b/tests/e2e/lib.sh index 6d5de4ec..e469c287 100755 --- a/tests/e2e/lib.sh +++ b/tests/e2e/lib.sh @@ -7,6 +7,8 @@ # A case never stops at the first failed check: it runs them all, prints one # [PASS] or [FAIL] line each, and `finish` exits 1 if any failed. Every check is # a count, a column or a content match; an exit code alone is never enough. +# A case that ran no check at all fails too, unless it set NO_CHECKS_REASON +# before `finish` (a case that has nothing to do on this run says why). set -uo pipefail : "${REPO:?run the cases through scripts/ci/e2e-run.sh}" @@ -18,11 +20,13 @@ CASE_NAME="$(basename "$0" .sh)" CASE_TMP="${E2E_WORK}/tmp/${CASE_NAME}" mkdir -p "$CASE_TMP" FAILS=0 +CHECKS=0 +NO_CHECKS_REASON="" STEP_RC="" STEP_LOG="${CASE_TMP}/step.log" -pass() { echo "[PASS] $*"; } -fail() { echo "[FAIL] $*"; FAILS=$((FAILS + 1)); } +pass() { echo "[PASS] $*"; CHECKS=$((CHECKS + 1)); } +fail() { echo "[FAIL] $*"; CHECKS=$((CHECKS + 1)); FAILS=$((FAILS + 1)); } # check : passes when the command succeeds. check() { @@ -180,11 +184,22 @@ setup_ancestry_panel() { echo "${CASE_NAME}: ancestry panel from ${from:-EBI}" >> "$E2E_NOTES" } +# finish: exit 1 when a check failed, or when no check ran and the case gave +# no NO_CHECKS_REASON; scripts/ci/e2e-run.sh reads the last line it prints. finish() { if [ "$FAILS" -gt 0 ]; then echo "${CASE_NAME}: ${FAILS} check(s) failed" exit 1 fi - echo "${CASE_NAME}: all checks passed" + if [ "$CHECKS" -eq 0 ]; then + if [ -n "$NO_CHECKS_REASON" ]; then + echo "${CASE_NAME}: no checks run: ${NO_CHECKS_REASON}" + echo "${CASE_NAME}: no checks run: ${NO_CHECKS_REASON}" >> "$E2E_NOTES" + exit 0 + fi + echo "${CASE_NAME}: no check ran, so the case proves nothing (set NO_CHECKS_REASON when that is meant)" + exit 1 + fi + echo "${CASE_NAME}: all ${CHECKS} checks passed" exit 0 } diff --git a/tests/e2e/prs-3-panel-measure.sh b/tests/e2e/prs-3-panel-measure.sh index ce43535b..8e6e76a1 100644 --- a/tests/e2e/prs-3-panel-measure.sh +++ b/tests/e2e/prs-3-panel-measure.sh @@ -65,7 +65,7 @@ if [ "${1:-}" = --download ]; then fi if ! measure_on; then - echo "the 1000 Genomes panel measurement runs on an E2E run dispatched on a branch other than main, or with PGSC_MEASURE=1." + NO_CHECKS_REASON="the 1000 Genomes panel measurement runs on an E2E run dispatched on a branch other than main, or with PGSC_MEASURE=1" finish fi command -v nextflow >/dev/null || { fail "nextflow is not on PATH"; finish; } diff --git a/tests/test_annotsv_postfilter.py b/tests/test_annotsv_postfilter.py index 0e84c9b2..4335f122 100755 --- a/tests/test_annotsv_postfilter.py +++ b/tests/test_annotsv_postfilter.py @@ -32,7 +32,10 @@ Manta deletions at one place, DHFFC 0.95 and 0.5: the first is class 3 (no depth support), the second stays 4; a duplication with DHBFC 1.1 is class 3; Delly RDCN 2 on a deletion is no support, RDCN 1 is; on - chrX a duplication with RDCN 2 counts; a record with no depth value is + chrX a duplication with RDCN 2 counts with no sex given; with --sex + and --par-bed a male's chrX deletion needs RDCN 0, a female's chrX + duplication 3 or more, and inside the PAR a duplication needs 3 for + any sex; a record with no depth value is "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 @@ -157,6 +160,37 @@ def run7(work, name, with_rules=True): return {r["ID"]: r for r in csv.DictReader(f, delimiter="\t")}, p.stderr +SEX_ROWS = "".join([ + row7("XDEL1", "X", 50000000, 50050000, "DEL", 4, "delly", "G20", "GT:RDCN", "1/1:1"), + row7("XDEL0", "X", 51000000, 51050000, "DEL", 4, "delly", "G21", "GT:RDCN", "1/1:0"), + row7("XDUP2", "X", 52000000, 52050000, "DUP", 4, "delly", "G22", "GT:RDCN", "0/1:2"), + row7("PARDUP2", "X", 1000000, 1050000, "DUP", 4, "delly", "G23", "GT:RDCN", "0/1:2"), + row7("PARDEL1", "X", 1000000, 1050000, "DEL", 4, "delly", "G24", "GT:RDCN", "0/1:1"), + row7("YDEL1", "Y", 20000000, 20050000, "DEL", 4, "delly", "G25", "GT:RDCN", "1/1:1"), + row7("YDUP2", "Y", 21000000, 21050000, "DUP", 4, "delly", "G26", "GT:RDCN", "0/1:2"), + # PAR1 ends at chrX:2,781,479 (assets/par_grch38.bed): the last base in + # it, then the first base past it. + row7("PAREDGE_IN", "X", 2781479, 2831479, "DUP", 4, "delly", "G27", "GT:RDCN", "0/1:2"), + row7("PAREDGE_OUT", "X", 2781480, 2831480, "DUP", 4, "delly", "G28", "GT:RDCN", "0/1:2"), + row7("BIGDUP", "20", 48000000, 48050000, "DUP", 4, "delly", "G29", "GT:RDCN", "0/1:1000000"), +]) + + +def run_sex(work, name, sex, par): + """{ID: written row} of the command line on SEX_ROWS with --sex and --par-bed.""" + c45 = os.path.join(work, name + ".tsv") + with open(c45, "w") as f: + f.write("\t".join(HEAD7) + "\n" + SEX_ROWS) + dst = os.path.join(work, name + "_checked.tsv") + p = subprocess.run([sys.executable, os.path.join(REPO, "bin", "annotsv_postfilter.py"), "--class45", c45, + "--sex", sex, "--par-bed", par, "--out", dst], capture_output=True, text=True) + if p.returncode != 0 or not os.path.isfile(dst): + print(f" (the command failed: {p.stderr.strip()[-300:]})") + return {}, p.stderr + with open(dst) as f: + return {r["ID"]: r for r in csv.DictReader(f, delimiter="\t")}, p.stderr + + def rules_section(work): got, err = run7(work, "rules") @@ -193,7 +227,7 @@ def col(i, c): check("Delly RDCN 2 on a deletion: no support; RDCN 1: support", (col("DELLY2", "Depth_support"), col("DELLY1", "Depth_support")) == ("no", "yes"), (got.get("DELLY2"), got.get("DELLY1"))) - check("a chrX duplication with RDCN 2 counts (the sex is not known here)", col("XDUP", "Depth_support") == "yes", + check("a chrX duplication with RDCN 2 counts when the sex is not given", col("XDUP", "Depth_support") == "yes", got.get("XDUP")) check("no depth value: not checked, class kept", (col("NODEPTH", "Depth_support"), col("NODEPTH", "ACMG_class")) == ("not checked", "4"), got.get("NODEPTH")) @@ -203,6 +237,25 @@ def col(i, c): check("the same call of 800 bp had no RDCN gate: not checked, class kept", (col("DELLY_SHORT", "Depth_support"), col("DELLY_SHORT", "ACMG_class")) == ("not checked", "4"), got.get("DELLY_SHORT")) + # Delly RDCN on chrX and chrY by sex and PAR (--sex, --par-bed), the rule + # cnvpytor_to_vcf.py applies when step 05 adds a Delly call. + par = os.path.join(REPO, "assets", "par_grch38.bed") + yes, no = ("yes", "4"), ("no", "3") + for sex, want_sex in (("male", {"XDEL1": no, "XDUP2": yes, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": no, "YDUP2": yes, "PAREDGE_IN": no, "PAREDGE_OUT": yes, + "BIGDUP": yes}), + ("female", {"XDEL1": yes, "XDUP2": no, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": yes, "YDUP2": no, "PAREDGE_IN": no, "PAREDGE_OUT": no, + "BIGDUP": yes}), + ("", {"XDEL1": yes, "XDUP2": yes, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": yes, "YDUP2": yes, "PAREDGE_IN": no, "PAREDGE_OUT": yes, + "BIGDUP": yes})): + sx, _ = run_sex(work, "sex_" + (sex or "none"), sex, par) + got_sex = {i: (sx.get(i, {}).get("Depth_support"), sx.get(i, {}).get("ACMG_class")) for i in want_sex} + check(f"sex {sex or 'not given'}: chrX and chrY RDCN against 1 copy (male), 2 (female), either (not " + f"given); inside the PAR, up to its last base, 2 copies; RDCN 1000000 is a gain", + got_sex == want_sex, got_sex) + moved = sorted(i for i, r in got.items() if r["ACMG_class"] != r["AnnotSV_class"]) 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) diff --git a/tests/test_annotsv_sv_inputs.py b/tests/test_annotsv_sv_inputs.py index e014e5ce..a153490d 100755 --- a/tests/test_annotsv_sv_inputs.py +++ b/tests/test_annotsv_sv_inputs.py @@ -226,6 +226,17 @@ def test_svlen(w): rc, out2, err2 = run("sv_svlen_fix.py", again) check("a second run changes nothing", rc == 0 and out2 == out, err2) + # One sample only: the reports read the GT and PS of the first sample + # column, so a second person's column (here the carrier, after a 0/0) + # must stop the step, not be read silently. + two = os.path.join(w, "two.vcf") + with open(two, "w") as f: + f.write(SV.replace("\ts1\n", "\ts1\ts2\n") + + "chr20\t1300000\tsym\tN\t\t50\tPASS\tEND=1310000;SVTYPE=DEL\tGT\t0/0\t0/1\n") + rc, out3, err3 = run("sv_svlen_fix.py", two) + check("a VCF with two sample columns stops (exit 2) and names them", + rc == 2 and "holds 2 samples (s1, s2)" in err3 and not records(out3), (rc, err3)) + BASE = ( "##fileformat=VCFv4.1\n##contig=\n##contig=\n" @@ -295,6 +306,39 @@ def test_extra_callers(w): "delly: 2 PASS DEL/DUP calls, 1 added" in err and "cnvpytor: 3 DEL/DUP calls with e-val < 0.01 and q0 < 0.5, 2 added" in err, err) + # A Manta record that is filtered or genotyped 0/0 is no call: it must not + # hide the Delly PASS deletion it overlaps. + for label, gt, filt in (("a 0/0", "0/0", "PASS"), ("a MinGQ", "0/1", "MinGQ"), ("a ./.", "./.", "PASS")): + b = os.path.join(w, "base_nocall.vcf") + with open(b, "w") as f: + f.write(BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\t{0}\tEND=2049900;SVTYPE=DEL;SVLEN=-49800\tGT\t{1}\n" + .format(filt, gt)) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", b, "--delly", delly) + ids = [r[2] for r in records(out)] + check("{0} Manta deletion stays and does not hide Delly's DEL00000002".format(label), + rc == 0 and "MantaDEL:9" in ids and "DEL00000002" in ids, (ids, err)) + with open(b, "w") as f: + f.write(BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\tPASS\tEND=2049900;SVTYPE=DEL;SVLEN=-49800\tGT\t0|1\n") + rc, out, err = run("cnvpytor_to_vcf.py", "--base", b, "--delly", delly) + check("control: a PASS 0|1 Manta deletion still covers Delly's DEL00000002", + rc == 0 and "DEL00000002" not in [r[2] for r in records(out)], err) + sites = "".join("\t".join(l.split("\t")[:8]) + "\n" for l in (BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\tPASS\tEND=2049900;SVTYPE=DEL\tGT\t0/0\n").splitlines()) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", "-", "--delly", delly, stdin=sites) + check("a sites-only base deletion covers Delly's DEL00000002 (no genotype to read)", + rc == 0 and "DEL00000002" not in [r[2] for r in records(out)], (records(out), err)) + + # A Delly PASS deletion genotyped 0/0 or ./. is no call either: it is not + # added, and CNVpytor's deletion at the same place is added instead. + for gt in ("0/0", "./."): + with gzip.open(delly, "wt") as f: + f.write(DELLY.replace("\t1/1:40\n", "\t{0}:40\n".format(gt))) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", base, "--delly", delly, "--cnvpytor", cnvp) + ids = [r[2] for r in records(out)] + check("a {0} Delly deletion is not added and does not hide CNVpytor_del0".format(gt), + rc == 0 and "DEL00000002" not in ids and "CNVpytor_del0" in ids, (ids, err)) + with gzip.open(delly, "wt") as f: + f.write(DELLY) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", "-", "--base-name", "sv_vcf", stdin=BASE) check("no other caller: stdin in, the base out with SVCALLER=sv_vcf", rc == 0 and "none given" in err and [info(r).get("SVCALLER") for r in records(out)] == ["sv_vcf"] * 2, err) diff --git a/tests/test_archaic_report.py b/tests/test_archaic_report.py index 087409cd..a54d508e 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -9,6 +9,8 @@ site no archaic genome was called at, or only called as reference with no modern variation, is left out; a position given twice by 1000 Genomes is left out; + 1b. strict: a lifted row outside the GRCh38 strict mask is kept once per + position with every archaic call marked masked (9 stays 9), not dropped; 2. genotypes: a masked archaic call (10 to 12) goes to IBDmix as its ALT count and its position into that genome's mask file; the stand-in columns carry the population's ALT count; your column is your ALT @@ -163,6 +165,30 @@ def table_info(out, extra=()): check("a mask with no interval on the chromosome stops the build, instead of emptying the table", r.returncode == 1 and "no interval" in r.stderr, r.stderr) +# --- 1b. strict: rows outside the GRCh38 strict mask are marked, not dropped ------- +# Lifted rows (13 columns, no header). Strict mask 1-based 101-300 on chr22: +# 100 is outside, 150 and 300 inside, 301 outside. A duplicate position keeps +# the first row. Codes 0, 1, 2 outside become 10, 11, 12; 9 stays 9; a code +# already masked is not marked twice. +rows = ["chr22\t100\tA\tG\t2\t0\t9\t0.5\t0\t0\t0\t0\t0", + "chr22\t150\tA\tG\t1\t1\t1\t0.5\t0\t0\t0\t0\t0", + "chr22\t301\tC\tT\t0\t10\t12\t0.5\t0\t0\t0\t0\t0", + "chr22\t301\tC\tT\t2\t2\t2\t0.5\t0\t0\t0\t0\t0", + "chr22\t300\tG\tA\t0\t0\t0\t0.5\t0\t0\t0\t0\t0"] +write(f"{TMP}/lifted.tsv", "\n".join(rows) + "\n") +write(f"{TMP}/strict.bed", "chr22\t100\t300\n") +r = run("strict", "--chrom", "22", "--strict", f"{TMP}/strict.bed", "--rows", f"{TMP}/lifted.tsv", "--out", f"{TMP}/marked.tsv") +check("strict exits 0", r.returncode == 0, r.stderr) +got = [l.rstrip("\n").split("\t") for l in open(f"{TMP}/marked.tsv")] +by = {int(c[1]): c[4:7] for c in got} +check("every position is kept, once", [int(c[1]) for c in got] == [100, 150, 301, 300], [c[1] for c in got]) +check("outside the strict mask 2, 0 and unknown become 12, 10 and 9", by[100] == ["12", "10", "9"], by.get(100)) +check("inside the strict mask the calls are untouched (last base 300 too)", by[150] == ["1", "1", "1"] and by[300] == ["0", "0", "0"], (by.get(150), by.get(300))) +check("a masked call stays as it is, and the first of a duplicate wins", by[301] == ["10", "10", "12"], by.get(301)) +write(f"{TMP}/strict_other.bed", "chr21\t0\t10\n") +r = run("strict", "--chrom", "22", "--strict", f"{TMP}/strict_other.bed", "--rows", f"{TMP}/lifted.tsv", "--out", f"{TMP}/marked2.tsv") +check("a strict mask with no interval on the chromosome stops the step", r.returncode != 0 and "no interval" in r.stderr, r.stderr) + # --- 2. genotypes ------------------------------------------------------------------ tdir = f"{TMP}/table" os.makedirs(tdir) diff --git a/tests/test_e2e_data_cache.sh b/tests/test_e2e_data_cache.sh index 05c47bcd..93348564 100755 --- a/tests/test_e2e_data_cache.sh +++ b/tests/test_e2e_data_cache.sh @@ -1,8 +1,9 @@ #!/usr/bin/env bash -# The data cache entries hla and pgs of scripts/ci/e2e-run.sh: --data-check -# passes a whole entry and fails a changed, missing or partial one (the e2e -# workflow saves and uses only an entry that passes), and --data-key names -# the releases; a data fill that runs past E2E_DATA_BUDGET stops with its +# The data cache entries hla, pgs and catalogs of scripts/ci/e2e-run.sh: +# --data-check passes a whole entry and fails a changed, missing or partial +# one (the e2e workflow saves and uses only an entry that passes), and +# --data-key names the releases; kill_tree stops a parent before its +# children; a data fill that runs past E2E_DATA_BUDGET stops with its # message (a planted slow mirror on 127.0.0.1, the sources blocked); and # .github/actions/free-disk names a cleanup that failed in # its summary (fake sudo, docker and df). No network, no docker. @@ -69,6 +70,17 @@ else echo "FAIL: the hla key does not match its inputs (${other:-?} ${pin:-?} ${ if "$RUN" --data-key pgs | grep -qE '^e2e-data-pgs-[0-9a-f]{16}$'; then echo "PASS: the pgs key is a hash" else echo "FAIL: the pgs key"; failed=1; fi +# catalogs: the downloads of the two table builders, each in its folder. +if "$RUN" --data-key catalogs | grep -qE '^e2e-data-catalogs-[0-9a-f]{16}$'; then echo "PASS: the catalogs key is a hash" +else echo "FAIL: the catalogs key"; failed=1; fi +n=$( { bash "${REPO}/scripts/ci/build-callability.sh" --list; python3 "${REPO}/scripts/ci/build-cn-intervals.py" --list; } | awk 'NF == 3 && $2 ~ /^https:/ && $3 ~ /^[0-9a-f]{64}$/' | wc -l | tr -d ' ') +if [ "$n" = 11 ]; then echo "PASS: the two builders list 11 pinned downloads (7 and 4)" +else echo "FAIL: the builders list ${n} pinned downloads, want 11"; failed=1; fi +mkdir -p "${T}/catalogs/callability" +expect "catalogs: an entry with no file fails" 1 catalogs +bash "${REPO}/scripts/ci/build-callability.sh" --list | head -n 1 | { read -r name _ _; echo planted > "${T}/catalogs/callability/${name}"; } +expect "catalogs: a file that does not match its pin fails" 1 catalogs + # --- the data fill: the mirror, its deadline and the budget ------------------ # A mirror on 127.0.0.1; every other host is sent to a closed port, so the # PGS Catalog, the source after the mirror, always fails. Its modes: @@ -160,6 +172,28 @@ if [ "$rc" != 0 ] && [ -z "$(find "${SL}/bad.cache" -name '*.txt.gz' 2> /dev/nul echo "PASS: a mirror whose bytes do not match the pins is refused and leaves no score in the entry (exit ${rc})" else echo "FAIL: the bad mirror gave exit ${rc}: $(find "${SL}/bad.cache" -type f 2> /dev/null | head -3) $(tail -n 3 "${SL}/bad.log")"; failed=1; fi +# --- kill_tree: the parent stops before its children ------------------------ +# A parent that starts its next command when its child dies (mfetch starts +# the source download when the mirror try ends) must not get to start it: a +# bash whose first sleep is killed would start a second one. Stopped +# children first, the second sleep ran on after kill_tree returned. +fn=$(sed -n '/^kill_tree() {/,/^}/p' "$RUN") +[ -n "$fn" ] || { echo "FAIL: no kill_tree function in e2e-run.sh"; exit 1; } +eval "$fn" +KT="${T}/kill-tree" +mkdir -p "$KT" +bash -c 'sleep 30 & first=$!; sleep 40 & wait "$first"; sleep 31 & echo $! > "$1/second"; wait' _ "$KT" & kp=$! +for _ in $(seq 50); do [ "$(pgrep -P "$kp" 2> /dev/null | wc -l)" -ge 2 ] && break; sleep 0.1; done +kill_tree "$kp" +wait "$kp" 2> /dev/null +sleep 1 +if [ -s "${KT}/second" ] && kill -0 "$(cat "${KT}/second")" 2> /dev/null; then + echo "FAIL: kill_tree let the parent start a new child after its first one was stopped"; failed=1 + kill "$(cat "${KT}/second")" 2> /dev/null +elif pgrep -P "$kp" > /dev/null 2>&1 || kill -0 "$kp" 2> /dev/null; then + echo "FAIL: kill_tree left the parent or its child running"; failed=1 +else echo "PASS: kill_tree stops the parent first, so it starts nothing after its child is stopped"; fi + # --- free-disk: a cleanup that fails is named in the summary ----------------- # The action's script, run as GitHub runs a composite bash step, with sudo, # docker and df faked on PATH: sudo fails for the planted path, docker's prune diff --git a/tests/test_e2e_finish.sh b/tests/test_e2e_finish.sh new file mode 100755 index 00000000..cb32cb94 --- /dev/null +++ b/tests/test_e2e_finish.sh @@ -0,0 +1,31 @@ +#!/usr/bin/env bash +# finish in tests/e2e/lib.sh: an e2e case that ran no check fails, unless it +# set NO_CHECKS_REASON, which finish prints and adds to the notes. A case +# with a passing check passes, one with a failed check fails. Planted cases, +# no docker, no network. +set -uo pipefail +REPO="$(cd "$(dirname "$0")/.." && pwd)" +T=$(mktemp -d) +trap 'rm -rf "$T"' EXIT +failed=0 +mkdir -p "${T}/work" "${T}/genome" "${T}/fixture" "${T}/cases" +# expect DESC EXIT BODY: a case of BODY between lib.sh and finish exits EXIT. +expect() { + local f="${T}/cases/planted-$RANDOM.sh" rc=0 + printf '. "%s/tests/e2e/lib.sh"\n%s\nfinish\n' "$REPO" "$3" > "$f" + : > "${T}/notes.md" + REPO=$REPO GENOME_DIR="${T}/genome" SAMPLE=HG002 FIXTURE_DIR="${T}/fixture" E2E_WORK="${T}/work" \ + E2E_NOTES="${T}/notes.md" bash "$f" > "${T}/out" 2>&1 || rc=$? + if [ "$rc" = "$2" ]; then echo "PASS: $1 (exit ${rc}: $(tail -n 1 "${T}/out"))" + else echo "FAIL: $1 (exit ${rc}, want $2): $(cat "${T}/out")"; failed=1; fi +} +expect "a case that ran no check fails" 1 ':' +if grep -q "no check ran, so the case proves nothing" "${T}/out"; then echo "PASS: and says it ran no check" +else echo "FAIL: no message for a case with no check"; failed=1; fi +expect "a case with no check and a reason passes" 0 'NO_CHECKS_REASON="nothing to measure on this run"' +if grep -q "no checks run: nothing to measure on this run" "${T}/notes.md"; then echo "PASS: and its reason is in the notes" +else echo "FAIL: the reason is not in the notes"; failed=1; fi +expect "a case with one passing check passes" 0 'check_eq "one" 1 1' +expect "a case with a failed check fails" 1 'check_eq "one" 1 2' +expect "a failed check fails the case even with a reason set" 1 'NO_CHECKS_REASON=x; check_eq "one" 1 2' +exit "$failed" diff --git a/tests/test_e2e_scope.sh b/tests/test_e2e_scope.sh index daf4d102..b7717c68 100755 --- a/tests/test_e2e_scope.sh +++ b/tests/test_e2e_scope.sh @@ -3,8 +3,8 @@ # required e2e check trusts it: a case it leaves out never runs. Each check # here is run against a planted fault first: a case with no paths row, a # second row, a row for no case, a base row with globs, a glob that matches -# no file. Then the picks: nothing for a path outside the scope, every case -# for shared wiring, the smoke cases plus the matching ones for a helper, +# no file, a glob that matches a file outside the e2e scope. Then the +# picks: nothing for a path outside the scope, every case for shared wiring, the smoke cases plus the matching ones for a helper, # the case whose file or row changed, the cases those need, never a base # case. Last, the real table must pass the check. set -euo pipefail @@ -57,6 +57,10 @@ broken "a glob that matches no file" "bin/no_such_helper.py matches no tracked f broken "a module glob that matches no file" "modules/local/nope/\\* matches no tracked file" "^lpa-1${TAB}" \ "lpa-1${TAB}modules/local/nope/*" broken "a row of one column" "has 1 columns" "^traits-2${TAB}" "traits-2" +broken "a glob that matches a file outside the e2e scope" "docs/testing.md matches docs/testing.md, which is outside the e2e scope" \ + "^traits-2${TAB}" "traits-2${TAB}bin/traits_report.py docs/testing.md" +broken "a wildcard glob that reaches outside the e2e scope" "docs/\\* matches docs/" "^traits-2${TAB}" \ + "traits-2${TAB}bin/traits_report.py docs/*" mkdir -p "${WORK}/nosmoke" cp "${E2E_CASE_DIR}"/*.sh "${WORK}/nosmoke/" rm "${WORK}/nosmoke/nextflow-from-fastq-2-nextflow.sh" diff --git a/tests/test_owngenome_recheck.py b/tests/test_owngenome_recheck.py index bb5b8a50..bb98fcb7 100644 --- a/tests/test_owngenome_recheck.py +++ b/tests/test_owngenome_recheck.py @@ -27,6 +27,20 @@ 12. check() itself names step 53 rerun-ready when the run used an older LOCITYPER_PANEL_VERSION than this checkout pins, and not when they match; its action names setup.sh --locityper-data, which the main setup leaves out. + 13. an image re-pinned at the same tag with a new digest is rerun-ready, for a pin + with no upstream lookup (step 50's) and for a looked-up one (PharmCAT); the + same reference on both sides is no-newer-release; + 14. a pin with no upstream lookup that the run's manifest lacks is unknown, never + no-newer-release, names no step and does not make the command exit 1; + 15. a pin read only by a step the run did not run is left out: an opt-in step not in + TOOLS (run-all.sh writes no row for it) and one recorded as skipped (no BAM); + a status file with no step rows at all leaves every pin in; + 16. a malformed Last-Modified of a PGS file, and an offline release with no tag, + are errors of that lookup, never a crash; + 17. --refresh-clinvar refuses (exit 3) while a run holds the ClinVar lock shared, + and holds it exclusive while setup.sh refreshes; + 18. the ClinVar lock block of run-all.sh and of step 06 waits while a refresh holds + the lock and goes on once it ends (needs flock; it fails in CI without it). Each check is then run against a planted mistake and must fail there. Run: python3 tests/test_owngenome_recheck.py @@ -272,6 +286,8 @@ def all_checks(m): cmp["redone"] == ["06", "07"] and cmp["predicted_not_redone"] == ["13"] and cmp["redone_not_predicted"] == ["07"] and cmp["not_redone_script_steps"] == ["14"], cmp) + new_checks(m, tmp, steps, out) + ld = os.path.join(tmp, "launch") os.makedirs(os.path.join(ld, ".nextflow")) with open(os.path.join(ld, ".nextflow", "history"), "w") as f: @@ -290,6 +306,218 @@ def all_checks(m): return out +def manifest_add(g, rows): + with open(os.path.join(g, "S1", "run_manifest.tsv"), "a") as f: + f.writelines(f"version\t{k}\t{v}\n" for k, v in rows) + + +def step_rows(g, rows): + """run_status.tsv step rows as run-all.sh writes them; a step it left out gets no row.""" + with open(os.path.join(g, "S1", "logs", "run_status.tsv"), "a") as f: + f.writelines(f"step\t{k}\t{v}\n" for k, v in rows) + + +def new_checks(m, tmp, steps, out): + # Image references built from the fixture's own, so no image name is a literal here. + ref = THIS_PINS["PHARMCAT_IMAGE"] + d1, d2 = ref + "@sha256:" + "a" * 64, ref + "@sha256:" + "b" * 64 + g = genome(os.path.join(tmp, "digest")) + manifest_add(g, [("MEGANE_IMAGE", d1)]) + new = by_item(m.check("S1", g, answers(m, g), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0")) + same = by_item(m.check("S1", g, answers(m, g), steps, dict(THIS_PINS, MEGANE_IMAGE=d1), "3.66.0")) + mg, ms = new.get("MEGANE_IMAGE", {}), same.get("MEGANE_IMAGE", {}) + out["a new digest at the same tag of a pin with no lookup is rerun-ready"] = ( + mg.get("class") == "rerun-ready" and "50" in mg.get("steps", []) and mg.get("installed_ref") == d2 + and ms.get("class") == "no-newer-release", (mg, ms)) + + gp = genome(os.path.join(tmp, "digest-pharmcat")) + rows = [x for x in open(os.path.join(gp, "S1", "run_manifest.tsv")) if not x.startswith("version\tPHARMCAT_IMAGE\t")] + open(os.path.join(gp, "S1", "run_manifest.tsv"), "w").writelines(rows + [f"version\tPHARMCAT_IMAGE\t{d1}\n"]) + ph = by_item(m.check("S1", gp, answers(m, gp), steps, dict(THIS_PINS, PHARMCAT_IMAGE=d2), "3.66.0"))["PHARMCAT_IMAGE"] + ph1 = by_item(m.check("S1", gp, answers(m, gp), steps, dict(THIS_PINS, PHARMCAT_IMAGE=d1), "3.66.0"))["PHARMCAT_IMAGE"] + text = io.StringIO() + m.print_check({"sample": "S1", "checked_utc": "x", "items": [ph], "affected_steps": []}, text) + out["a new digest at the same tag of a looked-up image is rerun-ready"] = ( + ph["class"] == "rerun-ready" and ph1["class"] == "no-newer-release" and "same tag, new image" in text.getvalue(), + (ph, ph1, text.getvalue())) + + gu = genome(os.path.join(tmp, "unknown")) + st = m.check("S1", gu, answers(m, gu), steps, dict(THIS_PINS, MEGANE_IMAGE=d1), "3.66.0") + mu = by_item(st).get("MEGANE_IMAGE", {}) + with contextlib.redirect_stdout(io.StringIO()): + ans = write(tmp, "unknown.json", json.load(open(os.path.join(FIX, "offline.json")))) + keep = m.check + m.check = lambda s_, g_, a_, st_, pins, hla: keep(s_, g_, a_, st_, dict(pins, MEGANE_IMAGE=d1), hla) + try: + rc = m.main(["--sample", "S1", "--genome-dir", gu, "--offline", ans]) + finally: + m.check = keep + out["a pin with no lookup the run's manifest lacks is unknown, exit 0"] = ( + mu.get("class") == "unknown" and "50" not in st["affected_steps"] and "50" not in st["rerun_ready_steps"] + and "older version" in mu.get("action", "") and rc == 0, (mu, st["affected_steps"], rc)) + + # A default run (no TOOLS): run-all.sh writes rows for the steps it ran and none for step 50. + gd = genome(os.path.join(tmp, "not-in-tools")) + manifest_add(gd, [("MEGANE_IMAGE", d1)]) + step_rows(gd, [("06", "planned"), ("06", "ok")]) + nd = m.check("S1", gd, answers(m, gd), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0") + out["a pin read only by an opt-in step not in TOOLS is left out"] = ( + "MEGANE_IMAGE" not in by_item(nd) and "50" not in nd["rerun_ready_steps"], by_item(nd).get("MEGANE_IMAGE")) + + gs = genome(os.path.join(tmp, "skipped")) + manifest_add(gs, [("MEGANE_IMAGE", d1)]) + step_rows(gs, [("06", "ok"), ("50", "skipped (no BAM)")]) + sk = m.check("S1", gs, answers(m, gs), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0") + gr = genome(os.path.join(tmp, "ran")) + manifest_add(gr, [("MEGANE_IMAGE", d1)]) + step_rows(gr, [("06", "ok"), ("50", "planned"), ("50", "ok")]) + rn = by_item(m.check("S1", gr, answers(m, gr), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0")) + out["a pin read only by a step the run skipped is left out"] = ( + "MEGANE_IMAGE" not in by_item(sk) and "50" not in sk["rerun_ready_steps"] + and rn.get("MEGANE_IMAGE", {}).get("class") == "rerun-ready", (by_item(sk).get("MEGANE_IMAGE"), rn.get("MEGANE_IMAGE"))) + + srv, port = serve_pgs("not a date") + try: + fr = m.load_freshness() + fr.RETRIES = 0 + on = m.Online(fr, {}, "o/r", pgs_api=f"http://127.0.0.1:{port}/rest/score", raw=f"http://127.0.0.1:{port}", + out=io.StringIO()) + try: + gm = genome(os.path.join(tmp, "lastmod")) + pg = by_item(m.check("S1", gm, Mixed(on, answers(m, gm)), steps, THIS_PINS, "3.66.0"))["pgs:PGS000001"] + got = (pg["class"] == "error" and "Last-Modified" in pg.get("error", ""), pg) + except Exception as e: # noqa: BLE001 the crash this check is about + got = (False, repr(e)) + finally: + srv.shutdown() + out["a malformed Last-Modified of a PGS file is error, not a crash"] = got + + try: + go = genome(os.path.join(tmp, "notag")) + so = m.check("S1", go, answers(m, go, release={"pins": {}}), steps, THIS_PINS, "3.66.0") + got = ("error" in so["pipeline_release"] and "tag" in so["pipeline_release"]["error"], so["pipeline_release"]) + except Exception as e: # noqa: BLE001 + got = (False, repr(e)) + out["an offline release with no tag is error, not a crash"] = got + + gl = genome(os.path.join(tmp, "cvlock")) + lock = os.path.join(gl, "clinvar", ".refresh.lock") + keep_sp, ran = m.subprocess, [] + m.subprocess = type("SP", (), {"run": staticmethod(lambda *a, **k: ran.append(a) or type("R", (), {"returncode": 0})())}) + e = io.StringIO() + try: + with open(lock, "a") as held: + fcntl.flock(held, fcntl.LOCK_SH | fcntl.LOCK_NB) + with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(e): + rc = m.main(["--sample", "S1", "--genome-dir", gl, "--offline", os.path.join(FIX, "offline.json"), + "--rerun", "--refresh-clinvar"]) + finally: + m.subprocess = keep_sp + out["--refresh-clinvar refuses while a run holds the ClinVar lock"] = ( + rc == 3 and ".refresh.lock" in e.getvalue() and not ran, (rc, e.getvalue(), ran)) + + seen = [] + + def setup_run(argv, *a, **k): + if "--refresh" in argv: + with open(lock, "a") as f: + try: + fcntl.flock(f, fcntl.LOCK_SH | fcntl.LOCK_NB) + seen.append("free") + except OSError: + seen.append("held") + return type("R", (), {"returncode": 0})() + m.subprocess = type("SP", (), {"run": staticmethod(setup_run)}) + try: + with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()): + rc = m.main(["--sample", "S1", "--genome-dir", gl, "--offline", os.path.join(FIX, "offline.json"), + "--rerun", "--refresh-clinvar"]) + with open(lock, "a") as f: # let go before the replay's run-all.sh takes it shared + fcntl.flock(f, fcntl.LOCK_SH | fcntl.LOCK_NB) + after = True + except OSError: + after = False + finally: + m.subprocess = keep_sp + out["--refresh-clinvar holds the ClinVar lock while setup.sh refreshes"] = ( + seen == ["held"] and after and rc == 0, (seen, after, rc)) + + +def lock_block(path): + """The ClinVar lock block of a script: from 'CVLOCK=' through its closing 'fi'.""" + x = re.search(r"^CVLOCK=.*?\n fi\nfi\n", open(path).read(), re.S | re.M) + return x.group(0) if x else "" + + +def lock_waits(block): + """block, run in bash while this process holds the lock exclusive: it must wait, then go on + once the lock is let go.""" + import subprocess + import time + tmp = tempfile.mkdtemp(prefix="recheck-cvlock.") + try: + g = os.path.join(tmp, "genome") + os.makedirs(os.path.join(g, "clinvar")) + lock = os.path.join(g, "clinvar", ".refresh.lock") + with open(lock, "a") as held: + fcntl.flock(held, fcntl.LOCK_EX | fcntl.LOCK_NB) + p = subprocess.Popen(["bash", "-c", "set -euo pipefail\n" + block + "echo through\n"], + env=dict(os.environ, G=g, GENOME_DIR=g), stdout=subprocess.PIPE, text=True) + time.sleep(1.5) + waited = p.poll() is None + try: + o, _ = p.communicate(timeout=20) + except subprocess.TimeoutExpired: + p.kill() + o, _ = p.communicate() + return (bool(block) and waited and p.returncode == 0 and "Waiting" in o and o.rstrip().endswith("through"), + (waited, p.returncode, o)) + finally: + shutil.rmtree(tmp, ignore_errors=True) + + +def serve_pgs(last_modified): + """A PGS Catalog score answer, and a HEAD of its file with this Last-Modified.""" + class H(http.server.BaseHTTPRequestHandler): + def do_GET(self): + port = self.server.server_address[1] + body = json.dumps({"id": "PGS000001", "ftp_harmonized_scoring_files": {"GRCh38": { + "positions": f"http://127.0.0.1:{port}/f/PGS000001_hmPOS_GRCh38.txt.gz"}}}).encode() + self.send_response(200) + self.send_header("Content-Type", "application/json") + self.send_header("Content-Length", str(len(body))) + self.end_headers() + self.wfile.write(body) + + def do_HEAD(self): + self.send_response(200) + self.send_header("Last-Modified", last_modified) + self.send_header("Content-Length", "0") + self.end_headers() + + def log_message(self, *a): + pass + srv = http.server.ThreadingHTTPServer(("127.0.0.1", 0), H) + threading.Thread(target=srv.serve_forever, daemon=True).start() + return srv, srv.server_address[1] + + +def from_source(old, new): + """bin/owngenome_recheck.py with old replaced by new (which must be in it), as a module.""" + import types + src = open(TOOL).read() + if src.count(old) != 1: + raise AssertionError(f"planted mistake not in the copy: {old[:60]!r}") + mod = types.ModuleType("owngenome_recheck_planted") + mod.__file__ = TOOL + keep, sys.dont_write_bytecode = sys.dont_write_bytecode, True + try: + exec(compile(src.replace(old, new), TOOL, "exec"), mod.__dict__) + finally: + sys.dont_write_bytecode = keep + return mod + + class Mixed: """The online release and PGS lookups, the offline answers for the rest.""" @@ -428,6 +656,32 @@ def regex_resume(d): return cmd if re.search(r"(^|\s)-resume(\s|$)", cmd) else cmd + " -resume" m.launch_dir_command = regex_resume ms["--launch-dir: an argument holding -resume is not the option"] = m + ms["a new digest at the same tag of a pin with no lookup is rerun-ready"] = from_source( + '"rerun-ready" if ran_ref != inst_ref else', '"rerun-ready" if ran != inst else') + ms["a new digest at the same tag of a looked-up image is rerun-ready"] = from_source( + "and ran_ref != inst_ref:\n return \"rerun-ready\"", "and False:\n return \"rerun-ready\"") + ms["a pin with no lookup the run's manifest lacks is unknown, exit 0"] = from_source( + 'cls = "unknown" if not', 'cls = "no-newer-release" if not') + ms["a pin read only by an opt-in step not in TOOLS is left out"] = from_source( + 'if s["state"] in ("skipped", "not planned")}', 'if s["state"] in ("skipped",)}') + m = load() + m.steps_not_run = lambda *a: set() + ms["a pin read only by a step the run skipped is left out"] = m + ms["a malformed Last-Modified of a PGS file is error, not a crash"] = from_source( + " try: # a malformed date is a failed lookup of this score, never a crash of the whole check\n" + " return email.utils.parsedate_to_datetime(lm).strftime(\"%Y%m%d\")\n" + " except (TypeError, ValueError, AttributeError, OverflowError):\n", + " return email.utils.parsedate_to_datetime(lm).strftime(\"%Y%m%d\")\n if False:\n") + ms["an offline release with no tag is error, not a crash"] = from_source( + 'tag, pins = e.get("tag"), e.get("pins") or {}\n if not isinstance(tag, str) or not tag or not isinstance(pins, dict):', + 'tag, pins = e["tag"], e.get("pins") or {}\n if not isinstance(pins, dict):') + m = load() + m.fcntl = type("F", (), {"LOCK_EX": 0, "LOCK_NB": 0, "flock": staticmethod(lambda *a: None)}) + ms["--refresh-clinvar refuses while a run holds the ClinVar lock"] = m + ms["--refresh-clinvar holds the ClinVar lock while setup.sh refreshes"] = from_source( + " fcntl.flock(lock, fcntl.LOCK_EX | fcntl.LOCK_NB)\n except OSError:\n" + " print(f\"ERROR: a run or step 06", + " pass\n except OSError:\n print(f\"ERROR: a run or step 06") return ms @@ -595,6 +849,19 @@ def main(): unfiltered = plugin_from_source(src.replace(stale, "if False:"), "owngenome_live_recheck_unfiltered") check("fails on a planted mistake: the live-site plugin leaves out a check older than the run's manifest", not plugin_stale_ok(unfiltered)[0], "passed on the mutant") + import shutil as _sh + for script in ("scripts/run-all.sh", "scripts/06-clinvar-screen.sh"): + block = lock_block(os.path.join(REPO, script)) + if not _sh.which("flock"): + check(f"the ClinVar lock block of {script} waits for a refresh (flock is not installed: " + f"{'FAIL in CI' if os.environ.get('CI') else 'not run here'})", not os.environ.get("CI")) + continue + ok, detail = lock_waits(block) + check(f"the ClinVar lock block of {script} waits for a refresh, then goes on", ok, detail) + planted = re.sub(r"(?m)^ flock -s -n 7 .*$", " :", block) + check(f"planted mistake is in the copy: the ClinVar lock block of {script}", planted != block) + check(f"fails on a planted mistake: the ClinVar lock block of {script} waits for a refresh (no flock call)", + not lock_waits(planted)[0], "passed on the mutant") results = all_checks(load()) for name, (ok, detail) in results.items(): check(name, ok, detail)