From d3314d807233d29b339af2edea3dc9e0d183893b Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Mon, 9 Nov 2026 10:56:30 +0100 Subject: [PATCH 1/3] test(chip): read step 40's 23andMe file back with its chrM sites A new e2e case runs step 40 on the fixture, reads the 23andMe file back through chip-to-vcf.sh and checks every called chrM site lands on its rCRS position with case 21's call. On main it fails: chip-to-vcf.sh moves the MT positions as if they were on hg19's chrM. Epic: pgp-pfw6 --- scripts/ci/e2e-case-paths.tsv | 1 + scripts/ci/e2e-case-times.tsv | 1 + tests/e2e/chip-mt-1-roundtrip.sh | 164 +++++++++++++++++++++++++++++++ 3 files changed, 166 insertions(+) create mode 100755 tests/e2e/chip-mt-1-roundtrip.sh diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 59ebfe45..084cfefc 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 +chip-mt-1-roundtrip assets/par_grch38.bed bin/export_consumer.py scripts/40-export-consumer.sh scripts/chip-to-vcf.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz 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 diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index 8901ff00..be90f627 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 - +chip-mt-1-roundtrip shard 60 - ci-e2e-2-locityper shard 158 - clinical-filter-mane-and-quality-1-tiers shard 14 - clinical-filter-mane-and-quality-2-vep-norm shard 2 - diff --git a/tests/e2e/chip-mt-1-roundtrip.sh b/tests/e2e/chip-mt-1-roundtrip.sh new file mode 100755 index 00000000..5eafeeda --- /dev/null +++ b/tests/e2e/chip-mt-1-roundtrip.sh @@ -0,0 +1,164 @@ +#!/usr/bin/env bash +# Step 40's 23andMe-layout file read back by scripts/chip-to-vcf.sh, chrM +# included. Step 40 writes dbSNP's rCRS positions for the MT sites; a vendor +# 23andMe file numbers them on hg19's chrM, which chip-to-vcf.sh moves to the +# rCRS. It knows step 40's file by its first line ("# This data file was +# generated by OwnGenome") and keeps those positions. +# 1. step 40 on a copy of case 21's gVCF, with the CI slice of the site table +# (tests/fixtures/consumer_export/sites.tsv.gz, 285 chrM sites); +# 2. chip-to-vcf.sh on the 23andMe file says it kept the MT positions, and +# every called chrM site lands on the table's rCRS position with the +# reference base there and the allele of case 21's call (hom-ref where +# case 21 called nothing); +# 3. control: the same file without its first line, as a vendor file has +# no such line, is moved through the chain's chrM block, so most called +# chrM sites no longer land on their rCRS position with that call. +# Before the fix, check 2 failed the same way. +. "$(dirname "$0")/lib.sh" + +G="$GENOME_DIR" +GVCF="${G}/${SAMPLE}/vcf/${SAMPLE}.g.vcf.gz" +check "case 21 wrote a gVCF" vcf_ok "${SAMPLE}/vcf/${SAMPLE}.g.vcf.gz" +SLICE="${REPO}/tests/fixtures/consumer_export/sites.tsv.gz" +SITES="${G}/reference/consumer/23andme_v5_sites_dbsnp${CONSUMER_DBSNP_BUILD}.tsv.gz" +OWN_SITES=0 +if [ ! -e "$SITES" ]; then + mkdir -p "$(dirname "$SITES")" + cp "$SLICE" "$SITES" + OWN_SITES=1 +fi +CHAIN_URL=https://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz +CHAIN_SHA256=5c0598e500ceb5a78c73086929e8ef993aec309bcafb595139b53d440b125a1d +E=HG002CMT + +# --- 1. step 40 on a copy of the gVCF ------------------------------------------ +rm -rf "${G:?}/${E}" +mkdir -p "${G}/${E}/vcf" +cp "$GVCF" "${G}/${E}/vcf/${E}.g.vcf.gz" +cp "${GVCF}.tbi" "${G}/${E}/vcf/${E}.g.vcf.gz.tbi" +run_step 40-export-consumer.sh "$E" male +check_step_exit 40-export-consumer.sh +python3 -c 'import sys, zipfile; z = zipfile.ZipFile(sys.argv[1]); open(sys.argv[2], "wb").write(z.read(z.namelist()[0]))' \ + "${G}/${E}/consumer_export/${E}_23andme_v5.zip" "${CASE_TMP}/export.txt" +check "the 23andMe file starts with step 40's first line" \ + has '^# This data file was generated by OwnGenome' "$(head -n 1 "${CASE_TMP}/export.txt")" +N_MT=$(awk -F'\t' '!/^#/ && $2 == "MT"' "${CASE_TMP}/export.txt" | wc -l | tr -d ' ') +check_ge "chrM sites in the 23andMe file" "$N_MT" 200 + +# --- the GRCh37 side chip-to-vcf.sh reads ---------------------------------------- +# A FASTA with the table's GRCh37 base at each site (N elsewhere), as in +# consumer-export-1, and UCSC's chain. +mkdir -p "${G}/reference_hg19" "${G}/liftover" +python3 - "$SITES" "${G}/reference_hg19/human_g1k_v37.fasta" <<'PY' +import gzip, sys +sites, fa = sys.argv[1:3] +order = [str(i) for i in range(1, 23)] + ["X", "Y", "MT"] +bases = {} +for l in gzip.open(sites, "rt"): + if l.startswith("#"): + continue + c = l.rstrip("\n").split("\t") + bases.setdefault(c[1], {})[int(c[2])] = c[3] +off = 0 +with open(fa, "w") as f, open(fa + ".fai", "w") as fai: + for chrom in [c for c in order if c in bases]: + n = max(bases[chrom]) + 1 + head = f">{chrom}\n" + f.write(head) + off += len(head) + fai.write(f"{chrom}\t{n}\t{off}\t60\t61\n") + seq = bytearray(b"N" * n) + for p, b in bases[chrom].items(): + seq[p - 1] = ord(b) + for i in range(0, n, 60): + f.write(seq[i:i + 60].decode() + "\n") + off += n + (n + 59) // 60 +PY +if curl -fsSL --retry 3 -o "${G}/liftover/hg19ToHg38.over.chain.gz" "$CHAIN_URL"; then + check_eq "UCSC's hg19ToHg38 chain has the pinned sha256" \ + "$(sha256sum "${G}/liftover/hg19ToHg38.over.chain.gz" | cut -d' ' -f1)" "$CHAIN_SHA256" +else + fail "could not download ${CHAIN_URL}" +fi +bcf query -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER[\t%GT]\n' "${SAMPLE}/vcf/${SAMPLE}.vcf.gz" 2>/dev/null | awk -F'\t' '$1 == "chrM"' > "${CASE_TMP}/calls.tsv" + +# mt_roundtrip NAME TXT: chip-to-vcf.sh on TXT as a 23andMe file of sample +# NAME, its log in CASE_TMP/NAME.log; prints "called got on_place bad": the +# called chrM sites of TXT, the chrM records with a call, those on the table's +# rCRS position with its reference base, and those of them whose allele is +# not case 21's call there. +mt_roundtrip() { + local n=$1 + rm -rf "${G:?}/${n}" + mkdir -p "${G}/${n}/raw" + cp "$2" "${G}/${n}/raw/${n}_raw.txt" + "${REPO}/scripts/chip-to-vcf.sh" "$n" 23andme > "${CASE_TMP}/${n}.log" 2>&1 + echo "chip-to-vcf.sh ${n}: exit $?" >&2 + tail -n 25 "${CASE_TMP}/${n}.log" >&2 + bcf query -f '%ID\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' "${n}/vcf/${n}.vcf.gz" 2>/dev/null | awk -F'\t' '$2 == "chrM"' > "${CASE_TMP}/${n}_chrM.tsv" + python3 - "${CASE_TMP}/${n}_chrM.tsv" "${CASE_TMP}/calls.tsv" "$SITES" "$2" <<'PY' +import gzip, re, sys +lifted, calls, sites, txt = sys.argv[1:5] +table = {} +for l in gzip.open(sites, "rt"): + if not l.startswith("#"): + c = l.rstrip("\n").split("\t") + if c[1] == "MT": + table[c[0]] = c +dv = {} +for l in open(calls): + chrom, pos, ref, alt, flt, gt = l.rstrip("\n").split("\t") + if flt == "PASS": + al = [ref] + alt.split(",") + dv[int(pos)] = {al[int(i)][0] for i in re.split(r"[/|]", gt) if i not in (".", "")} +called = sum(1 for l in open(txt) + if not l.startswith("#") and l.split("\t")[0] in table and l.rstrip("\n").split("\t")[3] != "--") +got = on_place = bad = 0 +for l in open(lifted): + rid, chrom, pos, ref, alt, gt = l.rstrip("\n").split("\t") + t = table.get(rid) + idx = [i for i in re.split(r"[/|]", gt) if i not in (".", "")] + if t is None or not idx: + continue + got += 1 + if int(pos) != int(t[4]) or ref != t[5]: + if got - on_place <= 5: + print(f"off place {rid}: rCRS {t[4]} {t[5]}, read back at chrM:{pos} {ref}", file=sys.stderr) + continue + on_place += 1 + al = [ref] + alt.split(",") + have = {al[int(i)][0] for i in idx} + want = dv.get(int(pos), {t[5]}) + if have != want: + bad += 1 + if bad <= 5: + print(f"mismatch {rid} chrM:{pos} read back {sorted(have)}, case 21 {sorted(want)}", file=sys.stderr) +print(called, got, on_place, bad) +PY +} + +# --- 2. step 40's file keeps its chrM positions -------------------------------------- +read -r called got on_place bad <<< "$(mt_roundtrip HG002CMR "${CASE_TMP}/export.txt")" +echo "step 40's file: ${called:-?} called chrM sites, ${got:-?} read back, ${on_place:-?} on their rCRS position, ${bad:-?} with another allele than case 21's call" +check "chip-to-vcf.sh says it kept the MT positions of step 40's file" \ + has "written by OwnGenome step 40, ${N_MT} site\(s\) already on rCRS positions, kept" "$(cat "${CASE_TMP}/HG002CMR.log")" +check "chip-to-vcf.sh moved no MT site of step 40's file" lacks 'moved from hg19 chrM' "$(cat "${CASE_TMP}/HG002CMR.log")" +check_ge "called chrM sites in step 40's file" "${called:-0}" 100 +check_eq "called chrM sites read back" "${got:-x}" "${called:-y}" +check_eq "called chrM sites read back on their rCRS position with its reference base" "${on_place:-x}" "${called:-y}" +check_eq "chrM sites read back with another allele than case 21's call" "${bad:-x}" 0 + +# --- 3. control: the same file without step 40's first line ---------------------------- +tail -n +2 "${CASE_TMP}/export.txt" > "${CASE_TMP}/vendor.txt" +read -r v_called v_got v_on v_bad <<< "$(mt_roundtrip HG002CMV "${CASE_TMP}/vendor.txt")" +echo "control, no first line: ${v_called:-?} called chrM sites, ${v_got:-?} read back, ${v_on:-?} on their rCRS position, ${v_bad:-?} with another allele" +check "control: chip-to-vcf.sh moves the MT sites of a file without step 40's first line" \ + has 'site\(s\) moved from hg19 chrM to rCRS positions' "$(cat "${CASE_TMP}/HG002CMV.log")" +check "control: fewer than half of the called chrM sites stay on their rCRS position with case 21's call (${v_on:-?} of ${v_called:-?})" \ + test "$(( (${v_on:-0} - ${v_bad:-0}) * 2 ))" -lt "${v_called:-0}" +printf '#### chrM round trip of step 40\n\n%s called chrM sites, %s read back on their rCRS position, %s mismatches; without the first line %s stay\n\n' \ + "${called:-?}" "${on_place:-?}" "${bad:-?}" "${v_on:-?}" >> "$E2E_NOTES" + +rm -rf "${G:?}/${E}" "${G:?}/HG002CMR" "${G:?}/HG002CMV" "${G:?}/reference_hg19" "${G:?}/liftover" +[ "$OWN_SITES" = 0 ] || rm -f "$SITES" +finish From 14f2fe45ec4d2072f945cabce18ee2960c10280d Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Mon, 9 Nov 2026 11:00:40 +0100 Subject: [PATCH 2/3] fix(chip): keep the rCRS chrM positions of a file step 40 wrote chip-to-vcf.sh moves a 23andMe file's MT positions from hg19's chrM to the rCRS. Step 40's file already has rCRS positions, so reading it back put each chrM site one or two bases off. The script now knows that file by its first line and keeps its MT positions; a vendor file has no such line and moves as before. Epic: pgp-pfw6 --- docs/40-export-consumer.md | 2 +- docs/chip-data-guide.md | 2 +- scripts/chip-to-vcf.sh | 14 +++++++++++--- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/docs/40-export-consumer.md b/docs/40-export-consumer.md index c281cf47..fa4ba9d3 100644 --- a/docs/40-export-consumer.md +++ b/docs/40-export-consumer.md @@ -81,4 +81,4 @@ On the CI fixture (the gVCF of the fixture's 15 regions, 7,623 sites of the CI t - The AncestryDNA file has the 23andMe v5 sites, not the sites of an AncestryDNA chip; the two chips share part of their sites. A service that matches on the sites it expects uses the ones the files have in common. - The 23andMe file has no internal `i` sites and none of the dropped sites, so it has fewer rows than a file from 23andMe. - In a duplicated region the list's GRCh37 position and dbSNP's GRCh38 position can sit on two different copies. The file then carries the list's position with the genotype read at dbSNP's. On the CI table this is one site of 7,338 off chrM (in SMN1/SMN2), found by the round trip below. -- `scripts/chip-to-vcf.sh` reads both files back. The e2e case `consumer-export-1` does that on the CI fixture, lifts the result to GRCh38 and checks it against the fixture's calls; chrM is left out of both files there. `chip-to-vcf.sh` treats the MT positions of a 23andMe-layout file as 23andMe numbers them, on hg19's chrM, and moves each one to the rCRS through the chain's chrM block. This file already writes rCRS positions, so its MT sites would land one or two bases off. The AncestryDNA file's MT positions are read as rCRS positions and keep their place, but the case checks neither. +- `scripts/chip-to-vcf.sh` reads both files back. The e2e case `consumer-export-1` does that on the CI fixture, lifts the result to GRCh38 and checks it against the fixture's calls, chrM left out. `chip-to-vcf.sh` treats the MT positions of a 23andMe file as 23andMe numbers them, on hg19's chrM, and moves each one to the rCRS through the chain's chrM block. This file already writes rCRS positions, so `chip-to-vcf.sh` keeps its MT positions when the first line is the one this step writes (`# This data file was generated by OwnGenome`); without that line they would land one or two bases off. The e2e case `chip-mt-1-roundtrip` reads the 23andMe file back and checks that every called chrM site lands on its rCRS position with the fixture's call, and that the same file without that line is moved, as a vendor file is. The AncestryDNA file's MT positions are read as rCRS positions and keep their place. diff --git a/docs/chip-data-guide.md b/docs/chip-data-guide.md index 8741c2be..031f53ce 100644 --- a/docs/chip-data-guide.md +++ b/docs/chip-data-guide.md @@ -212,7 +212,7 @@ echo "Done. VCF at: ${GENOME_DIR}/${SAMPLE}/vcf/${SAMPLE}.vcf.gz" **Expected output:** A bgzipped, tabix-indexed VCF on GRCh38 coordinates with ~600,000 SNPs. Typical breakdown: ~430K hom-ref, ~107K het, ~66K hom-alt, ~1K missing. Some variants (~0.3%) are rejected during liftover; ~900 have swapped REF/ALT between builds. -> **X/Y/MT chromosomes:** All chromosomes are converted. MT is renamed to chrM to match GRCh38 convention and keeps its positions: GRCh37's MT and GRCh38's chrM are both the rCRS, so it is not lifted (the hg19 chain maps hg19's own chrM, the Yoruba sequence, and would move m.1555 to 1553). MyHeritage and AncestryDNA MT positions are read as rCRS positions, as the GRCh37 build they name. 23andMe numbers its MT sites on hg19's chrM instead (its v5 site list puts m.3243 at 3244), so for a 23andMe file `chip-to-vcf.sh` first moves each MT position to the rCRS through the chrM block of `hg19ToHg38.over.chain.gz`, and drops a site that falls in one of that chain's gaps, printing how many it moved and dropped. The 23andMe-layout file of [step 40](40-export-consumer.md) writes rCRS MT positions, unlike 23andMe, so its MT sites read back through this script land one or two bases off; the round trip there leaves chrM out. Chip arrays cover very few mtDNA positions, so for mitochondrial haplogroup estimation, dedicated tools like [HaploGrep](https://haplogrep.i-med.ac.at/) that accept raw 23andMe files directly will give better results. +> **X/Y/MT chromosomes:** All chromosomes are converted. MT is renamed to chrM to match GRCh38 convention and keeps its positions: GRCh37's MT and GRCh38's chrM are both the rCRS, so it is not lifted (the hg19 chain maps hg19's own chrM, the Yoruba sequence, and would move m.1555 to 1553). MyHeritage and AncestryDNA MT positions are read as rCRS positions, as the GRCh37 build they name. 23andMe numbers its MT sites on hg19's chrM instead (its v5 site list puts m.3243 at 3244), so for a 23andMe file `chip-to-vcf.sh` first moves each MT position to the rCRS through the chrM block of `hg19ToHg38.over.chain.gz`, and drops a site that falls in one of that chain's gaps, printing how many it moved and dropped. The 23andMe-layout file of [step 40](40-export-consumer.md) writes rCRS MT positions, unlike 23andMe: `chip-to-vcf.sh` knows it by its first line (`# This data file was generated by OwnGenome`) and keeps those positions. Chip arrays cover very few mtDNA positions, so for mitochondrial haplogroup estimation, dedicated tools like [HaploGrep](https://haplogrep.i-med.ac.at/) that accept raw 23andMe files directly will give better results. > **23andMe / AncestryDNA:** Place your file at `${GENOME_DIR}/${SAMPLE}/raw/${SAMPLE}_raw.txt`. Both are tab-separated but not in the same layout: AncestryDNA needs its own pre-step above (two allele columns, numeric chromosome codes), which `chip-to-vcf.sh` runs for you. AncestryDNA reports two alleles on X, Y and MT, so its genotypes there are diploid; 23andMe reports one allele there for a male. diff --git a/scripts/chip-to-vcf.sh b/scripts/chip-to-vcf.sh index 2e43f1e2..fb788541 100755 --- a/scripts/chip-to-vcf.sh +++ b/scripts/chip-to-vcf.sh @@ -94,14 +94,21 @@ case "$FORMAT" in # list puts m.3243 at 3244. Each MT position moves to the rCRS through the # chrM chain of the hg19ToHg38 chain file (hg38's chrM is the rCRS). A # site inside one of that chain's gaps has no rCRS position and is - # dropped, with its count printed. + # dropped, with its count printed. A file step 40 wrote + # (bin/export_consumer.py; its first line starts "# This data file was + # generated by OwnGenome") already has dbSNP's rCRS MT positions, so its + # MT sites stay where they are: moving them would put each one or two + # bases off. A vendor file has no such line and moves as above. + RCRS=0 + if head -n 1 "$RAW_TSV" | grep -q '^# This data file was generated by OwnGenome'; then RCRS=1; fi MT_BLOCKS="${RAW_DIR}/${SAMPLE}_chrM_blocks.tsv" gzip -dc "$CHAIN" | awk -v OFS='\t' ' $1 == "chain" { inm = ($3 == "chrM" && $8 == "chrM" && $5 == "+" && $10 == "+"); t = $6; q = $11; next } inm && $1 ~ /^[0-9]+$/ { print t, q, $1; t += $1 + $2; q += $1 + $3 }' > "$MT_BLOCKS" - awk -F'\t' -v OFS='\t' -v blocks="$MT_BLOCKS" ' + awk -F'\t' -v OFS='\t' -v blocks="$MT_BLOCKS" -v rcrs="$RCRS" ' BEGIN { while ((getline l < blocks) > 0) { split(l, b, "\t"); n++; T[n] = b[1]; Q[n] = b[2]; S[n] = b[3] } } - $2 == "MT" { + $2 == "MT" && rcrs { kept++ } + $2 == "MT" && !rcrs { mt++ p = $3 + 0; to = 0 for (i = 1; i <= n; i++) if (p > T[i] && p <= T[i] + S[i]) { to = Q[i] + p - T[i]; break } @@ -110,6 +117,7 @@ case "$FORMAT" in } { print } END { + if (kept) printf " 23andMe MT: written by OwnGenome step 40, %d site(s) already on rCRS positions, kept\n", kept > "/dev/stderr" if (mt && !n) printf " 23andMe MT: no chrM chain in the chain file; %d MT site(s) dropped\n", mt > "/dev/stderr" else if (mt) printf " 23andMe MT: %d site(s) moved from hg19 chrM to rCRS positions, %d in a chain gap dropped\n", mt - gap, gap > "/dev/stderr" }' "$CHIP_TSV" > "${CHIP_TSV}.rcrs" From cc2b7776b82f377d217bedff89e4746b63d31767 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Mon, 9 Nov 2026 16:36:17 +0100 Subject: [PATCH 3/3] fix(chip): know step 40's 23andMe file when an editor added a byte order mark A step 40 file saved again by an editor that writes a UTF-8 byte order mark lost its first-line marker, so chip-to-vcf.sh moved its chrM sites off their rCRS positions and kept the marked comment line as data. Drop the mark before both reads. The e2e case reads such a file back, and its control now also asserts the conversion exits 0 and reads chrM back. --- scripts/chip-to-vcf.sh | 7 +++++-- tests/e2e/chip-mt-1-roundtrip.sh | 22 +++++++++++++++++++--- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/scripts/chip-to-vcf.sh b/scripts/chip-to-vcf.sh index fb788541..776180d2 100755 --- a/scripts/chip-to-vcf.sh +++ b/scripts/chip-to-vcf.sh @@ -88,7 +88,10 @@ case "$FORMAT" in exit 1 fi # rsid, chromosome, position, genotype already; drop the comment header. - grep -v '^#' "$RAW_TSV" | tr -d '\r' > "$CHIP_TSV" + # A file saved again by an editor that adds a UTF-8 byte order mark keeps + # its first line a comment: the mark is dropped before both reads. + BOM=$'\xEF\xBB\xBF' + LC_ALL=C sed "1s/^${BOM}//" "$RAW_TSV" | grep -v '^#' | tr -d '\r' > "$CHIP_TSV" # 23andMe numbers its MT sites on hg19's chrM (the Yoruba sequence, # 16,571 bases), not on the rCRS that stage 1 reads REF from: its v5 site # list puts m.3243 at 3244. Each MT position moves to the rCRS through the @@ -100,7 +103,7 @@ case "$FORMAT" in # MT sites stay where they are: moving them would put each one or two # bases off. A vendor file has no such line and moves as above. RCRS=0 - if head -n 1 "$RAW_TSV" | grep -q '^# This data file was generated by OwnGenome'; then RCRS=1; fi + if head -n 1 "$RAW_TSV" | LC_ALL=C sed "s/^${BOM}//" | grep -q '^# This data file was generated by OwnGenome'; then RCRS=1; fi MT_BLOCKS="${RAW_DIR}/${SAMPLE}_chrM_blocks.tsv" gzip -dc "$CHAIN" | awk -v OFS='\t' ' $1 == "chain" { inm = ($3 == "chrM" && $8 == "chrM" && $5 == "+" && $10 == "+"); t = $6; q = $11; next } diff --git a/tests/e2e/chip-mt-1-roundtrip.sh b/tests/e2e/chip-mt-1-roundtrip.sh index 5eafeeda..3cf97f06 100755 --- a/tests/e2e/chip-mt-1-roundtrip.sh +++ b/tests/e2e/chip-mt-1-roundtrip.sh @@ -10,6 +10,8 @@ # every called chrM site lands on the table's rCRS position with the # reference base there and the allele of case 21's call (hom-ref where # case 21 called nothing); +# 2b. the same file with a UTF-8 byte order mark in front, as an editor +# that adds one saves it, is still known as step 40's file; # 3. control: the same file without its first line, as a vendor file has # no such line, is moved through the chain's chrM block, so most called # chrM sites no longer land on their rCRS position with that call. @@ -92,8 +94,10 @@ mt_roundtrip() { rm -rf "${G:?}/${n}" mkdir -p "${G}/${n}/raw" cp "$2" "${G}/${n}/raw/${n}_raw.txt" - "${REPO}/scripts/chip-to-vcf.sh" "$n" 23andme > "${CASE_TMP}/${n}.log" 2>&1 - echo "chip-to-vcf.sh ${n}: exit $?" >&2 + local status=0 + "${REPO}/scripts/chip-to-vcf.sh" "$n" 23andme > "${CASE_TMP}/${n}.log" 2>&1 || status=$? + echo "$status" > "${CASE_TMP}/${n}.status" + echo "chip-to-vcf.sh ${n}: exit ${status}" >&2 tail -n 25 "${CASE_TMP}/${n}.log" >&2 bcf query -f '%ID\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' "${n}/vcf/${n}.vcf.gz" 2>/dev/null | awk -F'\t' '$2 == "chrM"' > "${CASE_TMP}/${n}_chrM.tsv" python3 - "${CASE_TMP}/${n}_chrM.tsv" "${CASE_TMP}/calls.tsv" "$SITES" "$2" <<'PY' @@ -140,6 +144,7 @@ PY # --- 2. step 40's file keeps its chrM positions -------------------------------------- read -r called got on_place bad <<< "$(mt_roundtrip HG002CMR "${CASE_TMP}/export.txt")" echo "step 40's file: ${called:-?} called chrM sites, ${got:-?} read back, ${on_place:-?} on their rCRS position, ${bad:-?} with another allele than case 21's call" +check_eq "chip-to-vcf.sh exits 0 on step 40's file" "$(cat "${CASE_TMP}/HG002CMR.status")" 0 check "chip-to-vcf.sh says it kept the MT positions of step 40's file" \ has "written by OwnGenome step 40, ${N_MT} site\(s\) already on rCRS positions, kept" "$(cat "${CASE_TMP}/HG002CMR.log")" check "chip-to-vcf.sh moved no MT site of step 40's file" lacks 'moved from hg19 chrM' "$(cat "${CASE_TMP}/HG002CMR.log")" @@ -148,10 +153,21 @@ check_eq "called chrM sites read back" "${got:-x}" "${called:-y}" check_eq "called chrM sites read back on their rCRS position with its reference base" "${on_place:-x}" "${called:-y}" check_eq "chrM sites read back with another allele than case 21's call" "${bad:-x}" 0 +# --- 2b. the same file with a UTF-8 byte order mark --------------------------------- +{ printf '\357\273\277'; cat "${CASE_TMP}/export.txt"; } > "${CASE_TMP}/bom.txt" +read -r b_called b_got b_on b_bad <<< "$(mt_roundtrip HG002CMB "${CASE_TMP}/bom.txt")" +echo "byte order mark: ${b_called:-?} called chrM sites, ${b_got:-?} read back, ${b_on:-?} on their rCRS position, ${b_bad:-?} with another allele" +check "chip-to-vcf.sh keeps the MT positions of step 40's file with a byte order mark" \ + has "written by OwnGenome step 40, ${N_MT} site\(s\) already on rCRS positions, kept" "$(cat "${CASE_TMP}/HG002CMB.log")" +check_eq "byte order mark: called chrM sites read back on their rCRS position" "${b_on:-x}" "${called:-y}" +check_eq "byte order mark: chrM sites with another allele than case 21's call" "${b_bad:-x}" 0 + # --- 3. control: the same file without step 40's first line ---------------------------- tail -n +2 "${CASE_TMP}/export.txt" > "${CASE_TMP}/vendor.txt" read -r v_called v_got v_on v_bad <<< "$(mt_roundtrip HG002CMV "${CASE_TMP}/vendor.txt")" echo "control, no first line: ${v_called:-?} called chrM sites, ${v_got:-?} read back, ${v_on:-?} on their rCRS position, ${v_bad:-?} with another allele" +check_eq "control: chip-to-vcf.sh exits 0" "$(cat "${CASE_TMP}/HG002CMV.status")" 0 +check_ge "control: called chrM records read back" "${v_got:-0}" 1 check "control: chip-to-vcf.sh moves the MT sites of a file without step 40's first line" \ has 'site\(s\) moved from hg19 chrM to rCRS positions' "$(cat "${CASE_TMP}/HG002CMV.log")" check "control: fewer than half of the called chrM sites stay on their rCRS position with case 21's call (${v_on:-?} of ${v_called:-?})" \ @@ -159,6 +175,6 @@ check "control: fewer than half of the called chrM sites stay on their rCRS posi printf '#### chrM round trip of step 40\n\n%s called chrM sites, %s read back on their rCRS position, %s mismatches; without the first line %s stay\n\n' \ "${called:-?}" "${on_place:-?}" "${bad:-?}" "${v_on:-?}" >> "$E2E_NOTES" -rm -rf "${G:?}/${E}" "${G:?}/HG002CMR" "${G:?}/HG002CMV" "${G:?}/reference_hg19" "${G:?}/liftover" +rm -rf "${G:?}/${E}" "${G:?}/HG002CMR" "${G:?}/HG002CMB" "${G:?}/HG002CMV" "${G:?}/reference_hg19" "${G:?}/liftover" [ "$OWN_SITES" = 0 ] || rm -f "$SITES" finish