From f49f87c3b383ea40d52c35ce54973aac4ef1e75a Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 10:01:59 +0100 Subject: [PATCH 01/10] fix(archaic): keep the SNV call when an insertion sits at the same position A VCF split with bcftools norm -m- can hold an SNV and an insertion (or a multi-allelic record) at one position. Step 43 kept whichever came last, so an insertion after the SNV turned the call into "unknown". The SNV now wins in either order, and two different SNVs at one position are unknown. Bead: pgp-bgr1 --- bin/archaic_report.py | 15 +++++++++++++-- tests/test_archaic_report.py | 24 +++++++++++++++++++++++- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/bin/archaic_report.py b/bin/archaic_report.py index ad5ce190..2044fa77 100755 --- a/bin/archaic_report.py +++ b/bin/archaic_report.py @@ -290,9 +290,20 @@ def cmd_genotypes(a): pending = next(vcf, None) continue # The sample's calls on this chromosome, merged with the table by position. - calls = {} + # A position can come twice, as in a VCF split with bcftools norm -m-: + # an SNV beside an insertion or a multi-allelic record keeps the SNV, + # whichever comes first, and two different SNVs make it unknown (9). + calls, clash = {}, set() while pending is not None and pending[0] == chrom: - calls[pending[1]] = pending[2:] + pos, call = pending[1], pending[2:] + prev = calls.get(pos) + if pos in clash: + pass + elif prev is None or (prev[1] is None and call[1] is not None): + calls[pos] = call + elif call[1] is not None and call[1] != prev[1]: + calls[pos] = (prev[0], None, 9) + clash.add(pos) pending = next(vcf, None) n = first = last = dropped = carried = 0 num = chrom[3:] diff --git a/tests/test_archaic_report.py b/tests/test_archaic_report.py index 1539acf8..d76b0054 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -11,7 +11,9 @@ 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 or a multi-allelic record at the + same position, and two different SNVs there 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 +135,26 @@ 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: the SNV is kept beside an +# insertion or a multi-allelic record, in either order; two SNVs are unknown. +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" + "chr22\t1000\t.\tA\tG\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t2000\t.\tC\tT\t50\tPASS\t.\tGT\t1/1\n" + "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" + "chr22\t3000\t.\tG\tA\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t4000\t.\tT\tC\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t4000\t.\tT\tG\t50\tPASS\t.\tGT\t0/1\n" + "chr22\t4000\t.\tT\tC\t50\tPASS\t.\tGT\t0/1\n") +r = run("genotypes", "--table-dir", tdir, "--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:]} +check("an insertion before the SNV at one position keeps the SNV call", g3.get(1000, [""] * 8)[7] == "1", g3.get(1000)) +check("an insertion after the SNV at one position keeps the SNV call", g3.get(2000, [""] * 8)[7] == "2", g3.get(2000)) +check("a multi-allelic record beside the SNV keeps the SNV call", g3.get(3000, [""] * 8)[7] == "1", g3.get(3000)) +check("two different SNVs at one position are unknown (9), also when a third record follows", + g3.get(4000, [""] * 8)[7] == "9", g3.get(4000)) # --- 3. report --------------------------------------------------------------------- idir = f"{TMP}/ibd" From 74085939e37b375ce035247de9412562cd67bcc1 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 10:06:40 +0100 Subject: [PATCH 02/10] ci: run steps 43 and 44 whole on public HG002 and measure what they need A new dispatch workflow (also run on a branch push that changes it or bin/archaic_report.py) runs on GitHub-hosted runners and public data only. archaic-hg002 builds the archaic table in four parts and runs step 43 on GIAB's HG002 genotypes; archaic-share joins the parts and prints HG002's shares and the positions the SNV-over-insertion fix touches. An empty region list is the red control and must fail. ancient-footprint installs the full AADR panel and runs step 44, with the peak container memory, disk and time. Bead: pgp-bgr1 Bead: pgp-sejl --- .github/workflows/public-evidence.yml | 447 ++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 .github/workflows/public-evidence.yml diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml new file mode 100644 index 00000000..94f150a0 --- /dev/null +++ b/.github/workflows/public-evidence.yml @@ -0,0 +1,447 @@ +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-hg002 setup.sh --archaic and step 43 on HG002's GIAB v4.2.1 +# genotypes, the 22 autosomes in four parts (each part +# builds its own table: about 165 GB of downloads in +# all), then archaic-share joins the parts into one +# report per population and prints HG002's Neanderthal +# and Denisovan shares, and how many of HG002's +# positions carry an SNV and another record at once. +# The part "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. +# +# A dispatch-only workflow cannot run before it is on main, so a push that +# changes this file or bin/archaic_report.py runs it too (not on main, where a +# dispatch reruns it). Public repo: the GitHub-hosted runner is free here. +on: + push: + branches-ignore: [main] + paths: + - '.github/workflows/public-evidence.yml' + - 'bin/archaic_report.py' + 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-hg002: + name: archaic-hg002 (${{ matrix.part }}) + runs-on: ubuntu-latest + timeout-minutes: 350 + strategy: + fail-fast: false + matrix: + include: + # About 700 Mb of GRCh37 each, so the downloads split evenly. + - {part: a, chroms: '1 2 3', expect: pass} + - {part: b, chroms: '4 5 6 7', expect: pass} + - {part: c, chroms: '8 9 10 11 12', expect: pass} + - {part: d, chroms: '13 14 15 16 17 18 19 20 21 22', expect: pass} + - {part: empty, chroms: '', expect: fail} + 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" + + - name: setup.sh --archaic and step 43 on HG002 (chromosomes ${{ matrix.chroms || 'none' }}) + id: run + continue-on-error: ${{ matrix.expect == 'fail' }} + env: + CHROMS: ${{ matrix.chroms }} + PART: ${{ matrix.part }} + THREADS: '4' + run: | + set -euo pipefail + . ./versions.env + G="${EV}/genome" OUT="${EV}/out/${PART}" + mkdir -p "$OUT" + # Whole GRCh37 chromosomes; the empty part writes an empty list. + : > "${EV}/regions37.bed" + for c in $CHROMS; do printf '%s\t0\t300000000\n' "$c" >> "${EV}/regions37.bed"; done + T0=$(date +%s) + ARCHAIC_REGIONS="${EV}/regions37.bed" scripts/setup.sh --archaic "$G" 2>&1 | tee "${OUT}/setup.log" + rc=${PIPESTATUS[0]} + [ "$rc" -eq 0 ] || { echo "setup.sh --archaic exited ${rc}"; exit 1; } + DATA="${G}/reference/archaic/${ARCHAIC_DATA_VERSION}" + 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 + ARCHAIC_POP=$pop GENOME_DIR="$G" scripts/43-archaic.sh HG002 2>&1 | tee "${OUT}/step43-${pop}.log" + rc=${PIPESTATUS[0]} + [ "$rc" -eq 0 ] || { echo "step 43 (${pop}) exited ${rc}"; exit 1; } + 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 where HG002 carries an SNV and another record, as + # GIAB writes the VCF and split with bcftools norm -m-, and how + # many of them are rows of the table (the rows the fix changes). + 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:] + def clashes(path): + """Positions with an SNV record and another record (the VCF is sorted).""" + hit, key, kinds = set(), None, set() + for line in gzip.open(path, "rt"): + if line.startswith("#"): + continue + c = line.split("\t", 6) + k = (c[0], int(c[1])) + if k != key: + key, kinds = k, set() + kinds.add(len(c[3]) == 1 and len(c[4]) == 1) + if len(kinds) == 2: + hit.add(key) + return hit + sets = {"raw": clashes(raw), "split": clashes(split)} + on_table = {k: 0 for k in sets} + for t in glob.glob(os.path.join(data, "chr*.tsv.gz")): + 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.expect == 'fail' + env: + OUTCOME: ${{ steps.run.outcome }} + run: | + echo "outcome of the run on an empty region list: ${OUTCOME}" + [ "$OUTCOME" = failure ] + + - name: Upload the part + if: matrix.expect == 'pass' + 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 "| part | 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##*-} | $(kv setup_seconds) | $(awk -v b="$(kv table_bytes)" 'BEGIN {printf "%.2f", b / 1e9}') | $(kv table_rows) |" + done + 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) + band=$(awk -v x="$nea" 'BEGIN {print (x >= 1 && x <= 2) ? "inside" : "outside"}') + { + 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. The Neanderthal share is ${band} the published 1-2% for people outside Africa." + echo "" + } >> "$GITHUB_STEP_SUMMARY" + done + { + echo "HG002 positions with an SNV and another record at once, 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: peak memory.peak and memory.swap.peak of any running + # container, and peak growth of the disk under G, every 2 s, into + # M/NAME.peak, until stopped. + measure() { + local out="${M}/$1.peak" mem=0 swap=0 disk=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"); [ "$v" -le "$mem" ] || mem=$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\nswap_bytes\t%s\ndisk_kb\t%s\n' "$mem" "$swap" "$disk" > "$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 container memory | peak container swap | peak disk growth | setup.sh or step RSS (/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 GB", $2 / 1048576}' "${M}/${n}.time") + echo "| ${n} | $(( $(kv seconds) / 60 )) min | $(awk -v b="$(kv mem_bytes)" 'BEGIN {printf "%.2f GB", b / 1073741824}') | $(awk -v b="$(kv swap_bytes)" 'BEGIN {printf "%.2f GB", b / 1073741824}') | $(awk -v k="$(kv disk_kb)" 'BEGIN {printf "%.1f GB", 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 GB", b / 1e9}'), $(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 From 0f623c0b90fbabc9497632a2442d430b9962325b Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 10:29:33 +0100 Subject: [PATCH 03/10] fix(archaic): keep the carried SNV beside the 0/0 half of a split site bcftools norm -m- writes a 0/1 at A>G,T as A>G 0/1 and A>T 0/0. The same-position rule took that as two different SNVs and made the site unknown. pick_call now keeps the record that carries its ALT, lets a usable record of the same SNV win over a filtered one, and reads two 0/0 halves as the reference, in either order. Bead: pgp-bgr1 --- bin/archaic_report.py | 43 ++++++++++++++++++++-------- tests/test_archaic_report.py | 54 ++++++++++++++++++++++++++---------- 2 files changed, 70 insertions(+), 27 deletions(-) diff --git a/bin/archaic_report.py b/bin/archaic_report.py index 2044fa77..33553f1a 100755 --- a/bin/archaic_report.py +++ b/bin/archaic_report.py @@ -248,6 +248,31 @@ 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. 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]: + return call if prev[2] == 9 and call[2] != 9 else prev + 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,20 +315,14 @@ def cmd_genotypes(a): pending = next(vcf, None) continue # The sample's calls on this chromosome, merged with the table by position. - # A position can come twice, as in a VCF split with bcftools norm -m-: - # an SNV beside an insertion or a multi-allelic record keeps the SNV, - # whichever comes first, and two different SNVs make it unknown (9). calls, clash = {}, set() while pending is not None and pending[0] == chrom: - pos, call = pending[1], pending[2:] - prev = calls.get(pos) - if pos in clash: - pass - elif prev is None or (prev[1] is None and call[1] is not None): - calls[pos] = call - elif call[1] is not None and call[1] != prev[1]: - calls[pos] = (prev[0], None, 9) - clash.add(pos) + 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/tests/test_archaic_report.py b/tests/test_archaic_report.py index d76b0054..c1b76b31 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -12,8 +12,10 @@ 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; in a split VCF an - SNV keeps its call beside an insertion or a multi-allelic record at the - same position, and two different SNVs there are unknown (9); + 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 @@ -135,26 +137,48 @@ 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: the SNV is kept beside an -# insertion or a multi-allelic record, in either order; two SNVs are unknown. +# 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" + "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" + "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" + "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" + "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") -r = run("genotypes", "--table-dir", tdir, "--vcf", f"{TMP}/split.vcf", "--pseudo", "10", "--out-dir", f"{TMP}/gt3") + "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:]} -check("an insertion before the SNV at one position keeps the SNV call", g3.get(1000, [""] * 8)[7] == "1", g3.get(1000)) -check("an insertion after the SNV at one position keeps the SNV call", g3.get(2000, [""] * 8)[7] == "2", g3.get(2000)) -check("a multi-allelic record beside the SNV keeps the SNV call", g3.get(3000, [""] * 8)[7] == "1", g3.get(3000)) -check("two different SNVs at one position are unknown (9), also when a third record follows", - g3.get(4000, [""] * 8)[7] == "9", g3.get(4000)) +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))])) # --- 3. report --------------------------------------------------------------------- idir = f"{TMP}/ibd" From b27a9ed200dc7092f5135a29dc1068b67ca185c0 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 10:29:36 +0100 Subject: [PATCH 04/10] ci(public-evidence): count each part's own positions and check why the control fails Each part counted the same-position records of the whole VCF, so the joined total was four times too high, and two SNVs at one position were not counted. The red control now also requires the empty-region message, and the disk figure is named for what df measures. Bead: pgp-bgr1 --- .github/workflows/public-evidence.yml | 47 +++++++++++++++++---------- 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 94f150a0..62e1f2d0 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -108,9 +108,9 @@ jobs: : > "${EV}/regions37.bed" for c in $CHROMS; do printf '%s\t0\t300000000\n' "$c" >> "${EV}/regions37.bed"; done T0=$(date +%s) - ARCHAIC_REGIONS="${EV}/regions37.bed" scripts/setup.sh --archaic "$G" 2>&1 | tee "${OUT}/setup.log" - rc=${PIPESTATUS[0]} - [ "$rc" -eq 0 ] || { echo "setup.sh --archaic exited ${rc}"; exit 1; } + 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 DATA="${G}/reference/archaic/${ARCHAIC_DATA_VERSION}" 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)" \ @@ -121,9 +121,9 @@ jobs: 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 - ARCHAIC_POP=$pop GENOME_DIR="$G" scripts/43-archaic.sh HG002 2>&1 | tee "${OUT}/step43-${pop}.log" - rc=${PIPESTATUS[0]} - [ "$rc" -eq 0 ] || { echo "step 43 (${pop}) exited ${rc}"; exit 1; } + 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") @@ -133,31 +133,37 @@ jobs: [ "$chroms" -eq "$(echo "$CHROMS" | wc -w)" ] || { echo "step 43 (${pop}) analysed ${chroms} of the part's chromosomes"; exit 1; } done - # The positions where HG002 carries an SNV and another record, as - # GIAB writes the VCF and split with bcftools norm -m-, and how - # many of them are rows of the table (the rows the fix changes). + # 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 with an SNV record and another record (the VCF is sorted).""" - hit, key, kinds = set(), None, set() + """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, kinds = k, set() - kinds.add(len(c[3]) == 1 and len(c[4]) == 1) - if len(kinds) == 2: + 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 glob.glob(os.path.join(data, "chr*.tsv.gz")): + 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()): @@ -180,9 +186,12 @@ jobs: if: matrix.expect == 'fail' 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.expect == 'pass' @@ -236,7 +245,7 @@ jobs: } >> "$GITHUB_STEP_SUMMARY" done { - echo "HG002 positions with an SNV and another record at once, and how many are table rows:" + 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 "|---|---|---|" @@ -301,7 +310,9 @@ jobs: export GENOME_DIR="$G" mkdir -p "$M" # measure NAME: peak memory.peak and memory.swap.peak of any running - # container, and peak growth of the disk under G, every 2 s, into + # container, and peak growth of the used space of the filesystem that + # holds G (df: other writers count, a file made and deleted within + # 2 s does not), every 2 s, into # M/NAME.peak, until stopped. measure() { local out="${M}/$1.peak" mem=0 swap=0 disk=0 base v id f @@ -425,7 +436,7 @@ jobs: 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 container memory | peak container swap | peak disk growth | setup.sh or step RSS (/usr/bin/time) |" + echo "| run | wall time | peak container memory | peak container swap | peak growth of used disk (df) | setup.sh or step RSS (/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"; } From 9f71e527a7456b7ca67b8462729b9f773704c2ca Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 10:36:04 +0100 Subject: [PATCH 05/10] ci(public-evidence): one job per autosome, with the archaic files fetched side by side The archaic files come from one server at a few MB/s per connection, so four parts of about 40 GB each would not finish inside a job's six hours. Each autosome now runs in its own job (eight at a time), fetches its three files in parallel under the names setup.sh resumes, and takes the 1000 Genomes sites file from one shared job instead of downloading it 22 times. Bead: pgp-bgr1 --- .github/workflows/public-evidence.yml | 104 +++++++++++++++++++------- 1 file changed, 78 insertions(+), 26 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 62e1f2d0..13fae44c 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -4,14 +4,18 @@ name: Public evidence # 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, the 22 autosomes in four parts (each part -# builds its own table: about 165 GB of downloads in -# all), then archaic-share joins the parts into one -# report per population and prints HG002's Neanderthal -# and Denisovan shares, and how many of HG002's -# positions carry an SNV and another record at once. -# The part "empty" is the red control: an empty region +# genotypes, one job per autosome (about 165 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 @@ -41,20 +45,43 @@ env: 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: - include: - # About 700 Mb of GRCh37 each, so the downloads split evenly. - - {part: a, chroms: '1 2 3', expect: pass} - - {part: b, chroms: '4 5 6 7', expect: pass} - - {part: c, chroms: '8 9 10 11 12', expect: pass} - - {part: d, chroms: '13 14 15 16 17 18 19 20 21 22', expect: pass} - - {part: empty, chroms: '', expect: fail} + 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: @@ -92,29 +119,51 @@ jobs: samtools dict -o /r/GRCh38_no_alt_analysis_set.dict /r/GRCh38_no_alt_analysis_set.fasta ls -l "$R" - - name: setup.sh --archaic and step 43 on HG002 (chromosomes ${{ matrix.chroms || 'none' }}) + - 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.expect == 'fail' }} + continue-on-error: ${{ matrix.part == 'empty' }} env: - CHROMS: ${{ matrix.chroms }} PART: ${{ matrix.part }} THREADS: '4' run: | set -euo pipefail . ./versions.env G="${EV}/genome" OUT="${EV}/out/${PART}" - mkdir -p "$OUT" - # Whole GRCh37 chromosomes; the empty part writes an empty list. + 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 - DATA="${G}/reference/archaic/${ARCHAIC_DATA_VERSION}" + 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" + >> "${OUT}/setup.tsv" cat "${OUT}/setup.tsv" mkdir -p "${G}/HG002/vcf" @@ -183,7 +232,7 @@ jobs: cat "$OUT/affected.tsv" - name: The red control failed, as it must - if: matrix.expect == 'fail' + if: matrix.part == 'empty' env: OUTCOME: ${{ steps.run.outcome }} PART: ${{ matrix.part }} @@ -194,7 +243,7 @@ jobs: grep -F 'regions37.bed is not a file' "${EV}/out/${PART}/setup.log" - name: Upload the part - if: matrix.expect == 'pass' + if: matrix.part != 'empty' uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 with: name: archaic-part-${{ matrix.part }} @@ -219,12 +268,15 @@ jobs: { echo "### Step 43 on HG002, whole genome (GIAB v4.2.1 genotypes)" echo "" - echo "| part | setup seconds | table GB | table rows |" - 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##*-} | $(kv setup_seconds) | $(awk -v b="$(kv table_bytes)" 'BEGIN {printf "%.2f", b / 1e9}') | $(kv table_rows) |" + 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 From 2f91b41f3f971a97c1885122cb4f3cf16d0d46da Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 11:44:00 +0100 Subject: [PATCH 06/10] ci(public-evidence): report the anonymous memory beside the cgroup peak The first full run put the panel install's cgroup peak at 11.7 GiB of its 12 GiB limit. memory.peak counts the file cache of the 7.1 GB genotype file, which the kernel drops before it kills a task, so it overstates what smartpca holds. The sampler now also records the anonymous memory from memory.stat and the container's memory.max. Bead: pgp-sejl --- .github/workflows/public-evidence.yml | 32 ++++++++++++++++----------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 13fae44c..2c8f4c20 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -361,24 +361,29 @@ jobs: G="${EV}/genome" M="${EV}/measure" export GENOME_DIR="$G" mkdir -p "$M" - # measure NAME: peak memory.peak and memory.swap.peak of any running - # container, and peak growth of the used space of the filesystem that - # holds G (df: other writers count, a file made and deleted within - # 2 s does not), every 2 s, into - # M/NAME.peak, until stopped. + # 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 swap=0 disk=0 base v id f + 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"); [ "$v" -le "$mem" ] || mem=$v + v=$(cat "${f}/memory.peak"); [ "$v" -le "$mem" ] || { mem=$v; limit=$(cat "${f}/memory.max"); } + v=$(awk '$1 == "anon" {print $2}' "${f}/memory.stat"); [ "${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\nswap_bytes\t%s\ndisk_kb\t%s\n' "$mem" "$swap" "$disk" > "$out" + 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 } @@ -488,16 +493,17 @@ jobs: 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 container memory | peak container swap | peak growth of used disk (df) | setup.sh or step RSS (/usr/bin/time) |" - 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 GB", $2 / 1048576}' "${M}/${n}.time") - echo "| ${n} | $(( $(kv seconds) / 60 )) min | $(awk -v b="$(kv mem_bytes)" 'BEGIN {printf "%.2f GB", b / 1073741824}') | $(awk -v b="$(kv swap_bytes)" 'BEGIN {printf "%.2f GB", b / 1073741824}') | $(awk -v k="$(kv disk_kb)" 'BEGIN {printf "%.1f GB", k / 1048576}') | ${rss} |" + 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 GB", b / 1e9}'), $(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." + 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" From a6e77528db82c77bb872d5b04a0545e035035e48 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 13:31:40 +0100 Subject: [PATCH 07/10] docs(ancient): add step 44 to the shared docs, with its measured memory, disk and time The AADR install and step 44 ran whole on public HG002 in the ancient-footprint job of public-evidence.yml (run 41621716972), so the docs give what that run measured in place of the estimates. Bead: pgp-sejl --- docs/00-reference-setup.md | 4 ++++ docs/44-ancient.md | 6 ++++-- docs/hardware-requirements.md | 5 ++++- docs/interpreting-results.md | 6 ++++++ docs/nextflow.md | 4 ++++ docs/testing.md | 11 ++++++++++- 6 files changed, 32 insertions(+), 4 deletions(-) diff --git a/docs/00-reference-setup.md b/docs/00-reference-setup.md index a8cd22fe..2f61f421 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, ~7.1 GB) + +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, of which smartpca itself used 7.6 GiB. 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/44-ancient.md b/docs/44-ancient.md index 958fc4bd..51de8109 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: its own memory peaked at 7.6 GiB (anonymous memory, sampled every 2 seconds). 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 502fb4f5..5567de28 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 | step 44, `setup.sh --ancient-panel` | ### Annotation Databases (Optional, for Steps 30-31) @@ -125,8 +126,10 @@ 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 | -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 panel's smartpca held 7.6 GiB of its own memory, and its 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 af93d33c..e114ea75 100644 --- a/docs/interpreting-results.md +++ b/docs/interpreting-results.md @@ -656,6 +656,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 d67aa8a8..a28a8a6f 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 081508b5..88ee2035 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -132,11 +132,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 about 165 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 | @@ -160,6 +162,13 @@ The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub` [`scripts/ci/settle-doubts.sh`](https://github.com/GeiserX/Personal-Genome-Pipeline/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; run it by hand (`workflow_dispatch`), and a push to a branch other than main that changes the workflow or `bin/archaic_report.py` runs it too. + +- `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 each container's peak memory, 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). From b013ebff721205ae75c73b7a8ec7c11c6afec024 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 16:01:21 +0100 Subject: [PATCH 08/10] docs(archaic): give step 43's measured install, run time and HG002 share The whole-genome run of step 43 on public HG002 (run 41621716972) downloaded 150 GB, not 165, and chromosome 2 alone is 13.8 GB. HG002's Neanderthal share came out at 4.1 to 4.3%, about twice the quoted 2%: the step counts a stretch on either copy, the quoted figure is of one copy. The page now says so. Bead: pgp-bgr1 --- .github/workflows/public-evidence.yml | 2 +- docs/43-archaic.md | 11 ++++++++--- docs/testing.md | 2 +- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 2c8f4c20..9d7c735c 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -8,7 +8,7 @@ name: Public evidence # downloaded once (its md5 checked) and indexed, for # every chromosome job. # archaic-hg002 setup.sh --archaic and step 43 on HG002's GIAB v4.2.1 -# genotypes, one job per autosome (about 165 GB of +# 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 diff --git a/docs/43-archaic.md b/docs/43-archaic.md index bb20b16c..6ade3096 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 the 2% quoted above, and the two are not the same measure: the published share is of one copy of the genome, while this step counts a stretch when either of your two copies carries it, which for rare stretches is close to the sum of the two copies. Compare your share with other people's shares from this step, not with the 2%. ## Notes diff --git a/docs/testing.md b/docs/testing.md index 88ee2035..645288bd 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -137,7 +137,7 @@ This is every step script that no e2e case runs, plus the full-size index build | Step | Why no e2e case runs it | What checks it instead | |---|---|---| | 26b (fine-grained ancestry) on the full panel | the install reads about 19 GB from gnomAD's bucket | case `ancestry-fine-grained-1` runs it on a reduced panel of the fixture's regions plus chr22:20-35 Mb | -| 43 (archaic ancestry) on the whole genome | `setup.sh --archaic` reads about 165 GB of archaic genomes | case `archaic-ancestry-1` runs it on a slice; the `archaic-hg002` jobs of `public-evidence.yml` run setup and step 43 on HG002's 22 autosomes, by hand | +| 43 (archaic ancestry) on the whole genome | `setup.sh --archaic` reads 150 GB of archaic genomes | case `archaic-ancestry-1` runs it on 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 | From da91f7adcee3c78bc1a358a7bd395f57d1fc3b1c Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 17:23:39 +0100 Subject: [PATCH 09/10] fix(archaic): settle two usable codes of one SNV as unknown, and keep the docs to what was measured The archaic share gap against Chen et al. is stated as not yet explained. public-evidence runs by hand only, and its sampler no longer stops quietly when a container exits between reads. The panel's disk and memory figures say what the sampler measured. --- .github/workflows/public-evidence.yml | 15 +++++---------- bin/archaic_report.py | 7 +++++-- docs/00-reference-setup.md | 4 ++-- docs/43-archaic.md | 2 +- docs/44-ancient.md | 2 +- docs/testing.md | 2 +- tests/test_archaic_report.py | 3 ++- 7 files changed, 17 insertions(+), 18 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 9d7c735c..01c02b9d 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -22,15 +22,9 @@ name: Public evidence # its containers (the cgroup's memory.peak), the peak # disk use and the time. # -# A dispatch-only workflow cannot run before it is on main, so a push that -# changes this file or bin/archaic_report.py runs it too (not on main, where a -# dispatch reruns it). Public repo: the GitHub-hosted runner is free here. +# It runs by hand only (workflow_dispatch). Public repo: the GitHub-hosted +# runner is free here. on: - push: - branches-ignore: [main] - paths: - - '.github/workflows/public-evidence.yml' - - 'bin/archaic_report.py' workflow_dispatch: permissions: @@ -376,8 +370,9 @@ jobs: 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"); [ "$v" -le "$mem" ] || { mem=$v; limit=$(cat "${f}/memory.max"); } - v=$(awk '$1 == "anon" {print $2}' "${f}/memory.stat"); [ "${v:-0}" -le "$anon" ] || anon=$v + 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 diff --git a/bin/archaic_report.py b/bin/archaic_report.py index 33553f1a..fd551702 100755 --- a/bin/archaic_report.py +++ b/bin/archaic_report.py @@ -255,7 +255,8 @@ 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. Of two different SNVs, the + 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).""" @@ -264,7 +265,9 @@ def pick_call(prev, 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]: - return call if prev[2] == 9 and call[2] != 9 else prev + 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 diff --git a/docs/00-reference-setup.md b/docs/00-reference-setup.md index 0821a2ba..2c7d47f4 100644 --- a/docs/00-reference-setup.md +++ b/docs/00-reference-setup.md @@ -402,9 +402,9 @@ 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, ~7.1 GB) +## 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, of which smartpca itself used 7.6 GiB. See [Step 44](44-ancient.md#runtime). +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 diff --git a/docs/43-archaic.md b/docs/43-archaic.md index 6ade3096..2c096a8a 100644 --- a/docs/43-archaic.md +++ b/docs/43-archaic.md @@ -81,7 +81,7 @@ Measured on public HG002 by the `archaic-hg002` jobs of `public-evidence.yml`, o - 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 the 2% quoted above, and the two are not the same measure: the published share is of one copy of the genome, while this step counts a stretch when either of your two copies carries it, which for rare stretches is close to the sum of the two copies. Compare your share with other people's shares from this step, not with the 2%. +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 51de8109..785913ca 100644 --- a/docs/44-ancient.md +++ b/docs/44-ancient.md @@ -69,7 +69,7 @@ Both reports have an Ancient Populations section: the ten nearest groups, the ca 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: its own memory peaked at 7.6 GiB (anonymous memory, sampled every 2 seconds). 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). +- 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/testing.md b/docs/testing.md index 4ecc962c..e098f575 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -166,7 +166,7 @@ The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub` ## The public-evidence workflow -`.github/workflows/public-evidence.yml` runs two opt-in steps whole, on public data only, and writes what they measured to the job summary. It is not a gate; run it by hand (`workflow_dispatch`), and a push to a branch other than main that changes the workflow or `bin/archaic_report.py` runs it too. +`.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 each container's peak memory, 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). diff --git a/tests/test_archaic_report.py b/tests/test_archaic_report.py index c1b76b31..f60ebc21 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -178,7 +178,8 @@ def image(var): 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", 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" From 0a7867ec2cd6251c09e8edcbc256bd2157a11d60 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Tue, 3 Nov 2026 18:28:19 +0100 Subject: [PATCH 10/10] docs(archaic): drop the published-band verdict from the evidence summary, and give the AADR download size again The step 43 docs say not to compare the share with published figures, so the job summary no longer does. The optional resources table lists downloads, so the AADR row gives the download again. The sampler keeps one peak per step, and testing.md says so. --- .github/workflows/public-evidence.yml | 3 +-- docs/hardware-requirements.md | 2 +- docs/testing.md | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 01c02b9d..8d4318ab 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -284,9 +284,8 @@ jobs: [ "$(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) - band=$(awk -v x="$nea" 'BEGIN {print (x >= 1 && x <= 2) ? "inside" : "outside"}') { - 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. The Neanderthal share is ${band} the published 1-2% for people outside Africa." + 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 diff --git a/docs/hardware-requirements.md b/docs/hardware-requirements.md index 6fa5522b..47847550 100644 --- a/docs/hardware-requirements.md +++ b/docs/hardware-requirements.md @@ -68,7 +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) | ~1.05 GB (9.6 GiB at the peak) | step 44, `setup.sh --ancient-panel` | +| 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) diff --git a/docs/testing.md b/docs/testing.md index e098f575..08c171cd 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -168,7 +168,7 @@ The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub` `.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 each container's peak memory, 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). +- `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