diff --git a/.github/workflows/public-evidence.yml b/.github/workflows/public-evidence.yml index 055335c6..0126a66f 100644 --- a/.github/workflows/public-evidence.yml +++ b/.github/workflows/public-evidence.yml @@ -20,6 +20,9 @@ name: Public evidence # one an SNV. # The job "empty" is the red control: an empty region # list must make it fail. +# archaic-same-person the method check: one 1000 Genomes person with +# Chen et al.'s released calls, step 43 on her high- and +# low-coverage genotypes, beside the released figure. # 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 @@ -38,6 +41,18 @@ on: description: "Run the ancient-footprint job too" type: boolean default: true + hg002: + description: "Run the archaic-hg002 and archaic-share jobs (the whole-genome HG002 run)" + type: boolean + default: true + same_person: + description: "1000 Genomes person for the archaic-same-person job (Chen et al. released calls exist for them)" + type: string + default: NA12878 + same_person_chromosomes: + description: "Autosomes for the archaic-same-person job, separated by spaces or commas, or none" + type: string + default: none permissions: contents: read @@ -97,6 +112,7 @@ jobs: retention-days: 3 archaic-hg002: + if: inputs.hg002 name: archaic-hg002 (${{ matrix.part }}) needs: archaic-sites runs-on: ubuntu-latest @@ -276,6 +292,180 @@ jobs: path: ${{ env.EV }}/out/${{ matrix.part }} if-no-files-found: error + archaic-same-person: + # The method check behind the whole-genome gate: one 1000 Genomes person + # whose IBDmix calls Chen et al. 2020 released, run through step 43 on two + # inputs for the same person (the NYGC 30x high-coverage genotypes, and the + # original low-coverage phase 3 data called again on GRCh38), beside the + # released figure for that person and for the 1000 Genomes group. + needs: archaic-sites + if: inputs.same_person_chromosomes != 'none' + runs-on: ubuntu-latest + timeout-minutes: 350 + env: + SAMPLE: ${{ inputs.same_person }} + CHROMS: ${{ inputs.same_person_chromosomes }} + CHEN_ID: 125XCvTpQbK6TPX84y4442ZaDLDbxWMCf + CHEN_SHA256: 454bfa32f6b3cf4b48f0466b601054f5ed141fe61fea3aff4d3284a7a1f4dc5b + HC_URL: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV + LC_URL: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL + 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" + read -ra want <<< "$(printf '%s' "$CHROMS" | tr ',\r\n' ' ')" + for c in "${want[@]}"; do + [[ "$c" =~ ^([1-9]|1[0-9]|2[0-2])$ ]] || { echo "same_person_chromosomes: '${c}' is not an autosome 1 to 22"; exit 1; } + done + [[ "$SAMPLE" =~ ^(HG|NA)[0-9]+$ ]] || { echo "same_person: '${SAMPLE}' is not a 1000 Genomes sample"; exit 1; } + # each chromosome once, in order + echo "LIST=$(printf '%s\n' "${want[@]}" | sort -nu | tr '\n' ' ' | sed 's/ $//')" >> "$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; } + 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 + + - uses: actions/download-artifact@9000827ccba6bdab643e8b6fd33ac0654aef8333 # v8.0.2 + with: + name: archaic-1kg-sites + path: sites + + - name: Chen et al.'s released calls (sha256 pinned) + run: | + set -euo pipefail + mkdir -p "${EV}/chen" + curl -fsSL --retry 5 -o "${EV}/chen/chen.txt" "https://drive.usercontent.google.com/download?id=${CHEN_ID}&export=download&confirm=t" + [ "$(sha256sum "${EV}/chen/chen.txt" | cut -d' ' -f1)" = "$CHEN_SHA256" ] || { echo "sha256 of Chen's file is not the pinned one"; head -c 300 "${EV}/chen/chen.txt"; exit 1; } + awk -v s="$SAMPLE" '$9 == s {n++} END {print s, n+0, "released segments"; exit n == 0}' "${EV}/chen/chen.txt" + + - name: setup.sh --archaic for the chromosomes + env: + THREADS: '4' + run: | + set -euo pipefail + . ./versions.env + G="${EV}/genome" + DATA="${G}/reference/archaic/${ARCHAIC_DATA_VERSION}" + mkdir -p "${DATA}/download" + mv sites/1kg_sites.vcf.gz sites/1kg_sites.vcf.gz.tbi "${DATA}/download/" + : > "${EV}/regions37.bed" + for c in $LIST; do printf '%s\t0\t300000000\n' "$c" >> "${EV}/regions37.bed"; done + for c in $LIST; do + 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 + done + ARCHAIC_REGIONS="${EV}/regions37.bed" scripts/setup.sh --archaic "$G" 2>&1 | tee "${EV}/setup.log" + [ "${PIPESTATUS[0]}" -eq 0 ] + test -s "${DATA}/DONE" + + - name: The person on both inputs, step 43 on each, and the comparison + env: + THREADS: '4' + run: | + set -euo pipefail + . ./versions.env + G="${EV}/genome" OUT="${EV}/out" + mkdir -p "$OUT" + # sub LABEL URL-PREFIX NAME-OF-FILE-PATTERN: the person's non-reference + # SNVs on the chromosomes, one VCF. + pull() { + local label=$1 url=$2 pat=$3 c parts="" + mkdir -p "${G}/${SAMPLE}${label}/vcf" + for c in $LIST; do + f="${pat//@/${c}}" + # The container has no CA bundle, so curl fetches the file here. + curl -fsSL --retry 5 -o "${EV}/in.vcf.gz" "${url}/${f}" + curl -fsSL --retry 5 -o "${EV}/in.vcf.gz.tbi" "${url}/${f}.tbi" + docker run --rm -u "$(id -u):$(id -g)" -v "${EV}:/e" "$BCFTOOLS_IMAGE" \ + bcftools view -s "$SAMPLE" -v snps --min-ac 1 -Oz -o "/e/${label}.raw.vcf.gz" /e/in.vcf.gz + # The 2019 recall names the contigs 22, not chr22: rename them + # (a name already chrN stays as it is). + printf '%s\tchr%s\n' "$c" "$c" > "${EV}/rename.txt" + docker run --rm -u "$(id -u):$(id -g)" -v "${EV}:/e" "$BCFTOOLS_IMAGE" \ + bcftools annotate --rename-chrs /e/rename.txt -Oz -o "/e/${label}.chr${c}.vcf.gz" "/e/${label}.raw.vcf.gz" + rm -f "${EV}/${label}.raw.vcf.gz" + rm -f "${EV}/in.vcf.gz" "${EV}/in.vcf.gz.tbi" + parts="${parts} /e/${label}.chr${c}.vcf.gz" + done + # shellcheck disable=SC2086 + docker run --rm -u "$(id -u):$(id -g)" -v "${EV}:/e" -v "${G}:/g" "$BCFTOOLS_IMAGE" \ + bcftools concat -Oz -o "/g/${SAMPLE}${label}/vcf/${SAMPLE}${label}.vcf.gz" $parts + docker run --rm -u "$(id -u):$(id -g)" -v "${G}:/g" "$BCFTOOLS_IMAGE" \ + bcftools index -t -f "/g/${SAMPLE}${label}/vcf/${SAMPLE}${label}.vcf.gz" + rm -f "${EV}/${label}".chr*.vcf.gz + echo "${label}: $(docker run --rm -v "${G}:/g" "$BCFTOOLS_IMAGE" bcftools view -H "/g/${SAMPLE}${label}/vcf/${SAMPLE}${label}.vcf.gz" | wc -l) SNV records" + } + pull hc "$HC_URL" '1kGP_high_coverage_Illumina.chr@.filtered.SNV_INDEL_SV_phased_panel.vcf.gz' + pull lc "$LC_URL" 'ALL.chr@.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz' + args=() + for label in hc lc; do + for pop in ALL EUR; do + ARCHAIC_POP=$pop GENOME_DIR="$G" scripts/43-archaic.sh "${SAMPLE}${label}" 2>&1 | tee "${OUT}/step43-${label}-${pop}.log" + cp "${G}/${SAMPLE}${label}/archaic/${SAMPLE}${label}_archaic_segments.tsv" "${OUT}/${label}-${pop}.segments.tsv" + cp "${G}/${SAMPLE}${label}/archaic/${SAMPLE}${label}_archaic.tsv" "${OUT}/${label}-${pop}.summary.tsv" + mkdir -p "${OUT}/ibdmix-${label}-${pop}" + cp "${G}/${SAMPLE}${label}/archaic/work/ibdmix/"chr*.txt "${G}/${SAMPLE}${label}/archaic/work/ibdmix/spans.tsv" "${OUT}/ibdmix-${label}-${pop}/" + args+=(--mine "${label}-${pop}=${OUT}/${label}-${pop}.segments.tsv") + done + done + grp=$(awk -v s="$SAMPLE" '$9 == s {print $8; exit}' "${EV}/chen/chen.txt") + python3 -I bin/archaic_compare.py compare --chen-file "${EV}/chen/chen.txt" --sample "$SAMPLE" --group "$grp" \ + --chroms "$LIST" "${args[@]}" | tee "${OUT}/comparison.md" + for label in hc lc; do + { + echo "IBDmix output of ${label} (ALL frequencies), per archaic genome:" + echo "" + python3 -I bin/archaic_compare.py genomes --chen-file "${EV}/chen/chen.txt" --sample "$SAMPLE" \ + --ibdmix-dir "${OUT}/ibdmix-${label}-ALL" --spans "${OUT}/ibdmix-${label}-ALL/spans.tsv" + echo "" + } | tee -a "${OUT}/genomes.md" + done + { + echo "### Same person: ${SAMPLE} (${grp}), chromosomes ${LIST}" + echo "" + echo "hc = NYGC 30x high-coverage genotypes (GRCh38); lc = the original low-coverage phase 3 data called again on GRCh38. ALL and EUR are the 1000 Genomes frequencies step 43 used." + echo "" + cat "${OUT}/comparison.md" + echo "" + cat "${OUT}/genomes.md" + } >> "$GITHUB_STEP_SUMMARY" + + - name: Upload the segments and the comparison + if: always() + uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 + with: + name: archaic-same-person + path: ${{ env.EV }}/out + if-no-files-found: error + archaic-share: needs: [archaic-sites, archaic-hg002] runs-on: ubuntu-latest @@ -334,13 +524,15 @@ jobs: nea=$(kv neanderthal_percent) mb=$(awk -v b="$(kv neanderthal_bp)" 'BEGIN {printf "%.1f", b / 1e6}') # Chen et al. 2020 (Cell 180:677), released per-person calls for - # Europeans: the union of the Altai and Vindija segments, LOD 4, - # 50 kb or longer, the cutoffs of step 43 (docs/43-archaic.md). + # the 503 Europeans: 42.0 to 62.0 Mb, median 51.6 Mb, made with the + # Altai genome alone and a refined call set. Step 43 takes the + # union of Altai and Vindija without the refinements, so its total + # is larger by design (docs/43-archaic.md, same-person comparison). if [ "$WANT" -eq 22 ]; then - where=$(awk -v b="$(kv neanderthal_bp)" 'BEGIN {print (b < 49.9e6 ? "below" : (b > 82.9e6 ? "above" : "inside"))}') - cmp="HG002's ${mb} Mb of Neanderthal sequence is ${where} the published range: Chen et al. 2020 released per-person calls with the same method and cutoffs (Altai and Vindija union, LOD 4, 50 kb or longer), and Europeans carry 49.9 to 82.9 Mb in them (see docs/43-archaic.md)." + ratio=$(awk -v b="$(kv neanderthal_bp)" 'BEGIN {printf "%.1f", b / 51.6e6}') + cmp="HG002's ${mb} Mb of Neanderthal sequence is ${ratio} times the median (51.6 Mb) of the per-person totals Chen et al. 2020 released for 503 Europeans (42.0 to 62.0 Mb). Their method uses the Altai genome alone and a refined call set; step 43 takes the union of Altai and Vindija, so it finds more by design (see docs/43-archaic.md)." else - cmp="HG002 has ${mb} Mb of Neanderthal sequence on these ${WANT} autosomes. The published range of 49.9 to 82.9 Mb (Chen et al. 2020, Europeans, Altai and Vindija union, LOD 4, 50 kb or longer) is for the 22 autosomes, so it is not compared here." + cmp="HG002 has ${mb} Mb of Neanderthal sequence on these ${WANT} autosomes. The per-person totals of Chen et al. 2020 (503 Europeans: 42.0 to 62.0 Mb) are for the 22 autosomes, so they are not compared here." fi { 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. ${cmp}" diff --git a/bin/archaic_compare.py b/bin/archaic_compare.py new file mode 100644 index 00000000..914f10f1 --- /dev/null +++ b/bin/archaic_compare.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python3 +"""Compare step 43's Neanderthal totals with Chen et al. 2020's released calls. + +Used by the archaic jobs of public-evidence.yml only (not part of the +pipeline). Standard library only. + +The released calls are Chen et al.'s file "Neanderthal sequence in 1000 +genome.50kb.txt" (chr, start, end, LOD, maxLOD, size, pop, anc, ID, tab +separated), merged per person and chromosome. + + compare the same person run through step 43 on one or more inputs (a + segments TSV), beside the released figure, and how much of each + lies under the released calls. + genomes the raw IBDmix output of one step 43 run, per archaic genome + (Altai, Vindija, Denisova) at LOD 4 and 50 kb, and the union, each + beside the released calls. + +Run: python3 bin/archaic_compare.py compare --chen-file F --sample NA12878 \ + --group EUR --chroms "21 22" --mine label=SEGMENTS.tsv [--mine label=...] +""" +import argparse +import collections +import os +import statistics +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import archaic_report # noqa: E402 (the same reader and cutoffs as step 43) + + +def norm(c): + """The chromosome name without a leading chr.""" + return c[3:] if c.startswith("chr") else c + + +def merged_bp(ivs): + """Total length of the union of (start, end) intervals.""" + total, cur = 0, None + for s, e in sorted(ivs): + if cur and s <= cur[1]: + cur[1] = max(cur[1], e) + else: + if cur: + total += cur[1] - cur[0] + cur = [s, e] + return total + (cur[1] - cur[0] if cur else 0) + + +def intersect_bp(a, b): + """Base pairs under both interval lists (each merged first).""" + def merged(ivs): + """The intervals with overlaps joined.""" + out = [] + for s, e in sorted(ivs): + if out and s <= out[-1][1]: + out[-1][1] = max(out[-1][1], e) + else: + out.append([s, e]) + return out + a, b, i, j, total = merged(a), merged(b), 0, 0, 0 + while i < len(a) and j < len(b): + lo, hi = max(a[i][0], b[j][0]), min(a[i][1], b[j][1]) + if hi > lo: + total += hi - lo + if a[i][1] < b[j][1]: + i += 1 + else: + j += 1 + return total + + +def read_chen(path): + """({person: {chrom: [(start, end)]}}, {person: group}).""" + calls = collections.defaultdict(lambda: collections.defaultdict(list)) + group = {} + with open(path) as f: + for line in f: + c = line.split() + if len(c) < 9 or c[0] == "chr": + continue + calls[c[8]][norm(c[0])].append((int(c[1]), int(c[2]))) + group[c[8]] = c[7] + return calls, group + + +def read_mine(path): + """{chrom: [(start, end)]} of the Neanderthal rows of a segments TSV.""" + out = collections.defaultdict(list) + with open(path) as f: + head = f.readline().rstrip("\n").split("\t") + if head[:5] != ["chrom", "start", "end", "length_bp", "ancestry"]: + sys.exit(f"{path}: not a step 43 segments file") + for line in f: + c = line.rstrip("\n").split("\t") + if c[4] == "Neanderthal": + out[norm(c[0])].append((int(c[1]), int(c[2]))) + return out + + +def chrom_order(c): + """Sort key: numbered chromosomes first, in order.""" + return (0, int(c)) if c.isdigit() else (1, c) + + +def cmd_compare(a): + """Print the per-chromosome table and the totals of the compare command.""" + calls, group = read_chen(a.chen_file) + if a.sample not in calls: + sys.exit(f"{a.sample} is not in {a.chen_file}") + peers = [p for p, g in group.items() if g == a.group] + if not peers: + sys.exit(f"no person of group {a.group} in {a.chen_file}") + ivs = collections.defaultdict(lambda: collections.defaultdict(list)) + labels = [] + for item in a.mine: + label, _, path = item.partition("=") + if label not in labels: + labels.append(label) + for c, v in read_mine(path).items(): + ivs[label][c].extend(v) + mine = {l: {c: merged_bp(v) for c, v in ivs[l].items()} for l in labels} + chroms = sorted({norm(c) for c in a.chroms.replace(",", " ").split()}, key=chrom_order) + if not chroms: + sys.exit("--chroms names no chromosome") + + def chen_bp(p, c): + """One person's merged base pairs on one chromosome.""" + return merged_bp(calls[p].get(c, [])) + + print(f"| chromosome | Chen {a.sample} (bp) | Chen {a.group} median | Chen {a.group} max | " + + " | ".join(f"{l} (bp)" for l in labels) + " |") + print("|---|---|---|---|" + "---|" * len(labels)) + tot_me, tot_med, tot_max = 0, 0.0, 0 + tot_mine = dict.fromkeys(labels, 0) + for c in chroms: + vals = [chen_bp(p, c) for p in peers] + me, med, mx = chen_bp(a.sample, c), statistics.median(vals), max(vals) + tot_me, tot_med, tot_max = tot_me + me, tot_med + med, tot_max + mx + cells = [] + for l in labels: + v = mine[l].get(c, 0) + tot_mine[l] += v + cells.append(str(v)) + print(f"| {c} | {me} | {med:.0f} | {mx} | " + " | ".join(cells) + " |") + print(f"| all | {tot_me} | {tot_med:.0f} | {tot_max} | " + " | ".join(str(tot_mine[l]) for l in labels) + " |") + print() + # The person-level spread over the same chromosomes. + sums = sorted(sum(chen_bp(p, c) for c in chroms) for p in peers) + print(f"Chen et al., {len(peers)} {a.group} people, these {len(chroms)} chromosomes: " + f"lowest {sums[0]}, median {statistics.median(sums):.0f}, highest {sums[-1]} bp; " + f"{a.sample} {tot_me} bp.") + for l in labels: + ratio = tot_mine[l] / tot_me if tot_me else float("nan") + both = sum(intersect_bp(ivs[l].get(c, []), calls[a.sample].get(c, [])) for c in chroms) + print(f"{l}: {tot_mine[l]} bp, {ratio:.2f} times the released {a.sample} figure; " + f"{both} bp under the released calls, {tot_mine[l] - both} bp only here, {tot_me - both} bp only released.") + + +def cmd_genomes(a): + """Print the per-archaic-genome table of the genomes command.""" + calls, _ = read_chen(a.chen_file) + if a.sample not in calls: + sys.exit(f"{a.sample} is not in {a.chen_file}") + chroms = [] + with open(a.spans) as f: + next(f) + chroms = [norm(line.split("\t")[0]) for line in f] + per = collections.defaultdict(lambda: collections.defaultdict(list)) + for c in chroms: + for n in archaic_report.ARCHAICS: + segs = archaic_report.keep(archaic_report.read_segments(os.path.join(a.ibdmix_dir, f"chr{c}.{n}.txt"))) + per[n][c] += [(s, e) for _, s, e, _ in segs] + per["Altai+Vindija"] = {c: per["Altai"][c] + per["Vindija"][c] for c in chroms} + released = sum(merged_bp(calls[a.sample].get(c, [])) for c in chroms) + print(f"| IBDmix output | bp | under the released {a.sample} calls | only here | only released |") + print("|---|---|---|---|---|") + for n in ("Altai", "Vindija", "Denisova", "Altai+Vindija"): + mine = sum(merged_bp(per[n][c]) for c in chroms) + both = sum(intersect_bp(per[n][c], calls[a.sample].get(c, [])) for c in chroms) + print(f"| {n} | {mine} | {both} | {mine - both} | {released - both} |") + print(f"\nReleased {a.sample}: {released} bp on these {len(chroms)} chromosomes.") + + +def main(argv=None): + """Parse the command line and run the chosen command.""" + p = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + sub = p.add_subparsers(dest="cmd", required=True) + c = sub.add_parser("compare") + c.add_argument("--chen-file", required=True) + c.add_argument("--sample", required=True) + c.add_argument("--group", default="EUR") + c.add_argument("--chroms", required=True, help="the chromosomes analysed, separated by spaces or commas") + c.add_argument("--mine", action="append", required=True, metavar="LABEL=SEGMENTS.tsv") + c.set_defaults(fn=cmd_compare) + g = sub.add_parser("genomes") + g.add_argument("--chen-file", required=True) + g.add_argument("--sample", required=True) + g.add_argument("--ibdmix-dir", required=True) + g.add_argument("--spans", required=True) + g.set_defaults(fn=cmd_genomes) + a = p.parse_args(argv) + a.fn(a) + + +if __name__ == "__main__": + main() diff --git a/docs/43-archaic.md b/docs/43-archaic.md index 7dbd38a6..c25d4055 100644 --- a/docs/43-archaic.md +++ b/docs/43-archaic.md @@ -63,7 +63,7 @@ The first version of the masks dropped the masked calls instead (code 9, or no r One of their exclusions is not applied here: the CpG mask, which needs a GRCh37 reference and the chimpanzee, orangutan and macaque SNP sets, none of which this pipeline installs. What that leaves out is written down in the Notes below. -The filters cut hard: on chr22 the strict mask covers 25.2 Mb of 50.8, the Vindija FilterBed 19.6 Mb, and segmental duplications take out 4.3 Mb. Without them the step called sequence from positions the archaic genomes cannot be read at, and HG002's Neanderthal total came out about twice the published figure (see Runtime below). The filters are part of the data version: a change to a pin or to a mask means `snpad2016-1kgp3-v3` and a rebuild. `setup.sh --archaic` also rebuilds a table whose build code changed (its stamp holds the sha256 of the table half of `bin/archaic_report.py`), so a table built before the masked codes, or before the strict-mask rows were kept, is rebuilt on the next `setup.sh --archaic`. Until then step 43 reads it as before: a table with no masked code gives empty masks. +The filters cut hard: on chr22 the strict mask covers 25.2 Mb of 50.8, the Vindija FilterBed 19.6 Mb, and segmental duplications take out 4.3 Mb. Without them the step called sequence from positions the archaic genomes cannot be read at, and HG002's Neanderthal total was higher still: the filters took it from 115.2 Mb to 107.5 Mb, which is still 2.1 times the published median (see the same-person section under Runtime). The filters are part of the data version: a change to a pin or to a mask means `snpad2016-1kgp3-v3` and a rebuild. `setup.sh --archaic` also rebuilds a table whose build code changed (its stamp holds the sha256 of the table half of `bin/archaic_report.py`), so a table built before the masked codes, or before the strict-mask rows were kept, is rebuilt on the next `setup.sh --archaic`. Until then step 43 reads it as before: a table with no masked code gives empty masks. The table is built on GRCh37 and lifted once to GRCh38 with Picard `LiftoverVcf`. A site that does not lift, or lifts to another chromosome, is dropped. The archaic files are downloaded one chromosome at a time (150 GB in all, at most 13.8 GB at once, for chromosome 2) and deleted once the chromosome is built, so an interrupted install resumes at the chromosome it stopped in. `ARCHAIC_REGIONS` (a GRCh37 BED, chromosomes named 1 to 22) builds only those regions. @@ -99,9 +99,26 @@ 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. -That run was made before the masks, from a table that held every lifted site. It gave HG002 115.2 Mb of Neanderthal sequence over the 22 autosomes, about twice what Chen et al. report for Europeans with the same method and cutoffs (about 51 Mb per person, and 49.9 to 82.9 Mb across their released per-person calls). The cause was the missing masks: the table carried the positions each archaic genome's own filters exclude, and IBDmix scored them. On chr22, 43% of that run's Neanderthal sequence lay outside the 1000 Genomes GRCh38 strict mask alone. +That run was made before the masks, from a table that held every lifted site. It gave HG002 115.2 Mb of Neanderthal sequence over the 22 autosomes, about twice what Chen et al. report for Europeans (a median of 51.6 Mb per person, and 42.0 to 62.0 Mb across the released per-person calls of their 503 Europeans; see the same-person comparison below for what their method is). The cause was the missing masks: the table carried the positions each archaic genome's own filters exclude, and IBDmix scored them. On chr22, 43% of that run's Neanderthal sequence lay outside the 1000 Genomes GRCh38 strict mask alone. -With the masks dropping sites, a run on chr22 (run 42089578470, before the masked codes) gave HG002 1.3 Mb of Neanderthal sequence with the `ALL` frequencies and 1.7 Mb with `EUR`, still at the top of the published spread below. Keeping the masked sites as IBDmix's mask, as the published workflow does, lets the mismatches there end a segment; the chr22 figure on this table is printed by the e2e case on every run. The e2e case compares it with the spread of Chen et al.'s own released calls for the 497 Europeans of CEU, TSI, GBR, FIN and IBS (51 kb at the lowest, 541 kb at the median, 1,132 kb at the 95th percentile, 1,716 kb at the highest). The pre-mask run gave 1,748 kb on chr22, more than all 497 of them. With the masked codes but the strict-mask rows still dropped, the fork runs gave 1,177 to 1,190 kb. Keeping those rows marked, the last difference from the published workflow in how masked sites are handled, gave 1,138 kb on chr22 (run 42335592918). One chromosome of one person varies widely, so the case gates on a regression guard of 1,200 kb (and on the lowest published value, 51 kb), and prints the published figures beside the measured one. It does not claim the figure is inside the published range. The whole-genome total (49.9 to 82.9 Mb across the published per-person calls) is checked by the release run of `public-evidence.yml`, which installs the 150 GB of archaic data once per release, not by the e2e case; that run decides whether the share is in range. The case prints where HG002 falls inside the spread; it does not gate on the median or the 95th percentile, for two reasons: one person's chromosome is one draw from that distribution, and the published calls come from imputed low-coverage genotypes, which lose rare alleles and so lose archaic ones, while HG002 is a high-coverage benchmark. +With the masks dropping sites, a run on chr22 (run 42089578470, before the masked codes) gave HG002 1.3 Mb of Neanderthal sequence with the `ALL` frequencies and 1.7 Mb with `EUR`, still at the top of the published spread below. Keeping the masked sites as IBDmix's mask, as the published workflow does, lets the mismatches there end a segment; the chr22 figure on this table is printed by the e2e case on every run. The e2e case compares it with the spread of Chen et al.'s own released calls for the 497 Europeans of CEU, TSI, GBR, FIN and IBS (51 kb at the lowest, 541 kb at the median, 1,132 kb at the 95th percentile, 1,716 kb at the highest). The pre-mask run gave 1,748 kb on chr22, more than all 497 of them. With the masked codes but the strict-mask rows still dropped, the fork runs gave 1,177 to 1,190 kb. Keeping those rows marked, the last difference from the published workflow in how masked sites are handled, gave 1,138 kb on chr22 (run 42335592918). One chromosome of one person varies widely, so the case gates on a regression guard of 1,200 kb (and on the lowest published value, 51 kb), and prints the published figures beside the measured one. It does not claim the figure is inside the published range. The whole-genome total is measured by the release run of `public-evidence.yml`, which installs the 150 GB of archaic data once per release, not by the e2e case, and printed beside the published 42.0 to 62.0 Mb; it is not gated on it, for the reasons in the same-person comparison below. The case prints where HG002 falls inside the spread; it does not gate on the median or the 95th percentile, because one person's chromosome is one draw from that distribution. The same-person comparison below shows the input's coverage is not why the step finds more than the published calls do. + +### Same person as Chen et al.: the method differs, not the data + +The whole-genome run on main (run 42596685257, after the masks) gave HG002 107.5 Mb of Neanderthal sequence with the `ALL` frequencies and 109.4 Mb with `EUR`, 2.1 times the median of the published per-person calls. The open question was whether that is the input (HG002 is a high-coverage benchmark, the published calls come from the low-coverage 1000 Genomes data, which loses rare alleles) or the method. The `archaic-same-person` job of `public-evidence.yml` answers it with one person both sides have: NA12878 (CEU), whose IBDmix calls Chen et al. released (file `Neanderthal sequence in 1000 genome.50kb.txt` of their Google Drive folder, sha256 pinned in the job). It runs step 43 on her NYGC 30x high-coverage genotypes and on the original low-coverage phase 3 data called again on GRCh38, over chromosomes 19 to 22 (run https://github.com/geiresx/OwnGenome/actions/runs/42645379472): + +| Neanderthal bp, chr19 to 22 | Altai+Vindija, `ALL` | Altai+Vindija, `EUR` | +|---|---|---| +| step 43, 30x high-coverage input | 6,027,711 | 5,667,159 | +| step 43, low-coverage input | 6,051,348 | 5,586,970 | +| Chen et al., same person | 1,886,257 | 1,886,257 | +| Chen et al., 503 Europeans: lowest, median, highest | 846,101 / 2,165,624 / 3,987,644 | | + +- **The input is not the cause.** The two inputs differ by 0.4% with `ALL` and 1.4% with `EUR`. Step 43 finds 3.0 to 3.2 times her released figure on either, so the hypothesis that high coverage inflates the total is refuted. +- **The released calls are on GRCh37.** (This overlap was computed offline from the run's artifact and the UCSC chain; no workflow step reproduces it.) Lifted to GRCh38 with the UCSC chain, her released calls cover 1,885,100 bp, and 98% of them (1,852,663 bp) lie inside step 43's Altai+Vindija segments, 95% (1,794,145 bp) inside the Altai segments alone. Step 43 finds everything the released calls find, plus 4.2 Mb more on these four chromosomes. So the difference is calls the published set does not contain, not calls it has and step 43 misses. +- **Where the extra comes from.** The paper's method differs from step 43 in ways step 43 does not reproduce: (1) it uses the Altai Neanderthal only (and the Altai Denisovan), where step 43 takes the union of Altai and Vindija; on the same input the Altai segments alone are 3,819,442 bp and Vindija's alone 5,001,594 bp, against 6,027,711 for the union; (2) it runs IBDmix per 1000 Genomes subgroup (CEU), where step 43 uses the frequencies of all 2,504 people or of the continent; (3) it refines the Neanderthal call set by removing regions called Denisovan in Africans and Neanderthal anywhere, and regions in the top 0.1% of derived-allele proportion in the Neanderthal genome, which it reports lowers each population by several Mb; (4) it masks CpG sites, which step 43 does not (see Notes). With the Altai genome alone the total is still 2.0 times her released figure, which leaves the refinements and the CpG mask. + +What this means for the figure step 43 reports. It is the union of the Altai and Vindija segments at LOD 4 and 50 kb, without the paper's call-set refinements, so it is larger than the published per-person totals by design: 2.0 times on the Altai genome alone and 3.2 times on the union for NA12878 over four chromosomes, and 2.1 times the published median for HG002 over the genome. The published 42.0 to 62.0 Mb (median 51.6 Mb) is a figure for the paper's method; it is not a range step 43's number should fall into, and the earlier range of 49.9 to 82.9 Mb was not that either (it was stated for the Altai and Vindija union, which the paper does not use). The step's number is an upper-end estimate of the Neanderthal sequence a person carries, and compares between people run through the same step, not with a published total. ## Notes diff --git a/docs/testing.md b/docs/testing.md index e3ef3253..0ecbfa15 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -197,7 +197,7 @@ gh workflow run public-evidence.yml -f chromosomes=22 -f ancient=false ``` - `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. When all 22 autosomes ran, the summary says whether HG002's Neanderthal total lies inside the 49.9 to 82.9 Mb that Chen et al.'s released per-person calls give Europeans with the same method and cutoffs; it prints the comparison and does not gate on it. 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). +- `archaic-hg002`: `setup.sh --archaic` and step 43 on HG002, one job per autosome; `archaic-share` joins them into one report per population. When all 22 autosomes ran, the summary prints HG002's total as a multiple of the 51.6 Mb median of Chen et al.'s 503 European per-person totals (42.0 to 62.0 Mb), not an inside or outside verdict, and does not gate on it. The `archaic-same-person` job takes the `hg002` and `same_person` inputs and compares step 43 with the released calls of one person both sides have. 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 diff --git a/tests/test_archaic_compare.py b/tests/test_archaic_compare.py new file mode 100644 index 00000000..ad09e9b7 --- /dev/null +++ b/tests/test_archaic_compare.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +"""bin/archaic_compare.py: merging, Chen's file format, and the comparison table.""" +import contextlib +import io +import os +import subprocess +import sys +import tempfile +import unittest + +HERE = os.path.dirname(os.path.abspath(__file__)) +SCRIPT = os.path.join(HERE, "..", "bin", "archaic_compare.py") +sys.path.insert(0, os.path.join(HERE, "..", "bin")) +import archaic_compare as ac # noqa: E402 + +CHEN = ( + "chr start end LOD maxLOD size pop anc ID \n" + "22\t100\t200\t9\t9\t100\tCEU\tEUR\tA\n" + "22\t150\t300\t9\t9\t150\tCEU\tEUR\tA\n" # overlaps the first: union 200 bp + "22\t1000\t1100\t9\t9\t100\tCEU\tEUR\tA\n" + "21\t10\t60\t9\t9\t50\tCEU\tEUR\tA\n" + "22\t100\t150\t9\t9\t50\tGBR\tEUR\tB\n" + "22\t100\t900\t9\t9\t800\tYRI\tAFR\tC\n" +) +MINE = ( + "chrom\tstart\tend\tlength_bp\tancestry\tmatched\tlod\tgenes\n" + "chr22\t100\t400\t300\tNeanderthal\tAltai\t5\t\n" + "chr22\t5000\t6000\t1000\tDenisovan\tDenisova\t5\t\n" +) + + +class Compare(unittest.TestCase): + def setUp(self): + self.d = tempfile.mkdtemp() + self.chen = os.path.join(self.d, "chen.txt") + self.mine = os.path.join(self.d, "mine.tsv") + open(self.chen, "w").write(CHEN) + open(self.mine, "w").write(MINE) + + def run_cmp(self, *extra): + buf = io.StringIO() + with contextlib.redirect_stdout(buf): + ac.main(["compare", "--chen-file", self.chen, "--sample", "A", "--group", "EUR", + "--chroms", "21 22", "--mine", f"x={self.mine}", *extra]) + return buf.getvalue() + + def test_merged_bp(self): + self.assertEqual(ac.merged_bp([(100, 200), (150, 300), (1000, 1100)]), 300) + self.assertEqual(ac.merged_bp([]), 0) + + def test_person_total_is_the_union(self): + out = self.run_cmp() + # chr22: 100-300 and 1000-1100 = 300; chr21: 50 + self.assertIn("| 22 | 300 |", out) + self.assertIn("| 21 | 50 |", out) + self.assertIn("| all | 350 |", out) + + def test_only_neanderthal_rows_of_mine_count(self): + out = self.run_cmp() + self.assertIn("x: 300 bp, 0.86 times the released A figure;", out) + + def test_overlap_with_the_released_calls(self): + out = self.run_cmp() + # mine chr22 100-400; released chr22 100-300 and 1000-1100: 200 bp shared + self.assertIn("200 bp under the released calls, 100 bp only here, 150 bp only released", out) + + def test_intersect_bp(self): + self.assertEqual(ac.intersect_bp([(0, 10), (20, 30)], [(5, 25)]), 10) + self.assertEqual(ac.intersect_bp([(0, 10)], []), 0) + + def test_group_spread_uses_the_group_only(self): + out = self.run_cmp() + self.assertIn("2 EUR people", out) # A and B, not the AFR person C + self.assertIn("highest 350", out) + + def test_unknown_sample_stops(self): + with self.assertRaises(SystemExit): + ac.main(["compare", "--chen-file", self.chen, "--sample", "Z", "--chroms", "22", + "--mine", f"x={self.mine}"]) + + def test_empty_group_stops(self): + with self.assertRaises(SystemExit) as r: + ac.main( + ["compare", "--chen-file", self.chen, "--sample", "A", "--group", "NOPE", "--chroms", "22", + "--mine", f"x={self.mine}"]) + self.assertIn("NOPE", str(r.exception)) + + def test_wrong_segments_file_stops(self): + open(self.mine, "w").write("a\tb\n") + with self.assertRaises(SystemExit): + self.run_cmp() + + def test_genomes_union_row(self): + hdr = "ID\tchrom\tstart\tend\tslod\n" + rows = { + "Altai": "SAMPLE\t22\t100000\t300000\t5\n", + "Vindija": "SAMPLE\t22\t1000000\t1100000\t5\n", + "Denisova": "SAMPLE\t22\t5000000\t5100000\t5\n", + } + for n, body in rows.items(): + open(os.path.join(self.d, f"chr22.{n}.txt"), "w").write(hdr + body) + spans = os.path.join(self.d, "spans.tsv") + open(spans, "w").write("chrom\tstart\tend\n22\t1\t9999\n") + chen = os.path.join(self.d, "chen_big.txt") + open(chen, "w").write( + "chr start end LOD maxLOD size pop anc ID\n" + "22\t100000\t300000\t9\t9\t200000\tCEU\tEUR\tA\n" + "22\t1000000\t1100000\t9\t9\t100000\tCEU\tEUR\tA\n") + buf = io.StringIO() + with contextlib.redirect_stdout(buf): + ac.main(["genomes", "--chen-file", chen, "--sample", "A", + "--ibdmix-dir", self.d, "--spans", spans]) + out = buf.getvalue() + # released A on chr22: 300,000 bp; Altai and Vindija together are exactly that + self.assertIn("| Altai+Vindija | 300000 | 300000 | 0 | 0 |", out) + self.assertIn("| Denisova | 100000 | 0 | 100000 | 300000 |", out) + self.assertIn("Released A: 300000 bp on these 1 chromosomes.", out) + + def test_cli(self): + r = subprocess.run([sys.executable, "-I", SCRIPT, "compare", "--chen-file", self.chen, "--sample", "A", + "--chroms", "22", "--mine", f"x={self.mine}"], capture_output=True, text=True) + self.assertEqual(r.returncode, 0, r.stderr) + + +if __name__ == "__main__": + unittest.main()