diff --git a/bin/collect_summary.py b/bin/collect_summary.py index d8de0d4e..4f59a0f8 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -763,6 +763,35 @@ def parascopy_cn_at(rows, chrom, pos): return {"cn": None, "text": "no row covers exon 7", "agcn": ".", "agcn_qual": ".", "agcn_filter": "."} +MARKER_FIELDS = ("marker", "rsid", "status", "alt_reads", "depth", "fraction", "copies", "message") + + +def read_marker_table(p, step, key): + """A copy-aware marker table (metric, value; bin/traits_report.py MARKERS): + step 46's KIV-2 4925G>A beside the count or step 35's SMN1 g.27134T>G + beside the psCN. A table that is not there or has no known status reads + as not assessed, with the reason in its message, so the report line + never goes missing beside the count.""" + m = traits_report.MARKERS[key] + what = m["name"] + out = dict.fromkeys(MARKER_FIELDS, ".") + out.update(marker=what, rsid=m["dbsnp"][0], status="not assessed") + if not os.path.isfile(p): + out["message"] = (f"{what}: not assessed, step {step} wrote no {os.path.basename(p)} (a run from before " + "this reader, or the check did not run).") + return out + try: + kv = dict(read_tsv_pairs(p)) + except OSError as e: + kv = {"message": f"{os.path.basename(p)} could not be read: {e}"} + if kv.get("status") not in traits_report.MARKER_STATUSES: + out["message"] = (f"{what}: not assessed, {os.path.basename(p)} has no status row " + f"({', '.join(traits_report.MARKER_STATUSES)}).") + return out + out.update({k: kv[k] for k in MARKER_FIELDS if kv.get(k)}) + return out + + def smn_reading(smn1): """What an SMN1 copy number means for SMA, in one line.""" if smn1 is None: @@ -881,6 +910,9 @@ def sec_parascopy(d, s): catalog = {r["table"]: r for r in asset_rows("paralog_loci.tsv")} if others else {} return p, {"smn1": smn1, "smn2": {"cn": smn2["cn"], "text": smn2["text"]}, "agcn": smn1["agcn"], "agcn_qual": smn1["agcn_qual"], "agcn_filter": smn1["agcn_filter"], + # SMN1 g.27134T>G beside the psCN (step 35, _smn1_marker.tsv) + "smn1_marker": read_marker_table(os.path.join(d, "paralogs", f"{s}_smn1_marker.tsv"), "35", + "SMN1_27134"), "rows": len(rows), # HBA1/HBA2 from depth (bin/hba_copy_number.py); None when # step 35 wrote no HBA table (a run from before the check) @@ -2786,6 +2818,9 @@ def sec_lpa(d, s): raise Unreadable(f"{os.path.basename(p)} has no status row (estimated, too little depth or not assessed)") data = {k: kv.get(k) or "." for k in ("status", "kiv2_copies", "ci95_low", "ci95_high", "class", "kiv2_depth", "flank_depth", "mapq0_share", "message", "advice")} + # KIV-2 4925G>A beside the count (step 46, lpa/_lpa_kiv2_4925.tsv) + data["null_allele"] = read_marker_table(os.path.join(os.path.dirname(p), f"{s}_lpa_kiv2_4925.tsv"), "46", + "LPA_KIV2_4925") variants = [] t = first_existing(d, [f"traits/{s}_traits.tsv"]) for r in (read_tsv(t) if t else []): diff --git a/bin/live_results.py b/bin/live_results.py index 74f77a53..f1a4d0aa 100644 --- a/bin/live_results.py +++ b/bin/live_results.py @@ -301,6 +301,8 @@ def a_lpa(S, s): d = S["lpa"]["data"] facts = [["KIV-2 repeats", f"{d.get('kiv2_copies')} (95% interval {d.get('ci95_low')} to {d.get('ci95_high')})" if d.get("status") == "estimated" else d.get("status")], ["Class", d.get("class")]] + if d.get("null_allele"): + facts.append(["KIV-2 4925G>A null allele", rr.marker_state(d["null_allele"])]) return {"facts": facts, "lines": rr.lpa_lines(d), "attention": d.get("flag") == "high"} @@ -357,15 +359,19 @@ def a_couple(S, s): def a_paralogs(S, s): - lines = [] + lines, facts = [], [] if (S.get("parascopy") or {}).get("state") in ("ok", "stale"): d = S["parascopy"]["data"] lines += rr.parascopy_lines(d) + rr.parascopy_locus_lines(d) + rr.cn_interval_lines(d.get("cn_intervals")) + if d.get("smn1_marker"): + # the silent-carrier marker beside the SMN1 psCN + facts = [["SMN1 copies", (d.get("smn1") or {}).get("text", ".")], + ["SMN1 g.27134T>G silent-carrier marker", rr.marker_state(d["smn1_marker"])]] rows = [] if (S.get("paralogs") or {}).get("state") in ("ok", "stale"): rows = [[g.get("gene"), g.get("copy_number"), g.get("status"), g.get("detail")] for g in S["paralogs"]["data"].get("genes") or []] - return {"lines": lines, "table": table(["Gene", "Copies", "Status", "Detail"], rows)} + return {"facts": facts, "lines": lines, "table": table(["Gene", "Copies", "Status", "Detail"], rows)} def a_pharmcat(S, s): diff --git a/bin/render_report.py b/bin/render_report.py index 0e17c2db..25a29f2f 100755 --- a/bin/render_report.py +++ b/bin/render_report.py @@ -119,6 +119,23 @@ def ancient_lines(d): "health on its own.") LPA_NOTE = ("The KIV-2 repeat count from read depth is a research estimate, not an Lp(a) level: the level is a " "blood test, and the count explains only part of the difference between people.") +LPA_NULL_CAVEAT = ("A null allele lowers the Lp(a) made from the copies it sits on; the Lp(a) blood test is the " + "measure.") +SMN1_MARKER_CAVEAT = ("The marker raises the chance that two SMN1 copies sit on one chromosome (2+0, a silent carrier); " + "it does not prove 2+0, and a carrier test in a clinical laboratory confirms it.") + + +def marker_state(m): + """A copy-aware marker's state in a few words (collect_summary.read_marker_table): + present or uncertain with the reads and the copies that carry it, not seen + with the reads, too little depth, or not assessed.""" + st = m.get("status") or "not assessed" + if st == "not assessed": + return st + out = f"{st}, {m.get('alt_reads', '.')} of {m.get('depth', '.')} reads carry it" + if m.get("copies") not in (None, "", "."): + out += f", about {m['copies']} copies" + return out def combined_cell(h): @@ -1799,7 +1816,11 @@ def paralog_lines(d): def parascopy_lines(d): """Step 35's lines, the same in both reports.""" smn1, smn2 = d.get("smn1") or {}, d.get("smn2") or {} - return [f"SMN1 copies: {smn1.get('text', '.')}", + mk = d.get("smn1_marker") + marker = ([f"SMN1 g.27134T>G silent-carrier marker ({mk.get('rsid') or '.'}): {marker_state(mk)}", + f"SMN1 marker reading: {mk.get('message', '.')}", + f"SMN1 marker caveat: {SMN1_MARKER_CAVEAT}"] if mk else []) + return [f"SMN1 copies: {smn1.get('text', '.')}"] + marker + [ f"SMN2 copies: {smn2.get('text', '.')}", f"SMN1 and SMN2 together (agCN) at exon 7: {d.get('agcn', '.')} (quality {d.get('agcn_qual', '.')}, " f"filter {d.get('agcn_filter', '.')})", @@ -2146,6 +2167,10 @@ def lpa_lines(d): else: out = [f"KIV-2 repeats: {d.get('status', '.')}"] out.append(f" {d.get('message', '.')}") + na = d.get("null_allele") + if na: + out += [f"KIV-2 4925G>A null allele ({(na.get('rsid') or '.').split(' ')[0]}): {marker_state(na)}", + f" {na.get('message', '.')}", f" {LPA_NULL_CAVEAT}"] for v in d.get("variants", []): n = v["risk_alleles"] out.append(f"{v['rsid']}: {v['genotype']} ({'not read' if n is None else str(n) + ' risk allele' + ('s' if n != 1 else '')})") @@ -2897,7 +2922,8 @@ def a(key, lines): if S.get("parascopy", MISSING)["state"] != "missing": a("parascopy", card(S["parascopy"], "SMN1/SMN2 Copy Number (Parascopy)", [stat("Status", done_badge(S["parascopy"]))] - + [stat(*(E(x) for x in line.split(": ", 1))) for line in parascopy_lines(S["parascopy"]["data"])] + # stat escapes the label itself; the value is escaped here + + [stat(lab, E(val)) for lab, val in (line.split(": ", 1) for line in parascopy_lines(S["parascopy"]["data"]))] + parascopy_loci_html(S["parascopy"]["data"]) + cn_intervals_html(S["parascopy"]["data"].get("cn_intervals")) + [small_note(PARASCOPY_NOTE)])) diff --git a/docs/24-html-report.md b/docs/24-html-report.md index 2ec5d054..7306b6b6 100644 --- a/docs/24-html-report.md +++ b/docs/24-html-report.md @@ -96,6 +96,7 @@ The cards, by chapter (cards not listed here are named in the chapter table abov - **Structural Variants** — Manta, Delly, CNVpytor and the SURVIVOR consensus (2+ callers) counts; then step 05's class 4 and 5 SVs (AnnotSV) with their caller support, exon-edge and frequency checks. A call step 05's rules moved to class 3 (a coding tandem repeat, an immunoglobulin or T-cell receptor locus, the MHC, a pseudogene named by a disease gene's alias, or read depth that does not show the change) is listed apart with AnnotSV's class, the rule and the reason, and the class 4-5 count leaves it out - **Repeat Expansions** — a table of every locus, the five key loci (HTT, FMR1, C9orf72, ATXN1, DMPK) first, with its repeat count, Stranger's status and notes; and counts of the loci outside the normal range (step 9b, Stranger), of those short reads could not size (not counted as normal), and the loci in an intermediate range (`bin/str_intermediate_ranges.tsv`). When Stranger flags RFC1, or another locus where the motif decides (`expansion_hunter/${SAMPLE}_str_motifs.tsv`), the card adds the verdict of the read-level motif count (step 9b). A locus the repeat catalog gives no thresholds (its `NoThresholds` field in `expansion_hunter/${SAMPLE}_eh_catalog.json`: MUC1, CEL, TYMS and PLIN4 with the STRchive catalog, among others) is marked not graded, with the catalog's reason, never normal. Any other locus outside Stranger's catalog has no status and is counted apart. Without step 9b the card says Stranger did not run; a Stranger file older than the ExpansionHunter result it came from is not read, and the card says to rerun step 9b. Repeat-expansion calls from short reads can be wrong at some loci - **Mitochondrial** — chrM PASS variants and heteroplasmic calls (allele fraction 0.05 to 0.95; below 5% NUMT reads and noise dominate) +- **Lp(a): LPA KIV-2 Repeat Count** — when step 46 ran: the KIV-2 count with its interval and class, the KIV-2 4925G>A null allele under the count with its caveat (`not assessed` without step 46's 4925G>A table), step 42's two LPA risk variants, the flag and the blood-test advice - **Mosaic Loss of Chromosome Y or X** — step 33's call on chrY (men) or chrX (women), the reason for it, and a note that this is a population-level association, not a diagnosis **Carrier** @@ -103,7 +104,7 @@ The cards, by chapter (cards not listed here are named in the chapter table abov - **Carrier List (this sample)** — the step 06 ClinVar hits in the recessive and X-linked genes of `assets/carrier_genes.tsv`, for this sample alone; hypomorph and pseudodeficiency alleles set aside - **Couple Carrier Screen** — when step 39 ran for this sample: the partner, the genes where both carry a change, the X-linked findings and the genes where only one partner does, each shared or X-linked gene with the chance for a child - **Paralog Genes (Paraphase)** — when step 35b ran: copy number and call per gene, SMN1 from Paraphase beside Parascopy, with whether they agree, when step 35 ran too, and the count of variants on phased haplotypes -- **SMN1/SMN2 Copy Number (Parascopy)** — when step 35 ran: SMN1 and SMN2 copies at exon 7 with quality and filter, the two together (agCN), the SMA reading, and, when step 35 wrote its HBA table, the count of alpha-globin genes from depth or the reason it was not assessed +- **SMN1/SMN2 Copy Number (Parascopy)** — when step 35 ran: SMN1 and SMN2 copies at exon 7 with quality and filter, the SMN1 g.27134T>G silent-carrier marker under the SMN1 copies with its caveat (`not assessed` without step 35's marker table), the two together (agCN), the SMA reading, and, when step 35 wrote its HBA table, the count of alpha-globin genes from depth or the reason it was not assessed - **Paralog Copy Number** — `not checked`, gray, when neither step 35 nor step 35b left a result, naming a step that failed **Medicines** diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md index 51a450e5..a45b3850 100644 --- a/docs/35-paralogs.md +++ b/docs/35-paralogs.md @@ -96,7 +96,7 @@ SMN2 holds the same 21 bases 875,420 bases upstream, so the reads cannot tell an | `too little depth` | under 10 reads over the two positions | | `not assessed` | the pileup could not be read, or the check failed | -The message reads the marker beside the SMN1 psCN at exon 7 (filter `PASS` and quality 20 or more, the rule of the report's SMN1 line): with two SMN1 copies it says whether the marker raises or lowers the chance of 2+0; with another count, or a psCN that is not resolved, it says the marker refines a two-copy result only. Every message ends with "This marker does not prove 2+0, and its absence does not rule 2+0 out." The reports do not read this table yet; it sits beside the SMN table in the same folder. +The message reads the marker beside the SMN1 psCN at exon 7 (filter `PASS` and quality 20 or more, the rule of the report's SMN1 line): with two SMN1 copies it says whether the marker raises or lowers the chance of 2+0; with another count, or a psCN that is not resolved, it says the marker refines a two-copy result only. Every message ends with "This marker does not prove 2+0, and its absence does not rule 2+0 out." Both reports and the live person page print the marker right under the SMN1 copies: its state (`present` or `uncertain` with the reads that carry G and the copies, `not seen` with the reads, `too little depth`, or `not assessed`), the reader's message, and the caveat that the marker raises the chance that two SMN1 copies sit on one chromosome, does not prove 2+0, and that a carrier test in a clinical laboratory confirms it. Without the table (a run from before this reader) or with a table that has no known status, the line says `not assessed` with the reason. The live card also lists the SMN1 copies and the marker state as facts. The Nextflow `HTML_REPORT` gets the table from the `PARASCOPY` task the same way. ## Runtime diff --git a/docs/42-traits.md b/docs/42-traits.md index 1a644b35..7c6affb7 100644 --- a/docs/42-traits.md +++ b/docs/42-traits.md @@ -115,8 +115,8 @@ The six rules after TAS2R38 need no phase. Their sites lie close together on one Two more markers of the planned list are not rules here, because each sits in sequence the genome holds more than once and a gVCF genotype at one position cannot read it. They are read from the BAM by the step that counts the copies, as a fraction of copies rather than a genotype (`MARKERS` in `bin/traits_report.py`, which pins each one's dbSNP b157 record and the GRCh38 sequence around every copy): -- **LPA KIV-2 4925G>A** (Coassin et al. 2017, PMID 28444229), a splice variant that lowers the Lp(a) made from the chromosome that carries it. GRCh38 holds the site in each of its six KIV-2 units and in KIV-1. [Step 46](46-lpa.md) reads it beside the KIV-2 count. -- **SMN1 g.27134T>G** (rs143838139; Luo et al. 2014, PMID 23788250), a marker of SMN1 duplication alleles that SMN2 holds too. [Step 35](35-paralogs.md) reads it beside the SMN1 copy number. It refines the chance of a silent 2+0 carrier for a two-copy result, and its line always says that it does not prove 2+0. +- **LPA KIV-2 4925G>A** (Coassin et al. 2017, PMID 28444229), a splice variant that lowers the Lp(a) made from the chromosome that carries it. GRCh38 holds the site in each of its six KIV-2 units and in KIV-1. [Step 46](46-lpa.md) reads it beside the KIV-2 count, and the reports print it under the count. +- **SMN1 g.27134T>G** (rs143838139; Luo et al. 2014, PMID 23788250), a marker of SMN1 duplication alleles that SMN2 holds too. [Step 35](35-paralogs.md) reads it beside the SMN1 copy number, and the reports print it under the SMN1 copies. It refines the chance of a silent 2+0 carrier for a two-copy result, and its line always says that it does not prove 2+0. Kim et al. (2003) found five TAS2R38 haplotypes worldwide. The rule reads the two common ones. A copy with AAV, AAI, PVI or a rarer one is `other_haplotype`, so the GRCh38 reference itself, which carries AAI (Ala49, Ala262, Ile296), is not read when every site matches it. The readings follow Bufe et al. (2005): PAV/PAV rated PTC and PROP the most bitter, PAV/AVI recognised them as bitter at nearly the same concentrations with a lower and more variable intensity, and AVI/AVI recognised PTC only at far higher concentrations. diff --git a/docs/46-lpa.md b/docs/46-lpa.md index 3c0c31f2..b1ead3ce 100644 --- a/docs/46-lpa.md +++ b/docs/46-lpa.md @@ -33,6 +33,7 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags. ## Input - `${SAMPLE}/aligned/${SAMPLE}_sorted.bam` with its `.bai` (`ALIGN_DIR` as for the other BAM steps), aligned to GRCh38 with `chr` contig names. +- A failure of the 4925G>A reader never costs the count: step 46 and the `LPA` process then write a 4925G>A table that says `not assessed, the check failed`, print a warning, and go on. - The aligner must report a read that fits several places equally well, with mapping quality 0. minimap2 (step 2) and BWA-MEM do. novoalign by default leaves such reads unaligned: then the step finds almost no MAPQ 0 reads over KIV-2 and reports the count as not assessed rather than a count that comes out low. ## Command @@ -105,7 +106,7 @@ What a short-read fraction cannot say: which chromosome or which repeat carries | `lpa/${SAMPLE}_lpa_kiv2_depth_q1.tsv` | the same, reads of MAPQ 1 or more | | `lpa/${SAMPLE}_lpa_kiv2_4925_pileup.tsv` | `samtools mpileup` at the seven 4925G>A positions | -Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. The reports do not read the 4925G>A table yet; step 46 prints its message beside the count. +Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. Right under the count, both reports and the live person page print the 4925G>A line: `present` or `uncertain` with the reads that carry it and the copies, `not seen` with the reads, `too little depth`, or `not assessed`; then the reader's message and the caveat that a null allele lowers the Lp(a) made from the copies it sits on and that the Lp(a) blood test is the measure. Without the 4925G>A table (a run from before this reader) or with a table that has no known status, the line says `not assessed` with the reason. The live card also lists the state as a fact. The Nextflow `HTML_REPORT` gets the table from the `LPA` process the same way. ## Runtime diff --git a/main.nf b/main.nf index 7d09f5ff..b52aaa1a 100644 --- a/main.nf +++ b/main.nf @@ -956,12 +956,14 @@ workflow { // Read depth over the KIV-2 exons against single-copy LPA exons, for each // sample with a BAM. Equivalent to scripts/46-lpa.sh. ch_lpa = Channel.empty() + ch_lpa_null = Channel.empty() ch_lpa_versions = Channel.empty() if (tools_list.contains('lpa')) { LPA_REGIONS() LPA_DEPTH(ch_bam, LPA_REGIONS.out.bed) LPA(LPA_DEPTH.out.depth) ch_lpa = LPA.out.table + ch_lpa_null = LPA.out.null_allele ch_lpa_versions = LPA_REGIONS.out.versions.mix(LPA_DEPTH.out.versions, LPA.out.versions) } @@ -1165,6 +1167,7 @@ workflow { .join(BAM_ANALYSIS.out.paraphase_vars.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.smn_copy_number.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.hba_copy_number.map { meta, f -> [meta.id, f] }, remainder: true) + .join(BAM_ANALYSIS.out.smn1_marker.map { meta, f -> [meta.id, f] }, remainder: true) // step 35's other catalog tables, a list per sample (flattened below) .join(BAM_ANALYSIS.out.parascopy_tables.map { meta, fs -> [meta.id, fs] }, remainder: true) .join(CLINICAL.out.prs_scores.map { meta, f -> [meta.id, f] }, remainder: true) @@ -1211,6 +1214,7 @@ workflow { .join(ch_ancient_info.map { meta, f -> [meta.id, f] }, remainder: true) .join(ch_ancient_svg.map { meta, f -> [meta.id, f] }, remainder: true) .join(ch_lpa.map { meta, f -> [meta.id, f] }, remainder: true) + .join(ch_lpa_null.map { meta, f -> [meta.id, f] }, remainder: true) .join(ch_vntyper.map { meta, f -> [meta.id, f] }, remainder: true) .join(ch_viral.map { meta, f -> [meta.id, f] }, remainder: true) .join(ch_exomiser.map { meta, f -> [meta.id, f] }, remainder: true) diff --git a/modules/local/html_report/main.nf b/modules/local/html_report/main.nf index a7842255..845f470f 100644 --- a/modules/local/html_report/main.nf +++ b/modules/local/html_report/main.nf @@ -25,8 +25,9 @@ (the confirmed mtDNA variants step 20's section names) and the paralog catalogs assets/paralog_loci.tsv and assets/paraphase_genes.tsv. Every __copy_number.tsv of step 35 goes to paralogs/, and so - do step 35's _cn_intervals.tsv, _cn_dosage.tsv and - _cyp21a2_psv.tsv. The sample's gVCF and its index, named in + do step 35's _cn_intervals.tsv, _cn_dosage.tsv, + _cyp21a2_psv.tsv and _smn1_marker.tsv; step 46's _lpa_kiv2.tsv + and _lpa_kiv2_4925.tsv go to lpa/. The sample's gVCF and its index, named in meta.gvcf because a samplesheet gVCF keeps the user's name, become vcf/.g.vcf.gz and its .tbi (or .csi), where the bash report finds step 03's gVCF. A run_manifest.tsv in the sample folder records the @@ -102,6 +103,7 @@ process HTML_REPORT { ${id}_sample_qc.tsv|${id}_read_count.tsv) place qc "\$f" ;; ${id}_cn_intervals.tsv) place paralogs "\$f" ;; ${id}_cn_dosage.tsv|${id}_cyp21a2_psv.tsv) place paralogs "\$f" ;; + ${id}_smn1_marker.tsv) place paralogs "\$f" ;; ${id}_prs_summary.tsv) place prs "\$f" ;; ${id}_ancestry.tsv) place ancestry "\$f" ;; ${id}_ancestry_fine.tsv|${id}_painting.png) place ancestry_fine "\$f" ;; @@ -126,6 +128,7 @@ process HTML_REPORT { ${id}_archaic.tsv) place archaic "\$f" ;; ${id}_ancient.tsv|${id}_ancient_info.tsv|${id}_ancient_pca.svg) place ancient "\$f" ;; ${id}_lpa_kiv2.tsv) place lpa "\$f" ;; + ${id}_lpa_kiv2_4925.tsv) place lpa "\$f" ;; ${id}_vntyper.tsv) place vntyper "\$f" ;; ${id}_locityper.tsv) place locityper "\$f" ;; ${id}_viral.tsv) place viral "\$f" ;; diff --git a/modules/local/lpa/main.nf b/modules/local/lpa/main.nf index ac3229c5..97ff4ea1 100644 --- a/modules/local/lpa/main.nf +++ b/modules/local/lpa/main.nf @@ -13,7 +13,8 @@ LPA bin/copy_number.py call --locus LPA_KIV2: the count, its 95% interval and the message (laid out by bin/lpa_kiv2.py), then bin/traits_report.py --marker LPA_KIV2_4925: the 4925G>A - reading beside it, _lpa_kiv2_4925.tsv + reading beside it, _lpa_kiv2_4925.tsv ('not assessed' + when that check fails; the count stays) Equivalent to: scripts/46-lpa.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -123,8 +124,14 @@ process LPA { script: """ copy_number.py call --locus LPA_KIV2 --bedcov ${depth} --bedcov-mapq1 ${depth_q1} --table ${meta.id}_lpa_kiv2.tsv - traits_report.py --marker LPA_KIV2_4925 --pileup ${pileup} --copy-table ${meta.id}_lpa_kiv2.tsv \\ - --marker-out ${meta.id}_lpa_kiv2_4925.tsv + # KIV-2 4925G>A beside the count (as scripts/46-lpa.sh [4/4]): a failure + # never costs the count; the table then says not assessed. + if ! traits_report.py --marker LPA_KIV2_4925 --pileup ${pileup} --copy-table ${meta.id}_lpa_kiv2.tsv \\ + --marker-out ${meta.id}_lpa_kiv2_4925.tsv; then + echo "WARNING: the KIV-2 4925G>A check failed; ${meta.id}_lpa_kiv2_4925.tsv says not assessed." >&2 + printf 'metric\\tvalue\\nmarker\\tLPA KIV-2 4925G>A\\nrsid\\trs1484919544 (first KIV-2 copy)\\nstatus\\tnot assessed\\nalt_reads\\t.\\ndepth\\t.\\nfraction\\t.\\ncopies\\t.\\nkiv2_copies\\t.\\nmessage\\t%s\\n' \\ + 'LPA KIV-2 4925G>A: not assessed, the check failed; see the log of step 46.' > ${meta.id}_lpa_kiv2_4925.tsv + fi cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/scripts/46-lpa.sh b/scripts/46-lpa.sh index 7bbfe417..def07ded 100755 --- a/scripts/46-lpa.sh +++ b/scripts/46-lpa.sh @@ -28,7 +28,8 @@ # lpa/_lpa_kiv2_depth.tsv samtools bedcov -j over the regions (all MAPQ) # lpa/_lpa_kiv2_depth_q1.tsv the same, reads of MAPQ 1 or more # lpa/_lpa_kiv2_4925.tsv metricvalue: the 4925G>A reading -# (status, reads, fraction, copies, message) +# (status, reads, fraction, copies, message); +# 'not assessed' when the check fails # lpa/_lpa_kiv2_4925_pileup.tsv samtools mpileup at the seven # copies of the 4925G>A site (every MAPQ) set -euo pipefail @@ -91,9 +92,16 @@ run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ --bedcov-mapq1 "$(cpath "$DEPTH_Q1")" --table "$(cpath "$TABLE")" echo "[4/4] KIV-2 4925G>A beside the count (bin/traits_report.py --marker LPA_KIV2_4925)..." -run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ +# A failure here never costs the count (as step 35's SMN1 marker): the +# 4925G>A table then says not assessed, and the step goes on. The LPA +# process writes the same table. +if ! run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ python3 /pgp-bin/traits_report.py --marker LPA_KIV2_4925 --pileup "$(cpath "$PILEUP")" \ - --copy-table "$(cpath "$TABLE")" --marker-out "$(cpath "$NULL")" + --copy-table "$(cpath "$TABLE")" --marker-out "$(cpath "$NULL")"; then + echo "WARNING: the KIV-2 4925G>A check failed; ${NULL##*/} says not assessed." >&2 + printf 'metric\tvalue\nmarker\tLPA KIV-2 4925G>A\nrsid\trs1484919544 (first KIV-2 copy)\nstatus\tnot assessed\nalt_reads\t.\ndepth\t.\nfraction\t.\ncopies\t.\nkiv2_copies\t.\nmessage\t%s\n' \ + 'LPA KIV-2 4925G>A: not assessed, the check failed; see the log of step 46.' > "$NULL" +fi [ -s "$TABLE" ] || { echo "ERROR: no table written (${TABLE})" >&2; exit 1; } echo "" diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 32f1f190..e87f1e39 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -96,8 +96,8 @@ genome-lookup-1 scripts/ci/fake-docker/lib.sh scripts/owngenome-lookup.sh genotype-rule-rows-1-str-extra assets/str_extra_loci.json assets/stub/db/catalog.json scripts/09-expansion-hunter.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh genotype-rules-1-phase assets/genotype_rules.tsv assets/traits.tsv bin/traits_report.py modules/local/traits/* modules/local/vcf_precheck/* scripts/42-traits.sh scripts/ci/check-rule-alleles.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh genotype-rules-2-g2 assets/genotype_rules.tsv bin/collect_summary.py bin/genotype_rules_indels.py bin/traits_report.py -genotype-rules-2-lpa-null bin/cn_loci.tsv bin/copy_number.py bin/lpa_kiv2.py bin/traits_report.py modules/local/lpa/* scripts/46-lpa.sh -genotype-rules-2-smn1-marker bin/copy_number.py bin/traits_report.py modules/local/parascopy/* scripts/35-paralogs.sh +genotype-rules-2-lpa-null bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/lpa_kiv2.py bin/render_report.py bin/traits_report.py modules/local/lpa/* scripts/24-html-report.sh scripts/46-lpa.sh scripts/generate-report.sh +genotype-rules-2-smn1-marker bin/collect_summary.py bin/copy_number.py bin/render_report.py bin/traits_report.py modules/local/parascopy/* scripts/24-html-report.sh scripts/35-paralogs.sh scripts/generate-report.sh haplogroup-explore-1-files assets/callability_genes.tsv assets/live/index.html assets/mito_confirmed.tsv assets/viewer/index.html bin/collect_summary.py bin/y_tree_place.py bin/yleaf_run.py modules/local/indexcov/* modules/local/mito_haplogroup/* modules/local/mito_variants/* modules/local/vcf_precheck/* modules/local/y_haplogroup/* scripts/12-mito-haplogroup.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/setup.sh scripts/validate-setup.sh hba-short-read-accuracy-1-nf-report 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/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/render_report.py bin/results_index.py modules/local/export_consumer/* modules/local/html_report/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz hla-second-caller-1-bash 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/validate-setup.sh diff --git a/tests/e2e/genotype-rules-2-lpa-null.sh b/tests/e2e/genotype-rules-2-lpa-null.sh index 104fbed5..c287305c 100755 --- a/tests/e2e/genotype-rules-2-lpa-null.sh +++ b/tests/e2e/genotype-rules-2-lpa-null.sh @@ -20,7 +20,12 @@ # The reading goes to lpa/_lpa_kiv2_4925.tsv; step 46's own table # keeps its rows. # - the LPA module's pileup and 4925G>A lines, cut from its main.nf, write the -# same pileup and table in their images, byte for byte, on the planted BAM. +# same pileup and table in their images, byte for byte, on the planted BAM; +# with a 4925G>A reader that crashes, its block writes a table that says +# not assessed and leaves the count; +# - steps 24 and generate-report.sh on the planted and the reference samples: +# both reports print the reading under the count with its caveat +# ('present' with the reads and copies, 'not seen'). . "$(dirname "$0")/lib.sh" G="$GENOME_DIR" T="$CASE_TMP" @@ -128,6 +133,23 @@ check_eq "reference: not seen" "$(null "$R" status)" "not seen" check_eq "reference: no read with T" "$(null "$R" alt_reads)" 0 check "control: the planted and the reference tables differ" bash -c '! cmp -s "$1" "$2"' _ \ "${G}/${P}/lpa/${P}_lpa_kiv2_4925.tsv" "${G}/${R}/lpa/${R}_lpa_kiv2_4925.tsv" +# --- both reports print the reading under the count (steps 24 and generate-report.sh) ----- +for s in "$P" "$R"; do + run_step 24-html-report.sh "$s"; check_step_exit 24-html-report.sh + run_step generate-report.sh "$s"; check_step_exit generate-report.sh +done +rep_has() { grep -qF -- "$2" "${G}/$1/$1_report.$3"; } +check "planted: the HTML report prints the null allele as present, with reads and copies" \ + rep_has "$P" "KIV-2 4925G>A null allele (rs1484919544): present, $(null "$P" alt_reads) of $(null "$P" depth) reads carry it, about $(null "$P" copies) copies" html +check "planted: the text report prints it under the count" \ + bash -c 'sed -n "/KIV-2 repeats (both chromosomes)/,/rs10455872/p" "$1" | grep -qF "KIV-2 4925G>A null allele (rs1484919544): present"' _ "${G}/${P}/${P}_report.txt" +check "planted: both reports carry the caveat" \ + bash -c 'grep -qF "A null allele lowers the Lp(a) made from the copies it sits on" "$1" && grep -qF "A null allele lowers the Lp(a) made from the copies it sits on" "$2"' _ \ + "${G}/${P}/${P}_report.html" "${G}/${P}/${P}_report.txt" +check "reference: both reports print not seen (control: the line follows the table)" \ + bash -c 'grep -qF "KIV-2 4925G>A null allele (rs1484919544): not seen, 0 of" "$1" && grep -qF "KIV-2 4925G>A null allele (rs1484919544): not seen, 0 of" "$2"' _ \ + "${G}/${R}/${R}_report.html" "${G}/${R}/${R}_report.txt" + echo "genotype-rules-2-lpa-null: HG002 KIV-2 4925G>A $(null "$U" status)," \ "$(null "$U" alt_reads) of $(null "$U" depth) reads (fraction $(null "$U" fraction))," \ "KIV-2 count $(metric "$U" kiv2_copies); planted $(null "$P" fraction), copies $(null "$P" copies)" >> "$E2E_NOTES" @@ -140,9 +162,9 @@ unesc() { sed -e 's/\\\\/\\/g' -e 's/\\\$/$/g' -e "s|\${bam}|${G}/${P}/aligned/$ -e "s|\${pileup}|${P}_lpa_kiv2_4925_pileup.tsv|g"; } grep -E '^ *traits_report\.py --marker-regions LPA_KIV2_4925 > lpa_sites\.txt$' "$MOD" | unesc > "${M}/sites.sh" sed -n '/^ *( while read -r r; do samtools mpileup /,/|| echo /p' "$MOD" | unesc > "${M}/pileup.sh" -sed -n '/^ *traits_report\.py --marker LPA_KIV2_4925 --pileup/,/--marker-out/p' "$MOD" | unesc > "${M}/null.sh" +sed -n '/^ *if ! traits_report\.py --marker LPA_KIV2_4925 --pileup/,/^ *fi$/p' "$MOD" | unesc > "${M}/null.sh" check_eq "the module's sites, pileup and 4925G>A lines were cut from main.nf" \ - "$(cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" | wc -l | tr -d ' ')" 5 + "$(cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" | wc -l | tr -d ' ')" 9 cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" cp "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "$M/" mod() { docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -v "${REPO}/bin:/pgp-bin:ro" \ @@ -154,5 +176,22 @@ check "the module writes step 46's 4925G>A table, byte for byte" \ cmp -s "${G}/${P}/lpa/${P}_lpa_kiv2_4925.tsv" "${M}/${P}_lpa_kiv2_4925.tsv" check "the module leaves step 46's table byte for byte" cmp -s "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "${M}/${P}_lpa_kiv2.tsv" +# The module's fallback, in its image: a 4925G>A reader that crashes leaves the +# count and a table that says not assessed (scripts/46-lpa.sh writes the same +# table; tests/test_lpa_fallback.sh checks that). Control: the planted table +# above was 'present'. +mkdir -p "${M}/broken" +printf '#!/bin/sh\necho "planted crash" >&2\nexit 1\n' > "${M}/broken/traits_report.py" +chmod +x "${M}/broken/traits_report.py" +cp "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "${M}/count_before.tsv" +docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -v "${REPO}/bin:/pgp-bin:ro" \ + -w "$M" "$PYTHON_IMAGE" sh -c "PATH=\"${M}/broken:/pgp-bin:\$PATH\"; set -e; . ./null.sh" 2>&1 | tail -n 2 +check_eq "the module with a crashing 4925G>A reader: not assessed" \ + "$(awk -F'\t' '$1 == "status" {print $2}' "${M}/${P}_lpa_kiv2_4925.tsv")" "not assessed" +check "the module with a crashing 4925G>A reader: the message says the check failed" \ + grep -q "not assessed, the check failed" "${M}/${P}_lpa_kiv2_4925.tsv" +check "the module with a crashing 4925G>A reader: the count table stays" \ + cmp -s "${M}/count_before.tsv" "${M}/${P}_lpa_kiv2.tsv" + rm -rf "${G:?}/${U}" "${G:?}/${P}" "${G:?}/${R}" "${M:?}" "${T:?}"/*.sam "${T:?}"/*.bam "${T:?}"/*.bai finish diff --git a/tests/e2e/genotype-rules-2-smn1-marker.sh b/tests/e2e/genotype-rules-2-smn1-marker.sh index 5cdbaf03..a5804dbe 100755 --- a/tests/e2e/genotype-rules-2-smn1-marker.sh +++ b/tests/e2e/genotype-rules-2-smn1-marker.sh @@ -18,7 +18,10 @@ # from scripts/35-paralogs.sh, with scripts/lib/common.sh's run_in), and # the PARASCOPY module's block, cut from its main.nf, writes the same table # in its image, byte for byte, on the planted BAM. -# - every message says the marker does not prove 2+0. +# - every message says the marker does not prove 2+0; +# - steps 24 and generate-report.sh with the planted and the reference +# tables: both reports print the marker under the SMN1 copies with its +# caveat ('present' with its reads, 'not seen'). . "$(dirname "$0")/lib.sh" SRC="https://storage.googleapis.com/deepvariant/case-study-testdata/HG002.novaseq.pcr-free.35x.dedup.grch38_no_alt.bam" @@ -143,5 +146,24 @@ docker run --rm --user "$(id -u):$(id -g)" -v "${GENOME_DIR}:/genome" -v "${REPO "$PARASCOPY_IMAGE" sh -c 'PATH="/pgp-bin:$PATH"; sh -e block.sh' 2>&1 | tail -n 3 check "the module's block writes step 35's table, byte for byte" cmp -s "$L" "${M}/${P}_smn1_marker.tsv" +# --- both reports print the marker under the SMN1 copies (steps 24 and generate-report.sh) +rm -rf "$M" "${W}"/*.bam "${W}"/*.bai +for v in planted reference; do + cp "${CASE_TMP}/${v}.tsv" "${W}/paralogs/${P}_smn1_marker.tsv" + run_step 24-html-report.sh "$P"; check_step_exit 24-html-report.sh + run_step generate-report.sh "$P"; check_step_exit generate-report.sh + cp "${W}/${P}_report.html" "${CASE_TMP}/${v}.html"; cp "${W}/${P}_report.txt" "${CASE_TMP}/${v}.txt" +done +check "planted: the HTML report prints the marker as present, with its reads, under the SMN1 copies" \ + grep -qF "SMN1 g.27134T>G silent-carrier marker (rs143838139)present, $(col "$L" alt_reads) of $(col "$L" depth) reads carry it" "${CASE_TMP}/planted.html" +check "planted: the text report prints it right under the SMN1 copies" \ + bash -c 'grep -A1 "SMN1 copies: " "$1" | grep -qF "SMN1 g.27134T>G silent-carrier marker (rs143838139): present"' _ "${CASE_TMP}/planted.txt" +check "planted: both reports carry the caveat" \ + bash -c 'grep -qF "it does not prove 2+0, and a carrier test in a clinical laboratory confirms it" "$1" && grep -qF "it does not prove 2+0, and a carrier test in a clinical laboratory confirms it" "$2"' _ \ + "${CASE_TMP}/planted.html" "${CASE_TMP}/planted.txt" +check "reference: both reports print not seen (control: the line follows the table)" \ + bash -c 'grep -qF "SMN1 g.27134T>G silent-carrier marker (rs143838139)not seen, 0 of" "$1" && grep -qF "SMN1 g.27134T>G silent-carrier marker (rs143838139): not seen, 0 of" "$2"' _ \ + "${CASE_TMP}/reference.html" "${CASE_TMP}/reference.txt" + rm -rf "$W" finish diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index 5ef81ccc..4e28b325 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -132,7 +132,10 @@ folder with neither step 35 nor step 35b prints "Paralog copy number: not checked (neither step 35 nor step 35b ran)" in both reports, a failed step 35 is named in that line, and Paraphase alone leaves the line out; - a table without the psCN columns is unreadable. + a table without the psCN columns is unreadable. The SMN1 g.27134T>G + silent-carrier marker (step 35's _smn1_marker.tsv) prints beside the + SMN1 copies in both reports with its caveat: present with its reads, + not assessed without the table or with no known status in it. 26. step 08's opt-in passes in the HLA section: with HLA_QC=true, one line per class I gene from the concordance table pgx_outside_calls.py writes (agree, disagree with both calls, not compared with the reason), and @@ -1624,6 +1627,38 @@ def parascopy_section(work, schema): and "not a carrier by copy number" in card_s and "CYP21A2 and the gene-dosage loci have no Parascopy model" in card_s and "step 35 counts them from read depth against flanking sequence" in card_s and "Not checked" not in html, card_s) + # SMN1 g.27134T>G beside the psCN (step 35's _smn1_marker.tsv, bin/traits_report.py smn1_marker) + caveat = render_report.SMN1_MARKER_CAVEAT + after_cn = " SMN1 copies: 2 (quality 33, filter PASS)\n SMN1 g.27134T>G silent-carrier marker (rs143838139): " + check("parascopy: without step 35's marker table the marker line says not assessed, beside the psCN, with why", + (pd.get("smn1_marker") or {}).get("status") == "not assessed" + and after_cn + "not assessed\n" in txt and "step 35 wrote no S_smn1_marker.tsv" in txt + and f" SMN1 marker caveat: {caveat}\n" in txt, txt) + import traits_report + mk = f"{d}/paralogs/S_smn1_marker.tsv" + traits_report.write_metric_table(mk, traits_report.smn1_marker((35, 5, 0), f"{d}/paralogs/S_smn_copy_number.tsv")) + summ_m, txt_m, html_m = render(d) + card_m = html_m.split("

SMN1/SMN2 Copy Number (Parascopy)

", 1)[-1].split("

", 1)[0] + check("parascopy: a marker table with 5 of 40 reads reads as present beside the psCN in the text report, " + "with the reader's message and the caveat", + after_cn + "present, 5 of 40 reads carry it, about 0.5 copies\n" in txt_m + and " SMN1 marker reading: SMN1 g.27134T>G (rs143838139): present" in txt_m + and f" SMN1 marker caveat: {caveat}\n" in txt_m, txt_m) + check("parascopy: the HTML card shows the marker, its message and the caveat", + ">SMN1 g.27134T>G silent-carrier marker (rs143838139)<" in card_m + and "present, 5 of 40 reads carry it" in card_m and render_report.E(caveat) in card_m + and "does not prove 2+0" in card_m, card_m) + import live_results + lv = live_results.a_paralogs(summ_m["sections"], {}) + check("live: the paralog card has the SMN1 copies and the marker as facts, and its lines with the caveat", + ["SMN1 g.27134T>G silent-carrier marker", "present, 5 of 40 reads carry it, about 0.5 copies"] in lv["facts"] + and ["SMN1 copies", "2 (quality 33, filter PASS)"] in lv["facts"] + and f"SMN1 marker caveat: {caveat}" in lv["lines"], lv) + put(mk, "metric\tvalue\nstatus\tmaybe\n") + _, txt_b, _ = render(d) + check("parascopy: a marker table with no known status reads as not assessed, never as a result", + after_cn + "not assessed\n" in txt_b and "S_smn1_marker.tsv has no status row" in txt_b, txt_b) + os.remove(mk) check("parascopy: HBA1/HBA2 listed as not assessed without step 35's HBA table or step 35b", any(x.startswith("HBA1/HBA2 copy number") for x in summ["not_assessed"]) and "HBA1/HBA2 copy number (alpha-thalassaemia deletions): neither step 35's read-depth count" in txt diff --git a/tests/test_lpa_fallback.sh b/tests/test_lpa_fallback.sh new file mode 100755 index 00000000..acb0f7e2 --- /dev/null +++ b/tests/test_lpa_fallback.sh @@ -0,0 +1,112 @@ +#!/usr/bin/env bash +# Step 46 [4/4]: a KIV-2 4925G>A reader that crashes never costs the KIV-2 +# count. scripts/46-lpa.sh then writes lpa/_lpa_kiv2_4925.tsv saying +# 'not assessed', exits 0 and keeps lpa/_lpa_kiv2.tsv; the LPA Nextflow +# process (modules/local/lpa/main.nf) writes the same fallback table, byte for +# byte. The e2e case genotype-rules-2-lpa-null runs the module's block with a +# crashing reader in its image. +# +# A fake `docker` on PATH stands in for the containers: the count call writes +# a table, the 4925G>A call exits 1 (or 0 for the control, writing a +# 'present' table), everything else succeeds. +# +# Controls: with the reader working the table is the reader's ('present'), so +# the not-assessed check can fail; with the fallback removed from a copy of +# the step, the step exits non-zero, so the exit check can fail. +set -euo pipefail + +REPO=$(cd "$(dirname "$0")/.." && pwd) +WORK=$(mktemp -d) +trap 'rm -rf "$WORK"' EXIT + +FAILS=0 +fail() { echo "FAIL: $*"; FAILS=$((FAILS + 1)); } +pass() { echo "ok: $*"; } + +mkdir -p "${WORK}/bin" +cat > "${WORK}/bin/docker" <<'FAKE' +#!/usr/bin/env bash +printf '%s\n' "$*" >> "$FAKE_DOCKER_LOG" +host() { printf '%s' "${GENOME_DIR}${1#/genome}"; } +args=("$@") +for ((i = 0; i < ${#args[@]}; i++)); do + case "${args[i]}" in + --table) printf 'metric\tvalue\nstatus\testimated\nkiv2_copies\t38.0\n' > "$(host "${args[i + 1]}")" ;; + --marker-out) + if [ "${FAKE_READER:-crash}" = crash ]; then echo "planted crash" >&2; exit 1; fi + printf 'metric\tvalue\nstatus\tpresent\n' > "$(host "${args[i + 1]}")" ;; + esac +done +exit 0 +FAKE +chmod +x "${WORK}/bin/docker" + +GD="${WORK}/genome" +mkdir -p "${GD}/S1/aligned" "${GD}/reference" +: > "${GD}/S1/aligned/S1_sorted.bam" +: > "${GD}/S1/aligned/S1_sorted.bam.bai" +: > "${GD}/reference/GRCh38_no_alt_analysis_set.fasta" +: > "${GD}/reference/GRCh38_no_alt_analysis_set.fasta.fai" +LOG="${WORK}/docker.log" +NULL="${GD}/S1/lpa/S1_lpa_kiv2_4925.tsv" +COUNT="${GD}/S1/lpa/S1_lpa_kiv2.tsv" + +# run_step