diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml new file mode 100644 index 00000000..8d4318ab --- /dev/null +++ b/.github/workflows/public-evidence.yml @@ -0,0 +1,510 @@ +name: Public evidence + +# Runs two opt-in steps whole, on public data only (GIAB's HG002 and the +# public archaic and ancient panels), and writes what they measured to the +# job summary. Not a gate. +# +# archaic-sites the 1000 Genomes sites file setup.sh --archaic reads, +# downloaded once (its md5 checked) and indexed, for +# every chromosome job. +# archaic-hg002 setup.sh --archaic and step 43 on HG002's GIAB v4.2.1 +# genotypes, one job per autosome (about 150 GB of +# archaic downloads in all; each job fetches its three +# files side by side first, and setup.sh resumes them +# and checks their md5), then archaic-share joins the +# 22 into one report per population and prints HG002's +# Neanderthal and Denisovan shares, and how many of +# HG002's positions carry two records, one an SNV. +# The job "empty" is the red control: an empty region +# list must make it fail. +# ancient-footprint setup.sh --ancient-panel on the full AADR 1240K files, +# then step 44 on HG002, each with the peak memory of +# its containers (the cgroup's memory.peak), the peak +# disk use and the time. +# +# It runs by hand only (workflow_dispatch). Public repo: the GitHub-hosted +# runner is free here. +on: + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: false + +env: + NCBI_REF: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids + GIAB_HG002: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark + +jobs: + archaic-sites: + runs-on: ubuntu-latest + timeout-minutes: 120 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - name: The 1000 Genomes phase 3 sites file and its index + run: | + set -euo pipefail + . ./versions.env + url=$(sed -n 's/^ARCHAIC_1KG_SITES_URL=//p' scripts/setup.sh) + mkdir -p sites + curl -fsSL --retry 5 -o sites/1kg_sites.vcf.gz "$url" + [ "$(md5sum sites/1kg_sites.vcf.gz | cut -d' ' -f1)" = "$ARCHAIC_1KG_SITES_MD5" ] || { echo "md5 of ${url} is not ARCHAIC_1KG_SITES_MD5"; exit 1; } + docker run --rm -u "$(id -u):$(id -g)" -v "${PWD}/sites:/s" "$BCFTOOLS_IMAGE" bcftools index -t /s/1kg_sites.vcf.gz + ls -l sites + - uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 + with: + name: archaic-1kg-sites + path: sites + if-no-files-found: error + compression-level: 0 + retention-days: 3 + + archaic-hg002: + name: archaic-hg002 (${{ matrix.part }}) + needs: archaic-sites + runs-on: ubuntu-latest + timeout-minutes: 350 + strategy: + fail-fast: false + # Eight at a time, so the other pull requests keep runners. + max-parallel: 8 + matrix: + part: [empty, '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', + '18', '19', '20', '21', '22'] + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + + - name: Free disk space and pick the work area + run: | + sudo rm -rf /usr/local/lib/android /usr/share/dotnet /opt/ghc /usr/local/.ghcup \ + /opt/hostedtoolcache/CodeQL /usr/local/share/boost + df -h / /mnt + work="${RUNNER_TEMP}/ev" + if [ "$(df --output=avail -k /mnt | tail -1)" -gt "$(df --output=avail -k / | tail -1)" ]; then + work=/mnt/ev + sudo mkdir -p "$work" && sudo chown "$(id -u):$(id -g)" "$work" + fi + echo "EV=${work}" >> "$GITHUB_ENV" + + - name: The GRCh38 reference, its index and dictionary + run: | + set -euo pipefail + . ./versions.env + R="${EV}/genome/reference" + mkdir -p "$R" + f=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna + curl -fsSL --retry 5 -o "${R}/md5.txt" "${NCBI_REF}/md5checksums.txt" + curl -fsSL --retry 5 -o "${R}/${f}.gz" "${NCBI_REF}/${f}.gz" + curl -fsSL --retry 5 -o "${R}/GRCh38_no_alt_analysis_set.fasta.fai" "${NCBI_REF}/${f}.fai" + want=$(awk -v n="./${f}.gz" '$2 == n {print $1}' "${R}/md5.txt") + [ -n "$want" ] && [ "$(md5sum "${R}/${f}.gz" | cut -d' ' -f1)" = "$want" ] || { echo "md5 of ${f}.gz is not NCBI's"; exit 1; } + want=$(awk -v n="./${f}.fai" '$2 == n {print $1}' "${R}/md5.txt") + [ -n "$want" ] && [ "$(md5sum "${R}/GRCh38_no_alt_analysis_set.fasta.fai" | cut -d' ' -f1)" = "$want" ] || { echo "md5 of ${f}.fai is not NCBI's"; exit 1; } + gzip -dc "${R}/${f}.gz" > "${R}/GRCh38_no_alt_analysis_set.fasta" + rm -f "${R}/${f}.gz" + docker run --rm -u "$(id -u):$(id -g)" -v "${R}:/r" "$SAMTOOLS_IMAGE" \ + samtools dict -o /r/GRCh38_no_alt_analysis_set.dict /r/GRCh38_no_alt_analysis_set.fasta + ls -l "$R" + + - uses: actions/download-artifact@9000827ccba6bdab643e8b6fd33ac0654aef8333 # v8.0.2 + with: + name: archaic-1kg-sites + path: sites + + - name: setup.sh --archaic and step 43 on HG002 (chromosome ${{ matrix.part }}) + id: run + continue-on-error: ${{ matrix.part == 'empty' }} + env: + PART: ${{ matrix.part }} + THREADS: '4' + run: | + set -euo pipefail + . ./versions.env + G="${EV}/genome" OUT="${EV}/out/${PART}" + CHROMS=$PART + [ "$PART" != empty ] || CHROMS="" + DATA="${G}/reference/archaic/${ARCHAIC_DATA_VERSION}" + mkdir -p "$OUT" "${DATA}/download" + mv sites/1kg_sites.vcf.gz sites/1kg_sites.vcf.gz.tbi "${DATA}/download/" + # The whole GRCh37 chromosome; the empty job writes an empty list. + : > "${EV}/regions37.bed" + for c in $CHROMS; do printf '%s\t0\t300000000\n' "$c" >> "${EV}/regions37.bed"; done + # The three archaic files side by side, under the names setup.sh's + # fetch resumes (it checks each against the pinned md5 list). + DL0=$(date +%s) + for c in $CHROMS; do + T0=$(date +%s) + for g in Altai Vindija33.19 Denisova; do + curl -fsSL --retry 5 -C - -o "${DATA}/download/${g}.chr${c}.vcf.gz.part" \ + "https://cdna.eva.mpg.de/neandertal/Vindija/VCF/${g}/chr${c}_mq25_mapab100.vcf.gz" & + done + wait + du -cb "${DATA}"/download/*.chr"${c}".vcf.gz.part | tail -1 | cut -f1 > "${EV}/dl_bytes" + echo "chr${c}: $(cat "${EV}/dl_bytes") bytes in $(( $(date +%s) - T0 )) s" + done + T0=$(date +%s) + if ! ARCHAIC_REGIONS="${EV}/regions37.bed" scripts/setup.sh --archaic "$G" 2>&1 | tee "${OUT}/setup.log"; then + echo "setup.sh --archaic failed"; exit 1 + fi + printf 'download_seconds\t%s\ndownload_bytes\t%s\n' "$(( T0 - DL0 ))" \ + "$(cat "${EV}/dl_bytes" 2>/dev/null || echo 0)" > "${OUT}/setup.tsv" + printf 'setup_seconds\t%s\ntable_bytes\t%s\ntable_rows\t%s\n' "$(( $(date +%s) - T0 ))" \ + "$(du -sb "$DATA" | cut -f1)" "$(for t in "${DATA}"/chr*.tsv.gz; do gzip -dc "$t" | tail -n +2; done | wc -l)" \ + >> "${OUT}/setup.tsv" + cat "${OUT}/setup.tsv" + + mkdir -p "${G}/HG002/vcf" + curl -fsSL --retry 5 -o "${G}/HG002/vcf/HG002.vcf.gz" "${GIAB_HG002}.vcf.gz" + curl -fsSL --retry 5 -o "${G}/HG002/vcf/HG002.vcf.gz.tbi" "${GIAB_HG002}.vcf.gz.tbi" + for pop in ALL EUR; do + if ! ARCHAIC_POP=$pop GENOME_DIR="$G" scripts/43-archaic.sh HG002 2>&1 | tee "${OUT}/step43-${pop}.log"; then + echo "step 43 (${pop}) failed"; exit 1 + fi + mkdir -p "${OUT}/${pop}" + cp "${G}/HG002/archaic/work/ibdmix/"chr*.txt "${G}/HG002/archaic/work/ibdmix/spans.tsv" "${OUT}/${pop}/" + sites=$(awk -F'\t' 'NR > 1 {n += $4} END {print n + 0}' "${OUT}/${pop}/spans.tsv") + chroms=$(awk 'NR > 1' "${OUT}/${pop}/spans.tsv" | wc -l) + echo "step 43 (${pop}): ${chroms} chromosomes, ${sites} sites" + [ "$sites" -gt 0 ] || { echo "step 43 (${pop}) analysed no site"; exit 1; } + [ "$chroms" -eq "$(echo "$CHROMS" | wc -w)" ] || { echo "step 43 (${pop}) analysed ${chroms} of the part's chromosomes"; exit 1; } + done + + # The positions of this part where HG002 has two records or more, one + # of them an SNV, as GIAB writes the VCF and split with bcftools + # norm -m-, and how many of them are rows of the table (the rows the + # same-position rule of archaic_report.py decides). + docker run --rm -u "$(id -u):$(id -g)" -v "${G}:/g" "$BCFTOOLS_IMAGE" \ + bcftools norm -m- -Oz -o /g/HG002/split.vcf.gz /g/HG002/vcf/HG002.vcf.gz + python3 -I - "${G}/HG002/vcf/HG002.vcf.gz" "${G}/HG002/split.vcf.gz" "$DATA" "$OUT/affected.tsv" <<'PY' + import glob, gzip, os, sys + raw, split, data, out = sys.argv[1:] + tables = sorted(glob.glob(os.path.join(data, "chr*.tsv.gz"))) + chroms = {os.path.basename(t)[:-7] for t in tables} # this part's chromosomes only + def clashes(path): + """Positions of this part with two records or more, one an SNV (the VCF is sorted).""" + hit, key, n, snvs = set(), None, 0, 0 + for line in gzip.open(path, "rt"): + if line.startswith("#"): + continue + c = line.split("\t", 6) + if c[0] not in chroms: + continue + k = (c[0], int(c[1])) + if k != key: + key, n, snvs = k, 0, 0 + n += 1 + snvs += len(c[3]) == 1 and len(c[4]) == 1 + if n >= 2 and snvs >= 1: + hit.add(key) + return hit + sets = {"raw": clashes(raw), "split": clashes(split)} + on_table = {k: 0 for k in sets} + for t in tables: + chrom = os.path.basename(t)[:-7] + want = {k: {p for c, p in v if c == chrom} for k, v in sets.items()} + if not any(want.values()): + continue + for line in gzip.open(t, "rt"): + c = line.split("\t", 2) + if c[0] != chrom: + continue + p = int(c[1]) + for k in want: + if p in want[k]: + on_table[k] += 1 + with open(out, "w") as f: + for k in sets: + f.write(f"{k}_positions\t{len(sets[k])}\n{k}_on_table\t{on_table[k]}\n") + PY + cat "$OUT/affected.tsv" + + - name: The red control failed, as it must + if: matrix.part == 'empty' + env: + OUTCOME: ${{ steps.run.outcome }} + PART: ${{ matrix.part }} + run: | + echo "outcome of the run on an empty region list: ${OUTCOME}" + [ "$OUTCOME" = failure ] + # and for that reason, not an earlier download or install failure + grep -F 'regions37.bed is not a file' "${EV}/out/${PART}/setup.log" + + - name: Upload the part + if: matrix.part != 'empty' + uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 + with: + name: archaic-part-${{ matrix.part }} + path: ${{ env.EV }}/out/${{ matrix.part }} + if-no-files-found: error + + archaic-share: + needs: archaic-hg002 + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: actions/download-artifact@9000827ccba6bdab643e8b6fd33ac0654aef8333 # v8.0.2 + with: + pattern: archaic-part-* + path: parts + - name: One report per population, over the 22 autosomes + run: | + set -euo pipefail + { + echo "### Step 43 on HG002, whole genome (GIAB v4.2.1 genotypes)" + echo "" + echo "| chromosome | archaic download GB | download seconds | setup seconds | table GB | table rows |" + echo "|---|---|---|---|---|---|" + for p in parts/archaic-part-*; do + kv() { awk -F'\t' -v k="$1" '$1 == k {print $2}' "${p}/setup.tsv"; } + echo "| ${p##*-} | $(awk -v b="$(kv download_bytes)" 'BEGIN {printf "%.2f", b / 1e9}') | $(kv download_seconds) | $(kv setup_seconds) | $(awk -v b="$(kv table_bytes)" 'BEGIN {printf "%.2f", b / 1e9}') | $(kv table_rows) |" + done + cat parts/archaic-part-*/setup.tsv | awk -F'\t' '{t[$1] += $2} + END {printf "| all | %.2f | %d | %d | %.2f | %d |\n", t["download_bytes"] / 1e9, t["download_seconds"], + t["setup_seconds"], t["table_bytes"] / 1e9, t["table_rows"]}' + echo "" + } >> "$GITHUB_STEP_SUMMARY" + for pop in ALL EUR; do + d="joined/${pop}" + mkdir -p "${d}/ibdmix" + cp parts/archaic-part-*/"${pop}"/chr*.txt "${d}/ibdmix/" + { printf 'chrom\tfirst\tlast\tsites\n'; for s in parts/archaic-part-*/"${pop}"/spans.tsv; do tail -n +2 "$s" | cut -f1-4; done; } > "${d}/spans.tsv" + python3 bin/archaic_report.py report --sample HG002 --ibdmix-dir "${d}/ibdmix" --spans "${d}/spans.tsv" \ + --pop "$pop" --out-dir "${d}/out" + kv() { awk -F'\t' -v k="$1" '$1 == k {print $2}' "${d}/out/HG002_archaic.tsv"; } + [ "$(kv analysed_chromosomes)" -eq 22 ] || { echo "${pop}: $(kv analysed_chromosomes) chromosomes, want 22"; exit 1; } + [ "$(kv analysed_sites)" -gt 0 ] || { echo "${pop}: no site analysed"; exit 1; } + nea=$(kv neanderthal_percent) + { + echo "**1000 Genomes ${pop} frequencies:** Neanderthal ${nea}% ($(kv neanderthal_segments) segments, longest $(kv neanderthal_longest_bp) bp), Denisovan $(kv denisovan_percent)% ($(kv denisovan_segments) segments), over $(kv analysed_chromosomes) autosomes, $(kv analysed_bp) bp, $(kv analysed_sites) sites. Compare this share only with other shares from this step, not with published figures (see docs/43-archaic.md)." + echo "" + } >> "$GITHUB_STEP_SUMMARY" + done + { + echo "HG002 positions with two records or more, one an SNV, and how many are table rows:" + echo "" + echo "| VCF | positions | on table rows |" + echo "|---|---|---|" + for k in raw split; do + pos=$(cat parts/archaic-part-*/affected.tsv | awk -F'\t' -v k="${k}_positions" '$1 == k {n += $2} END {print n + 0}') + tab=$(cat parts/archaic-part-*/affected.tsv | awk -F'\t' -v k="${k}_on_table" '$1 == k {n += $2} END {print n + 0}') + case "$k" in raw) what="as GIAB writes it" ;; *) what="split with bcftools norm -m-" ;; esac + echo "| ${what} | ${pos} | ${tab} |" + done + } >> "$GITHUB_STEP_SUMMARY" + cat "$GITHUB_STEP_SUMMARY" + + ancient-footprint: + runs-on: ubuntu-latest + timeout-minutes: 350 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + + - name: Free disk space, add swap and pick the work area + run: | + sudo rm -rf /usr/local/lib/android /usr/share/dotnet /opt/ghc /usr/local/.ghcup \ + /opt/hostedtoolcache/CodeQL /usr/local/share/boost + # Swap keeps the runner alive if the panel's smartpca passes its RAM; + # the job then reports the swap it used. + sudo fallocate -l 8G /mnt/swap-ev && sudo chmod 600 /mnt/swap-ev && sudo mkswap /mnt/swap-ev && sudo swapon /mnt/swap-ev + df -h / /mnt + free -m + nproc + [ -x /usr/bin/time ] || { sudo apt-get update -q && sudo apt-get install -yq time; } + work="${RUNNER_TEMP}/ev" + if [ "$(df --output=avail -k /mnt | tail -1)" -gt "$(df --output=avail -k / | tail -1)" ]; then + work=/mnt/ev + sudo mkdir -p "$work" && sudo chown "$(id -u):$(id -g)" "$work" + fi + echo "EV=${work}" >> "$GITHUB_ENV" + + - name: The GRCh38 reference and its index + run: | + set -euo pipefail + R="${EV}/genome/reference" + mkdir -p "$R" + f=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna + curl -fsSL --retry 5 -o "${R}/md5.txt" "${NCBI_REF}/md5checksums.txt" + curl -fsSL --retry 5 -o "${R}/${f}.gz" "${NCBI_REF}/${f}.gz" + curl -fsSL --retry 5 -o "${R}/GRCh38_no_alt_analysis_set.fasta.fai" "${NCBI_REF}/${f}.fai" + want=$(awk -v n="./${f}.gz" '$2 == n {print $1}' "${R}/md5.txt") + [ -n "$want" ] && [ "$(md5sum "${R}/${f}.gz" | cut -d' ' -f1)" = "$want" ] || { echo "md5 of ${f}.gz is not NCBI's"; exit 1; } + want=$(awk -v n="./${f}.fai" '$2 == n {print $1}' "${R}/md5.txt") + [ -n "$want" ] && [ "$(md5sum "${R}/GRCh38_no_alt_analysis_set.fasta.fai" | cut -d' ' -f1)" = "$want" ] || { echo "md5 of ${f}.fai is not NCBI's"; exit 1; } + gzip -dc "${R}/${f}.gz" > "${R}/GRCh38_no_alt_analysis_set.fasta" + rm -f "${R}/${f}.gz" + + - name: setup.sh --ancient-panel and step 44 on HG002, measured + env: + THREADS: '4' + run: | + set -euo pipefail + . ./versions.env + G="${EV}/genome" M="${EV}/measure" + export GENOME_DIR="$G" + mkdir -p "$M" + # measure NAME: every 2 s, into M/NAME.peak until stopped, the peak of + # any running container's memory.peak (its whole cgroup charge, file + # cache included, which the kernel drops before it kills a task at + # the limit), of the anonymous memory in its memory.stat (what the + # tool itself holds, sampled), of memory.swap.peak, and of the + # growth of used space on the filesystem that holds G (df: other + # writers count, a file made and deleted within 2 s does not); and + # the container's memory.max. + measure() { + local out="${M}/$1.peak" mem=0 anon=0 swap=0 disk=0 limit=0 base v id f + base=$(df --output=used -k "$G" | tail -1) + while :; do + for id in $(docker ps -q --no-trunc); do + for f in "/sys/fs/cgroup/system.slice/docker-${id}.scope" "/sys/fs/cgroup/docker/${id}"; do + [ -r "${f}/memory.peak" ] || continue + v=$(cat "${f}/memory.peak" 2>/dev/null || echo 0) + [ "$v" -le "$mem" ] || { mem=$v; limit=$(cat "${f}/memory.max" 2>/dev/null || echo "$limit"); } + v=$(awk '$1 == "anon" {print $2}' "${f}/memory.stat" 2>/dev/null || true); [ "${v:-0}" -le "$anon" ] || anon=$v + v=$(cat "${f}/memory.swap.peak" 2>/dev/null || echo 0); [ "$v" -le "$swap" ] || swap=$v + done + done + v=$(( $(df --output=used -k "$G" | tail -1) - base )); [ "$v" -le "$disk" ] || disk=$v + printf 'mem_bytes\t%s\nanon_bytes\t%s\nswap_bytes\t%s\ndisk_kb\t%s\nlimit\t%s\n' \ + "$mem" "$anon" "$swap" "$disk" "$limit" > "$out" + sleep 2 + done + } + # timed NAME CMD...: CMD under /usr/bin/time -v with the sampler beside it. + timed() { + local name=$1 rc pid T0 secs + shift + measure "$name" & + pid=$! + T0=$(date +%s) + set +e + /usr/bin/time -v -o "${M}/${name}.time" "$@" 2>&1 | tee "${M}/${name}.log" + rc=${PIPESTATUS[0]} + set -e + secs=$(( $(date +%s) - T0 )) + sleep 3 + kill "$pid" 2>/dev/null || true + wait "$pid" 2>/dev/null || true + printf '%s\t%s\n' seconds "${secs}" >> "${M}/${name}.peak" + cat "${M}/${name}.peak" + grep -E 'Maximum resident|Elapsed' "${M}/${name}.time" || true + [ "$rc" -eq 0 ] || { echo "${name} exited ${rc}"; exit 1; } + [ "$(awk -F'\t' '$1 == "mem_bytes" {print $2}' "${M}/${name}.peak")" -gt 0 ] \ + || { echo "${name}: the sampler saw no container memory; the measurement is void"; exit 1; } + } + + timed setup scripts/setup.sh --ancient-panel "$G" + PANEL="${G}/reference/ancient/aadr_${AADR_VERSION}" + [ -s "${PANEL}/DONE" ] || { echo "no finished panel at ${PANEL}"; exit 1; } + { + printf '%s\t%s\n' panel_bytes "$(du -sb "$PANEL" | cut -f1)" + printf '%s\t%s\n' panel_sites "$(awk -F'\t' '$1 == "sites_grch38" {print $2}' "${PANEL}/DONE")" + printf '%s\t%s\n' moderns "$(wc -l < "${PANEL}/moderns.ind")" + printf '%s\t%s\n' ancients "$(( $(wc -l < "${PANEL}/ancients.tsv") - 1 ))" + } >> "${M}/setup.peak" + + # HG002 as a gVCF at the panel's sites: GIAB's v4.2.1 genotypes, and + # a one-base reference block at every other panel site inside GIAB's + # benchmark regions, where GIAB states HG002 matches the reference. + mkdir -p "${G}/HG002/vcf" + curl -fsSL --retry 5 -o "${EV}/truth.vcf.gz" "${GIAB_HG002}.vcf.gz" + curl -fsSL --retry 5 -o "${EV}/truth.bed" "${GIAB_HG002}_noinconsistent.bed" + python3 -I - "${EV}/truth.vcf.gz" "${EV}/truth.bed" "${G}/reference/GRCh38_no_alt_analysis_set.fasta" \ + "${PANEL}/sites_grch38.tsv" "${EV}/hg002.unsorted.vcf" <<'PY' + import bisect, gzip, sys + truth, bed, ref, sites, out = sys.argv[1:] + fai = {l.split("\t")[0]: [int(x) for x in l.split("\t")[1:5]] for l in open(ref + ".fai")} + fa = open(ref, "rb") + def base(c, p): + n, off, lb, lw = fai[c] + fa.seek(off + ((p - 1) // lb) * lw + (p - 1) % lb) + return fa.read(1).decode().upper() + conf = {} + for l in open(bed): + c, s, e = l.split()[:3] + conf.setdefault(c, []).append((int(s), int(e))) + starts = {c: [s for s, _ in sorted(v)] for c, v in conf.items()} + ends = {c: [e for _, e in sorted(v)] for c, v in conf.items()} + def confident(c, p): + i = bisect.bisect_right(starts.get(c, []), p - 1) - 1 + return i >= 0 and ends[c][i] >= p + covered, head, body = {}, [], [] + for l in gzip.open(truth, "rt"): + if l.startswith("##"): + head.append(l) + continue + if l.startswith("#"): + head.append('##INFO=\n') + head.append('##ALT=\n') + head.append(l) + continue + f = l.split("\t", 5) + p = int(f[1]) + covered.setdefault(f[0], set()).update(range(p, p + len(f[3]))) + body.append(l) + blocks = 0 + with open(sites) as s: + next(s) + for l in s: + _, c, p = l.split("\t")[:3] + p = int(p) + if c not in fai or p in covered.get(c, ()) or not confident(c, p): + continue + b = base(c, p) + if b not in "ACGT": + continue + body.append(f"{c}\t{p}\t.\t{b}\t<*>\t.\t.\tEND={p}\tGT\t0/0\n") + blocks += 1 + with open(out, "w") as o: + o.writelines(head) + o.writelines(body) + print(f"gVCF: {len(body) - blocks} GIAB records, {blocks} reference blocks at panel sites") + PY + docker run --rm -u "$(id -u):$(id -g)" -v "${EV}:/ev" -v "${G}:/g" "$BCFTOOLS_IMAGE" bash -euo pipefail -c ' + bcftools sort -m 2G -T /ev/sort -Oz -o /g/HG002/vcf/HG002.g.vcf.gz /ev/hg002.unsorted.vcf + bcftools index -t /g/HG002/vcf/HG002.g.vcf.gz' + rm -f "${EV}/hg002.unsorted.vcf" + + timed step44 scripts/44-ancient.sh HG002 + INFO="${G}/HG002/ancient/HG002_ancient_info.tsv" + used=$(awk -F'\t' '$1 == "sites_used" {print $2}' "$INFO") + printf '%s\t%s\n' sites_used "${used}" >> "${M}/step44.peak" + [ "${used:-0}" -gt 0 ] || { echo "step 44 used no site"; exit 1; } + + { + echo "### setup.sh --ancient-panel (AADR ${AADR_VERSION}, full 1240K) and step 44 on HG002" + echo "" + echo "Runner: $(nproc) CPUs, $(free -g | awk '/^Mem:/ {print $2}') GB RAM, 8 GB swap added. THREADS=4." + echo "" + echo "| run | wall time | peak cgroup memory (cache included) | that container's limit | peak anonymous memory (sampled) | peak swap | peak growth of used disk (df) | RSS of the host script (/usr/bin/time) |" + echo "|---|---|---|---|---|---|---|---|" + for n in setup step44; do + kv() { awk -F'\t' -v k="$1" '$1 == k {v = $2} END {print v}' "${M}/${n}.peak"; } + rss=$(awk -F': ' '/Maximum resident/ {printf "%.2f GiB", $2 / 1048576}' "${M}/${n}.time") + gib() { awk -v b="$1" 'BEGIN {if (b ~ /^[0-9]+$/) printf "%.2f GiB", b / 1073741824; else print b}'; } + echo "| ${n} | $(( $(kv seconds) / 60 )) min | $(gib "$(kv mem_bytes)") | $(gib "$(kv limit)") | $(gib "$(kv anon_bytes)") | $(gib "$(kv swap_bytes)") | $(awk -v k="$(kv disk_kb)" 'BEGIN {printf "%.1f GiB", k / 1048576}') | ${rss} |" + done + echo "" + kv() { awk -F'\t' -v k="$1" '$1 == k {print $2}' "${M}/setup.peak"; } + echo "Installed panel: $(awk -v b="$(kv panel_bytes)" 'BEGIN {printf "%.2f GiB", b / 1073741824}'), $(kv panel_sites) GRCh38 sites, $(kv moderns) present-day people, $(kv ancients) ancient people. Step 44 used $(awk -F'\t' '$1 == "sites_used" {print $2}' "${M}/step44.peak") sites of HG002." + } >> "$GITHUB_STEP_SUMMARY" + cat "$GITHUB_STEP_SUMMARY" + + - name: Upload the measurements + if: always() + uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 + with: + name: ancient-footprint + path: ${{ env.EV }}/measure + if-no-files-found: ignore diff --git a/bin/archaic_report.py b/bin/archaic_report.py index ad5ce190..fd551702 100755 --- a/bin/archaic_report.py +++ b/bin/archaic_report.py @@ -248,6 +248,34 @@ def read_vcf(path): yield c[0], int(c[1]), c[3], c[4], code +CLASH = object() + + +def pick_call(prev, call): + """The call to keep when a position comes twice, as in a VCF split with + bcftools norm -m-: each is (REF, ALT, code), ALT None for an insertion or a + multi-allelic record. An SNV wins over such a record, in either order. Of + one SNV twice, a usable record wins over code 9, and two different usable + codes are CLASH. Of two different SNVs, the + one that carries its ALT wins over a 0/0 or unusable one (a split 0/1 at + A>G,T is A>G 0/1 and A>T 0/0); both 0/0 is the reference; both carrying + their ALT, or neither readable, is CLASH (unknown, 9).""" + if prev is None: + return call + if prev[1] is None or call[1] is None: + return call if prev[1] is None and call[1] is not None else prev + if call[1] == prev[1]: + if call[2] == 9 or call[2] == prev[2]: + return prev + return call if prev[2] == 9 else CLASH + has_prev, has_call = prev[2] in (1, 2), call[2] in (1, 2) + if has_prev != has_call: + return prev if has_prev else call + if not has_prev and prev[2] == 0 and call[2] == 0: + return (prev[0], None, 0) + return CLASH + + def read_table(path): """Rows of one chromosome of the GRCh38 table, with the header checked.""" with opener(path) as f: @@ -290,9 +318,14 @@ def cmd_genotypes(a): pending = next(vcf, None) continue # The sample's calls on this chromosome, merged with the table by position. - calls = {} + calls, clash = {}, set() while pending is not None and pending[0] == chrom: - calls[pending[1]] = pending[2:] + pos = pending[1] + if pos not in clash: + calls[pos] = pick_call(calls.get(pos), pending[2:]) + if calls[pos] is CLASH: + calls[pos] = (pending[2], None, 9) + clash.add(pos) pending = next(vcf, None) n = first = last = dropped = carried = 0 num = chrom[3:] diff --git a/docs/00-reference-setup.md b/docs/00-reference-setup.md index e1f5487e..2c7d47f4 100644 --- a/docs/00-reference-setup.md +++ b/docs/00-reference-setup.md @@ -402,6 +402,10 @@ The release is `ENSEMBL_GFF_RELEASE` in `versions.env`, and the download is chec Only needed for step 48, the symptom-driven ranking that runs when you give HPO terms. `./scripts/setup.sh --exomiser-data ${GENOME_DIR}` downloads Exomiser's data release `EXOMISER_DATA_VERSION` (hg38 and phenotype archives, 37 GB, each checked against its sha256 in `versions.env`), unpacks it into `${GENOME_DIR}/reference/exomiser-/` (58 GB; about 72 GB free while it runs), deletes the archives and pulls `EXOMISER_IMAGE`. An interrupted download resumes when the command runs again. See [Step 48](48-exomiser.md). +## AADR Ancient Panel (Optional, for Step 44, ~1.05 GB kept, 9.6 GiB at the peak) + +Only needed for step 44, the placement among ancient genomes. `./scripts/setup.sh --ancient-panel ${GENOME_DIR}` downloads the AADR release `AADR_VERSION` (1240K, the 7.1 GB `.geno` and three small files, each checked against its md5 in `versions.env`) and UCSC's hg19ToHg38 chain, then builds the panel in `${GENOME_DIR}/reference/ancient/aadr_/` with convertf and two smartpca runs. On a 4-CPU GitHub runner it took 61 minutes, peaked at 9.6 GiB of extra disk and kept 1.05 GiB; its smartpca container gets 12 GB, and the setup's containers held 7.6 GiB of anonymous memory at most (any setup container, convertf included). See [Step 44](44-ancient.md#runtime). + ## Docker Images — Pre-Pull All `setup.sh` pulls every image in `versions.env` except the lines marked `# optional`, and `scripts/setup.sh --pull-only` pulls the same list without the rest of setup. A failed pull is tried three times (`FETCH_TRIES` and `FETCH_WAIT` change that), and Docker's own message is printed: `toomanyrequests` is Docker Hub's rate limit, `manifest unknown` a tag that does not exist. `--parascopy-data`, `--yleaf-data` and `--megane-kmers` pull their step's optional image along with its data. Each image pinned by digest (`name:tag@sha256:...`) also gets a local tag, `localhost/pgp-pinned/:sha256-`: Docker gives an image pulled by digest no tag, and a dangling-image cleanup such as `docker image prune` would delete it ([Troubleshooting](troubleshooting.md#images-that-were-pulled-are-missing-after-a-docker-cleanup)). Tasks still run the digest reference. To pull every image in advance, the optional ones included, from the repository root in bash (then run `scripts/setup.sh --pull-only` once to tag the digest-pinned ones): diff --git a/docs/43-archaic.md b/docs/43-archaic.md index bb20b16c..2c096a8a 100644 --- a/docs/43-archaic.md +++ b/docs/43-archaic.md @@ -45,9 +45,9 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags. | 1000 Genomes phase 3 | `ALL.wgs.phase3_shapeit2_mvncall_integrated_v5c.20130502.sites.vcf.gz` (1.5 GB) | md5 from the release's `current.tree` | | UCSC | `hg19ToHg38.over.chain.gz` | sha256 | -MPI-EVA publishes md5 lists, not sha256 sums, for 66 files of about 2.5 GB each. Pinning the sha256 of the three lists pins every file's md5 without downloading 165 GB to hash it. +MPI-EVA publishes md5 lists, not sha256 sums, for 66 files of about 2.5 GB each. Pinning the sha256 of the three lists pins every file's md5 without downloading 150 GB to hash it. -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 (about 165 GB in all, under 13 GB at once) 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. +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. ## Input @@ -76,7 +76,12 @@ Both reports show the shares, the analysed span and the ten longest segments in ## Runtime -Not measured on a whole genome. IBDmix runs once per chromosome and archaic genome, `THREADS` chromosomes at a time. The install is bound by the download. +Measured on public HG002 by the `archaic-hg002` jobs of `public-evidence.yml`, one job per autosome on a GitHub runner with 4 CPUs and `THREADS=4` (run 41621716972, 2026-11-03): + +- Install: the archaic files were 150.5 GB in all, from 1.8 GB (chromosome 22) to 13.8 GB (chromosome 2). Downloading them took 2 h 32 min summed over the 22 chromosomes, from 3 to 32 minutes each, at the speed MPI-EVA served that day. Building the table took 3 h 15 min summed, from 2 to 20 minutes per chromosome. The finished table is 0.44 GB with 51,812,908 sites. +- Step 43, one population: 19 to 88 seconds per chromosome, 18 minutes summed over the 22. IBDmix runs once per chromosome and archaic genome, `THREADS` chromosomes at a time. + +The same run gave HG002 a Neanderthal share of 4.13% with the 1000 Genomes ALL frequencies and 4.26% with EUR, over 51,810,955 sites of the 22 autosomes, and a Denisovan share of 0.40% and 0.42%. That is about twice what Chen et al. report with the same method and cutoffs for Europeans (about 51 Mb per person; HG002 got 115.2 Mb here), and the cause is not yet known. Until it is, compare your share with other shares from this step, not with published figures. ## Notes diff --git a/docs/44-ancient.md b/docs/44-ancient.md index 958fc4bd..785913ca 100644 --- a/docs/44-ancient.md +++ b/docs/44-ancient.md @@ -67,8 +67,10 @@ Both reports have an Ancient Populations section: the ten nearest groups, the ca ## Runtime -- Setup: the 7.1 GB download, then convertf and two smartpca runs over the whole release. smartpca holds the packed genotypes of every person in memory, about 7 GB, so the setup container gets 12 GB. -- Step 44: minutes. bcftools reads the gVCF at the panel's sites through its index, and smartpca works on about 820 people plus you; it gets 8 GB. +Measured on public HG002 by the `ancient-footprint` job of `public-evidence.yml`, on a GitHub runner with 4 CPUs, 16 GB of RAM and `THREADS=4` (run 41621716972, 2026-11-03): + +- Setup: 61 minutes, the 7.1 GB download included, then convertf and two smartpca runs over the whole release. smartpca holds the packed genotypes of every person in memory. The highest anonymous memory of any setup container, convertf included, was 7.6 GiB (sampled every 2 seconds; the sampler does not record which container). The container's whole charge, file cache included, reached 11.6 GiB of its 12 GiB limit; the kernel drops that cache before it kills a task, and no swap was used. Used disk grew by 9.6 GiB at most while it ran; the installed panel keeps 1.05 GiB (1,150,106 sites on GRCh38, 819 present-day and 10,580 ancient people). +- Step 44: 7 minutes. bcftools reads the gVCF at the panel's sites through its index, and smartpca works on the 819 people plus you; it gets 8 GB and used 1.2 GiB. The gVCF the job builds from GIAB's HG002 genotypes gave 1,133,951 of the 1,150,106 sites. ## Notes diff --git a/docs/hardware-requirements.md b/docs/hardware-requirements.md index 38c51324..47847550 100644 --- a/docs/hardware-requirements.md +++ b/docs/hardware-requirements.md @@ -68,6 +68,7 @@ While a VEP cache installs, the tarball and the unpacked tree are on disk togeth | CNVpytor GC/mask files | ~0.09 GB | step 18 | | pypgx bundle | ~0.4 GB | step 32 | | Somatic resources (gnomAD AF-only VCF, panel of normals) | ~3 GB | step 29 | +| AADR ancient panel, 1240K (7.1 GB `.geno`; 1.05 GiB kept, 9.6 GiB of extra disk at the peak while it builds) | ~7.1 GB (the `.geno`) plus three small files | step 44, `setup.sh --ancient-panel` | ### Annotation Databases (Optional, for Steps 30-31) @@ -125,10 +126,12 @@ Every step runs in a Docker container with a hard `--memory` limit. The two entr | 25 (PRS): scoring files written as pgsc_calc reads them | no limit | `PRS_PREPARE` 1 / 16 GB | | 25 (PRS): score positions genotyped from the gVCF | 2 / 8 GB | `PRS_SCORE_SITES` 2 / 4 GB | | 25 (PRS): pgsc_calc, each of its tasks | `PGSC_MAX_MEMORY` (three quarters of the RAM) | `PRS` 4 / 24 GB, given to pgsc_calc as `--max_memory` | +| 44 (ancient populations): smartpca | `THREADS` (8) / 8 GB | `ANCIENT_PCA` 4 / 8 GB | +| `setup.sh --ancient-panel`: smartpca over the whole AADR release | `THREADS` (8) / 12 GB | not a pipeline task | | 50 (MEI): MEGAnE | `THREADS` (8) / 16 GB | `MEI` 8 / 16 GB | | 50 (MEI): AnnotSV on the insertions and the calls table | 4 / 8 GB for AnnotSV, 1 / 2 GB for the table | `MEI_REPORT` 4 / 8 GB with `--annotsv_annotations`, 1 / 2 GB without | -These are limits, not measured peaks. The peaks measured for each task in the end test are in the [runtime table](#runtime-per-step): on 30-40x genomes the alignment task (minimap2 and samtools together) peaked at 18.3-18.4 GB and DeepVariant at 15.9-16.5 GB. For step 25, on the nine scores measured before the four small lipid and blood-pressure scores joined the list, `PRS_PREPARE` peaked at 11.5 GB and pgsc_calc's first task, run on one thread, needed about 14 to 18 GB on a GitHub runner without a genome, and 12.7 to 13.5 GiB on the three end-test genomes ([step 25, Memory](25-prs.md#memory)). Your own run records pgsc_calc's task peaks in `prs/pgsc_calc/pipeline_info/execution_trace.txt`; the `PRS` row of the pipeline's trace is only the launcher. +These are limits, not measured peaks. The two step 44 rows were measured on a 4-CPU GitHub runner on public HG002 ([step 44, Runtime](44-ancient.md#runtime)): the setup's containers held 7.6 GiB of anonymous memory at most (the sampler does not record which container, convertf or smartpca), and the largest container's whole charge, file cache included, reached 11.6 of its 12 GiB; step 44's smartpca used 1.2 GiB. The peaks measured for each task in the end test are in the [runtime table](#runtime-per-step): on 30-40x genomes the alignment task (minimap2 and samtools together) peaked at 18.3-18.4 GB and DeepVariant at 15.9-16.5 GB. For step 25, on the nine scores measured before the four small lipid and blood-pressure scores joined the list, `PRS_PREPARE` peaked at 11.5 GB and pgsc_calc's first task, run on one thread, needed about 14 to 18 GB on a GitHub runner without a genome, and 12.7 to 13.5 GiB on the three end-test genomes ([step 25, Memory](25-prs.md#memory)). Your own run records pgsc_calc's task peaks in `prs/pgsc_calc/pipeline_info/execution_trace.txt`; the `PRS` row of the pipeline's trace is only the launcher. **Minimum system RAM:** 16 GB. `run-all.sh` passes the machine's RAM as `--max_memory`, so on a 16 GB machine the 32 GB tasks run capped at about 15 GB. A task that needs more is killed (exit 137); its one retry asks for double, and the cap cuts that back to the same 15 GB. With the step scripts, lower DeepVariant's shards and memory with `THREADS` and `DV_MEM`, for example `THREADS=4 DV_MEM=12g ./scripts/03-deepvariant.sh [male|female]`. The other scripts set a fixed `--memory`; edit it in the script to change it. Step 25's pgsc_calc is the exception that does not fit: its first task needs about 14 to 18 GB on the default scores and was killed under a 14.1 GB limit when the host has no swap, so a 16 GB machine needs swap for it ([step 25, Memory](25-prs.md#memory)). **Recommended:** 32 GB or more, so minimap2 and DeepVariant get close to the 32 GB they ask for. diff --git a/docs/interpreting-results.md b/docs/interpreting-results.md index c0ad6fed..435bc36e 100644 --- a/docs/interpreting-results.md +++ b/docs/interpreting-results.md @@ -657,6 +657,12 @@ One SMN1 copy suggests SMA carrier status; two do not rule it out (two copies on Expect many calls: a genome carries about 1,200 such insertions, nearly all common and harmless. Read the catalog first. **Seen** means an insertion of the same element was called inside the site's window; one copy of the FKTN or BBS1 insertion is carrier status for a recessive condition, and the TAF1 insertion on a male's single X is what causes X-linked dystonia-parkinsonism. **Not seen** is not proof of absence: short reads find an SVA's edges, not its body. **Not callable** means too few reads there (under 10x) to tell. Confirm anything that matters in a clinical laboratory. [Step 50](50-mei.md) explains the columns. +## Ancient Populations (Step 44, opt-in) + +**Where to look:** `${SAMPLE}/ancient/${SAMPLE}_ancient.tsv`, the 20 ancient groups nearest to you on a PCA of present-day West Eurasians, with their dates, places and number of people (`n`), and `${SAMPLE}/ancient/${SAMPLE}_ancient_pca.svg`, the plot with you marked. + +Closeness on the PCA is genetic similarity, not descent: being near an ancient group does not mean you descend from those people. A group with `n` of 1 or 2 has a noisy position. The axes describe West Eurasia only, so a genome with most of its ancestry from elsewhere lands at the edge and its nearest groups say little. `sites_used` in `${SAMPLE}_ancient_info.tsv` tells how many panel sites your gVCF covered; a whole genome covers nearly all of them (1,133,951 of the 1,150,106 that lift to GRCh38, for a gVCF built from the public GIAB HG002 genotypes). [Step 44](44-ancient.md) explains the columns. + ## Polygenic Risk Scores (Step 25) **Where to look:** `${SAMPLE}/prs/${SAMPLE}_prs_summary.tsv`, one row per condition, and the PRS section of the report. diff --git a/docs/nextflow.md b/docs/nextflow.md index 47368f1e..79071c26 100644 --- a/docs/nextflow.md +++ b/docs/nextflow.md @@ -66,6 +66,8 @@ nextflow run main.nf \ # (docs/49-chip.md) # --tools '...,archaic' + --archaic_data (setup.sh --archaic), --archaic_pop, --hla_genes for the # genes of each segment (docs/43-archaic.md) +# --tools '...,ancient' + --ancient_panel (setup.sh --ancient-panel); reads the gVCF, so a sample +# without one is skipped with a warning (docs/44-ancient.md) # --hpo HP:0001250,HP:0001263 + --exomiser_data (setup.sh --exomiser-data): Exomiser ranks the # variants against these symptoms; --exomiser_proband names # the sample when there are several (docs/48-exomiser.md) @@ -216,6 +218,7 @@ results/ │ ├── chip/ # _chip.tsv, _chip_depth.tsv and Mutect2's filtered calls: clonal │ │ # haematopoiesis in 13 genes (opt-in: chip) │ ├── archaic/ # Neanderthal and Denisovan share and segments (IBDmix, opt-in: archaic) +│ ├── ancient/ # The nearest ancient groups (table and text), their info and the PCA plot (AADR, opt-in: ancient) │ ├── hla/ # HLA typing │ ├── expansion_hunter/ # Repeat expansion calls, Stranger's annotation and, when Stranger flags RFC1, its motif count │ ├── telomere/ # Telomere length estimation @@ -325,6 +328,7 @@ The pipeline and the scripts run the same commands at the same tool releases, so | y_haplogroup | Yleaf's prediction table; on the same hardening run | | traits | `_traits.tsv` and `_traits_genotypes.tsv`, byte for byte; in case `traits-panel-1`, which runs step 42 and `--tools traits,html_report` on the same VCF and gVCF. `_genotype_rules.tsv` with WhatsHap's phase from a BAM, byte for byte, in case `genotype-rules-1-phase` | | archaic | `_archaic.tsv`, `_archaic_segments.tsv` and `_archaic.txt`, byte for byte; in case `archaic-ancestry-1`, which runs step 43 and `--tools archaic` on the same VCF and table | +| ancient | the ranked groups of `_ancient.tsv` (`rank`, `group`, `n`) and `sites_used`; in case `ancient-similarity-1`, which runs step 44 and `--tools ancient,html_report` on the same gVCF and panel | A difference fails the check unless the script lists it with its reason; none is listed today. Output file names and folders differ (the scripts write under `$GENOME_DIR//`, the pipeline under `//`, mapped in the script). diff --git a/docs/testing.md b/docs/testing.md index 7b6df129..08c171cd 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -134,11 +134,13 @@ A new check has to be seen failing once before it is trusted: run it against the ### What no e2e case runs -This is every step script that no e2e case runs, plus the full-size index build of step 02a. The first six rows need data, hardware or a binary a GitHub runner does not have. The three long-read rows need long reads, and the e2e fixture has short reads only. The others could run, but no case runs them yet. The last column says what checks each one instead. An image's smoke test is a row of `tests/smoke/commands.tsv` that `container-test.yml` runs on the fixture when the image changes and once a month; it runs the tool, not the step's script. +This is every step script that no e2e case runs, plus the full-size index build of step 02a. The first eight rows need data, hardware or a binary an e2e job does not have. The three long-read rows need long reads, and the e2e fixture has short reads only. The others could run, but no case runs them yet. The last column says what checks each one instead. An image's smoke test is a row of `tests/smoke/commands.tsv` that `container-test.yml` runs on the fixture when the image changes and once a month; it runs the tool, not the step's script. | Step | Why no e2e case runs it | What checks it instead | |---|---|---| | 26b (fine-grained ancestry) on the full panel | the install reads about 19 GB from gnomAD's bucket | case `ancestry-fine-grained-1` runs it on a reduced panel of the fixture's regions plus chr22:20-35 Mb | +| 43 (archaic ancestry) on the whole genome | `setup.sh --archaic` reads 150 GB of archaic genomes | case `archaic-ancestry-1` runs it on a slice; the `archaic-hg002` jobs of `public-evidence.yml` run setup and step 43 on HG002's 22 autosomes, by hand | +| 44 (ancient populations) on the full AADR panel | the install reads the 7.1 GB AADR release and runs smartpca on it for an hour | case `ancient-similarity-1` runs it on a planted panel; the `ancient-footprint` job of `public-evidence.yml` runs setup and step 44 on HG002, by hand | | 13 (VEP, offline) | the cache is 26 GB; the fixture's VEP subset was annotated once with `--database` instead | the image's smoke test, which annotates 50 records with `--database` | | 17 (CPSR) | the PCGR bundle (about 8 GB) plus a second VEP cache. CI never runs real CPSR | the image's smoke test, and the fake-docker cases for the command line | | 05 (AnnotSV) | the annotation data is 5.3 GB | the image's smoke test; the row that annotates the fixture's SV records with the full data runs monthly and on dispatch only | @@ -162,6 +164,13 @@ The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub` [`scripts/ci/settle-doubts.sh`](https://github.com/GeiserX/OwnGenome/blob/main/scripts/ci/settle-doubts.sh) answers questions a review could not settle by reading the code: whether GATK accepts a BAM without a read group, which Clair3 models the image ships, what TIDDIT does with only a BWA-MEM2 index, whether fastp trimming changes DeepVariant's accuracy against the GIAB truth, and so on. It prints one table row per question to the job summary, with the command and what was observed. It is not a gate. Run it by hand with `job: settle-doubts`; it also runs on a pull request, and on main after a merge, that changes the script. +## The public-evidence workflow + +`.github/workflows/public-evidence.yml` runs two opt-in steps whole, on public data only, and writes what they measured to the job summary. It is not a gate; it runs by hand only (`workflow_dispatch`). + +- `ancient-footprint`: `setup.sh --ancient-panel` on the full AADR release, then step 44 on a gVCF built from GIAB's HG002 v4.2.1 genotypes, with the peak memory of each step's containers (one maximum per step, not per container), the peak disk growth and the time. It fails if a command fails, if the sampler saw no container memory, or if step 44 used no site. What it measured is in [step 44, Runtime](44-ancient.md#runtime). +- `archaic-hg002`: `setup.sh --archaic` and step 43 on HG002, one job per autosome; `archaic-share` joins them into one report per population. A job fails if a command fails or step 43 analysed no site, and the job `empty` is its red control: an empty region list must make setup fail. What it measured is in [step 43, Runtime](43-archaic.md#runtime). + ## The prs-memory job It measures the memory of the two PRS tasks that read only the scoring files: `PRS_PREPARE` and pgsc_calc's `FORMAT_SCOREFILES` on the default score list (`tests/test_prs_memory.sh --measure`). One leg runs them under the memory `conf/base.config` gives them and must pass; the other runs them under what the labels gave before and must be killed. It runs only by hand (`job: prs-memory`) and is not a gate. What it measured is in [PRS: Memory](25-prs.md#memory). diff --git a/tests/test_archaic_report.py b/tests/test_archaic_report.py index 1539acf8..f60ebc21 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -11,7 +11,11 @@ 2. genotypes: the stand-in columns carry the population's ALT count; your column is your ALT count, 0 where your VCF has no record, 9 for a non-PASS record; a site where your VCF has another ALT is left out; a - chromosome your VCF has no record on is not analysed; + chromosome your VCF has no record on is not analysed; in a split VCF an + SNV keeps its call beside an insertion, a multi-allelic record or the + 0/0 half of a split site, in either order, a usable record of the same + SNV wins over a filtered one, and two carried SNVs at one position are + unknown (9); 3. report: Altai and Vindija segments merge into Neanderthal ones, a Denisova segment that overlaps one is not counted as Denisovan, segments under 50 kb or LOD 4 are dropped, the share is of the analysed span, and @@ -133,6 +137,49 @@ def image(var): write(f"{TMP}/two.vcf", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tA\tB\n") r = run("genotypes", "--table-dir", tdir, "--vcf", f"{TMP}/two.vcf", "--out-dir", f"{TMP}/gt2") check("a two-sample VCF stops the step", r.returncode != 0 and "one-sample" in r.stderr, r.stderr) +# One position twice, as bcftools norm -m- writes it. Rows: REF>ALT of the table. +sdir2 = f"{TMP}/table_split" +os.makedirs(sdir2) +split_rows = {1000: "A\tG", 2000: "C\tT", 3000: "G\tA", 4000: "T\tC", 5000: "A\tG", 6000: "C\tT", + 7000: "G\tA", 8000: "T\tC", 9000: "A\tC"} +with gzip.open(f"{sdir2}/chr22.tsv.gz", "wt") as f: + f.write(hdr + "".join(f"chr22\t{p}\t{ra}\t0\t0\t0\t0.2\t0\t0\t0\t0\t0\n" for p, ra in split_rows.items())) +write(f"{TMP}/split.vcf", "##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" + "chr22\t1000\t.\tA\tAT\t50\tPASS\t.\tGT\t0/1\n" # insertion, then the SNV + "chr22\t1000\t.\tA\tG\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t2000\t.\tC\tT\t50\tPASS\t.\tGT\t1/1\n" # the SNV, then an insertion + "chr22\t2000\t.\tC\tCA\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t3000\t.\tG\tA,GT\t50\tPASS\t.\tGT\t1/2\n" # a multi-allelic record, then the SNV + "chr22\t3000\t.\tG\tA\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t4000\t.\tT\tC\t50\tPASS\t.\tGT\t0/1\n" # two SNVs, both carried + "chr22\t4000\t.\tT\tG\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t4000\t.\tT\tC\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t5000\t.\tA\tT\t50\tPASS\t.\tGT\t0/0\n" # a split 0/1 at A>G,T: 0/0, then 0/1 + "chr22\t5000\t.\tA\tG\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t6000\t.\tC\tT\t50\tPASS\t.\tGT\t0/1\n" # 0/1, then the 0/0 half + "chr22\t6000\t.\tC\tG\t50\tPASS\t.\tGT\t0/0\n" + "chr22\t7000\t.\tG\tA\t50\tLowQual\t.\tGT\t0/1\n" # the same SNV, filtered, then usable + "chr22\t7000\t.\tG\tA\t50\tPASS\t.\tGT\t1/1\n" + "chr22\t8000\t.\tT\tA\t50\tPASS\t.\tGT\t0/0\n" # two 0/0 halves: the reference + "chr22\t8000\t.\tT\tG\t50\tPASS\t.\tGT\t0/0\n") +r = run("genotypes", "--table-dir", sdir2, "--vcf", f"{TMP}/split.vcf", "--pseudo", "10", "--out-dir", f"{TMP}/gt3") +check("genotypes on a split VCF exits 0", r.returncode == 0, r.stderr) +g3 = {int(l.split("\t")[1]): l.split("\t") for l in gzip.open(f"{TMP}/gt3/chr22.gt.gz", "rt").read().splitlines()[1:]} +code3 = lambda p: g3.get(p, ["left out"] * 8)[7] +check("an insertion before the SNV at one position keeps the SNV call", code3(1000) == "1", g3.get(1000)) +check("an insertion after the SNV at one position keeps the SNV call", code3(2000) == "2", g3.get(2000)) +check("a multi-allelic record beside the SNV keeps the SNV call", code3(3000) == "1", g3.get(3000)) +check("two carried SNVs at one position are unknown (9), also when a third record follows", code3(4000) == "9", g3.get(4000)) +check("a 0/0 half of a split site before the carried SNV keeps the SNV call", code3(5000) == "1", g3.get(5000)) +check("a 0/0 half after the carried SNV keeps the SNV call", code3(6000) == "1", g3.get(6000)) +check("a usable record of the same SNV wins over a filtered one", code3(7000) == "2", g3.get(7000)) +check("two 0/0 halves are the reference (0)", code3(8000) == "0", g3.get(8000)) +pc_code = lambda r: 9 if r is ar.CLASH else r[2] +check("pick_call gives the same code in either order", + all(pc_code(ar.pick_call(x, y)) == pc_code(ar.pick_call(y, x)) + for x, y in [(("A", "G", 1), ("A", None, 9)), (("A", "G", 1), ("A", "T", 0)), (("A", "G", 9), ("A", "G", 2)), + (("A", "G", 0), ("A", "T", 0)), (("A", "G", 1), ("A", "T", 2)), (("A", "G", 9), ("A", "T", 0)), + (("A", "G", 1), ("A", "G", 2))])) # --- 3. report --------------------------------------------------------------------- idir = f"{TMP}/ibd"