fix(archaic): keep the SNV call when an insertion sits at the same position, and measure what the ancient panel really needs

bin/archaic_report.py

@@ -248,6 +248,34 @@ def read_vcf(path):
             yield c[0], int(c[1]), c[3], c[4], code
 
 
+CLASH = object()
+
+
+def pick_call(prev, call):
+    """The call to keep when a position comes twice, as in a VCF split with
+    bcftools norm -m-: each is (REF, ALT, code), ALT None for an insertion or a
+    multi-allelic record. An SNV wins over such a record, in either order. Of
+    one SNV twice, a usable record wins over code 9, and two different usable
+    codes are CLASH. Of two different SNVs, the
+    one that carries its ALT wins over a 0/0 or unusable one (a split 0/1 at
+    A>G,T is A>G 0/1 and A>T 0/0); both 0/0 is the reference; both carrying
+    their ALT, or neither readable, is CLASH (unknown, 9)."""
+    if prev is None:
+        return call
+    if prev[1] is None or call[1] is None:
+        return call if prev[1] is None and call[1] is not None else prev
+    if call[1] == prev[1]:
+        if call[2] == 9 or call[2] == prev[2]:
+            return prev
+        return call if prev[2] == 9 else CLASH
+    has_prev, has_call = prev[2] in (1, 2), call[2] in (1, 2)
+    if has_prev != has_call:
+        return prev if has_prev else call
+    if not has_prev and prev[2] == 0 and call[2] == 0:
+        return (prev[0], None, 0)
+    return CLASH
+
+
 def read_table(path):
     """Rows of one chromosome of the GRCh38 table, with the header checked."""
     with opener(path) as f:
@@ -290,9 +318,14 @@ def cmd_genotypes(a):
                 pending = next(vcf, None)
             continue
         # The sample's calls on this chromosome, merged with the table by position.
-        calls = {}
+        calls, clash = {}, set()
         while pending is not None and pending[0] == chrom:
-            calls[pending[1]] = pending[2:]
+            pos = pending[1]
+            if pos not in clash:
+                calls[pos] = pick_call(calls.get(pos), pending[2:])
+                if calls[pos] is CLASH:
+                    calls[pos] = (pending[2], None, 9)
+                    clash.add(pos)
             pending = next(vcf, None)
         n = first = last = dropped = carried = 0
         num = chrom[3:]

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-<release>/` (58 GB; about 72 GB free while it runs), deletes the archives and pulls `EXOMISER_IMAGE`. An interrupted download resumes when the command runs again. See [Step 48](48-exomiser.md).
 
+## AADR Ancient Panel (Optional, for Step 44, ~1.05 GB kept, 9.6 GiB at the peak)
+
+Only needed for step 44, the placement among ancient genomes. `./scripts/setup.sh --ancient-panel ${GENOME_DIR}` downloads the AADR release `AADR_VERSION` (1240K, the 7.1 GB `.geno` and three small files, each checked against its md5 in `versions.env`) and UCSC's hg19ToHg38 chain, then builds the panel in `${GENOME_DIR}/reference/ancient/aadr_<release>/` with convertf and two smartpca runs. On a 4-CPU GitHub runner it took 61 minutes, peaked at 9.6 GiB of extra disk and kept 1.05 GiB; its smartpca container gets 12 GB, and the setup's containers held 7.6 GiB of anonymous memory at most (any setup container, convertf included). See [Step 44](44-ancient.md#runtime).
+
 ## Docker Images — Pre-Pull All
 
 `setup.sh` pulls every image in `versions.env` except the lines marked `# optional`, and `scripts/setup.sh --pull-only` pulls the same list without the rest of setup. A failed pull is tried three times (`FETCH_TRIES` and `FETCH_WAIT` change that), and Docker's own message is printed: `toomanyrequests` is Docker Hub's rate limit, `manifest unknown` a tag that does not exist. `--parascopy-data`, `--yleaf-data` and `--megane-kmers` pull their step's optional image along with its data. Each image pinned by digest (`name:tag@sha256:...`) also gets a local tag, `localhost/pgp-pinned/<name>:sha256-<digest>`: 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):

docs/43-archaic.md

@@ -45,9 +45,9 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags.
 | 1000 Genomes phase 3 | `ALL.wgs.phase3_shapeit2_mvncall_integrated_v5c.20130502.sites.vcf.gz` (1.5 GB) | md5 from the release's `current.tree` |
 | UCSC | `hg19ToHg38.over.chain.gz` | sha256 |
 
-MPI-EVA publishes md5 lists, not sha256 sums, for 66 files of about 2.5 GB each. Pinning the sha256 of the three lists pins every file's md5 without downloading 165 GB to hash it.
+MPI-EVA publishes md5 lists, not sha256 sums, for 66 files of about 2.5 GB each. Pinning the sha256 of the three lists pins every file's md5 without downloading 150 GB to hash it.
 
-The table is built on GRCh37 and lifted once to GRCh38 with Picard `LiftoverVcf`. A site that does not lift, or lifts to another chromosome, is dropped. The archaic files are downloaded one chromosome at a time (about 165 GB in all, under 13 GB at once) and deleted once the chromosome is built, so an interrupted install resumes at the chromosome it stopped in. `ARCHAIC_REGIONS` (a GRCh37 BED, chromosomes named 1 to 22) builds only those regions.
+The table is built on GRCh37 and lifted once to GRCh38 with Picard `LiftoverVcf`. A site that does not lift, or lifts to another chromosome, is dropped. The archaic files are downloaded one chromosome at a time (150 GB in all, at most 13.8 GB at once, for chromosome 2) and deleted once the chromosome is built, so an interrupted install resumes at the chromosome it stopped in. `ARCHAIC_REGIONS` (a GRCh37 BED, chromosomes named 1 to 22) builds only those regions.
 
 ## Input
 
@@ -76,7 +76,12 @@ Both reports show the shares, the analysed span and the ten longest segments in
 
 ## Runtime
 
-Not measured on a whole genome. IBDmix runs once per chromosome and archaic genome, `THREADS` chromosomes at a time. The install is bound by the download.
+Measured on public HG002 by the `archaic-hg002` jobs of `public-evidence.yml`, one job per autosome on a GitHub runner with 4 CPUs and `THREADS=4` (run 41621716972, 2026-11-03):
+
+- Install: the archaic files were 150.5 GB in all, from 1.8 GB (chromosome 22) to 13.8 GB (chromosome 2). Downloading them took 2 h 32 min summed over the 22 chromosomes, from 3 to 32 minutes each, at the speed MPI-EVA served that day. Building the table took 3 h 15 min summed, from 2 to 20 minutes per chromosome. The finished table is 0.44 GB with 51,812,908 sites.
+- Step 43, one population: 19 to 88 seconds per chromosome, 18 minutes summed over the 22. IBDmix runs once per chromosome and archaic genome, `THREADS` chromosomes at a time.
+
+The same run gave HG002 a Neanderthal share of 4.13% with the 1000 Genomes ALL frequencies and 4.26% with EUR, over 51,810,955 sites of the 22 autosomes, and a Denisovan share of 0.40% and 0.42%. That is about twice what Chen et al. report with the same method and cutoffs for Europeans (about 51 Mb per person; HG002 got 115.2 Mb here), and the cause is not yet known. Until it is, compare your share with other shares from this step, not with published figures.
 
 ## Notes
 

docs/44-ancient.md

@@ -67,8 +67,10 @@ Both reports have an Ancient Populations section: the ten nearest groups, the ca
 
 ## Runtime
 
-- Setup: the 7.1 GB download, then convertf and two smartpca runs over the whole release. smartpca holds the packed genotypes of every person in memory, about 7 GB, so the setup container gets 12 GB.
-- Step 44: minutes. bcftools reads the gVCF at the panel's sites through its index, and smartpca works on about 820 people plus you; it gets 8 GB.
+Measured on public HG002 by the `ancient-footprint` job of `public-evidence.yml`, on a GitHub runner with 4 CPUs, 16 GB of RAM and `THREADS=4` (run 41621716972, 2026-11-03):
+
+- Setup: 61 minutes, the 7.1 GB download included, then convertf and two smartpca runs over the whole release. smartpca holds the packed genotypes of every person in memory. The highest anonymous memory of any setup container, convertf included, was 7.6 GiB (sampled every 2 seconds; the sampler does not record which container). The container's whole charge, file cache included, reached 11.6 GiB of its 12 GiB limit; the kernel drops that cache before it kills a task, and no swap was used. Used disk grew by 9.6 GiB at most while it ran; the installed panel keeps 1.05 GiB (1,150,106 sites on GRCh38, 819 present-day and 10,580 ancient people).
+- Step 44: 7 minutes. bcftools reads the gVCF at the panel's sites through its index, and smartpca works on the 819 people plus you; it gets 8 GB and used 1.2 GiB. The gVCF the job builds from GIAB's HG002 genotypes gave 1,133,951 of the 1,150,106 sites.
 
 ## Notes
 

docs/hardware-requirements.md

@@ -68,6 +68,7 @@ While a VEP cache installs, the tarball and the unpacked tree are on disk togeth
 | CNVpytor GC/mask files | ~0.09 GB | step 18 |
 | pypgx bundle | ~0.4 GB | step 32 |
 | Somatic resources (gnomAD AF-only VCF, panel of normals) | ~3 GB | step 29 |
+| AADR ancient panel, 1240K (7.1 GB `.geno`; 1.05 GiB kept, 9.6 GiB of extra disk at the peak while it builds) | ~7.1 GB (the `.geno`) plus three small files | step 44, `setup.sh --ancient-panel` |
 
 ### Annotation Databases (Optional, for Steps 30-31)
 
@@ -125,10 +126,12 @@ Every step runs in a Docker container with a hard `--memory` limit. The two entr
 | 25 (PRS): scoring files written as pgsc_calc reads them | no limit | `PRS_PREPARE` 1 / 16 GB |
 | 25 (PRS): score positions genotyped from the gVCF | 2 / 8 GB | `PRS_SCORE_SITES` 2 / 4 GB |
 | 25 (PRS): pgsc_calc, each of its tasks | `PGSC_MAX_MEMORY` (three quarters of the RAM) | `PRS` 4 / 24 GB, given to pgsc_calc as `--max_memory` |
+| 44 (ancient populations): smartpca | `THREADS` (8) / 8 GB | `ANCIENT_PCA` 4 / 8 GB |
+| `setup.sh --ancient-panel`: smartpca over the whole AADR release | `THREADS` (8) / 12 GB | not a pipeline task |
 | 50 (MEI): MEGAnE | `THREADS` (8) / 16 GB | `MEI` 8 / 16 GB |
 | 50 (MEI): AnnotSV on the insertions and the calls table | 4 / 8 GB for AnnotSV, 1 / 2 GB for the table | `MEI_REPORT` 4 / 8 GB with `--annotsv_annotations`, 1 / 2 GB without |
 
-These are limits, not measured peaks. The peaks measured for each task in the end test are in the [runtime table](#runtime-per-step): on 30-40x genomes the alignment task (minimap2 and samtools together) peaked at 18.3-18.4 GB and DeepVariant at 15.9-16.5 GB. For step 25, on the nine scores measured before the four small lipid and blood-pressure scores joined the list, `PRS_PREPARE` peaked at 11.5 GB and pgsc_calc's first task, run on one thread, needed about 14 to 18 GB on a GitHub runner without a genome, and 12.7 to 13.5 GiB on the three end-test genomes ([step 25, Memory](25-prs.md#memory)). Your own run records pgsc_calc's task peaks in `prs/pgsc_calc/pipeline_info/execution_trace.txt`; the `PRS` row of the pipeline's trace is only the launcher.
+These are limits, not measured peaks. The two step 44 rows were measured on a 4-CPU GitHub runner on public HG002 ([step 44, Runtime](44-ancient.md#runtime)): the setup's containers held 7.6 GiB of anonymous memory at most (the sampler does not record which container, convertf or smartpca), and the largest container's whole charge, file cache included, reached 11.6 of its 12 GiB; step 44's smartpca used 1.2 GiB. The peaks measured for each task in the end test are in the [runtime table](#runtime-per-step): on 30-40x genomes the alignment task (minimap2 and samtools together) peaked at 18.3-18.4 GB and DeepVariant at 15.9-16.5 GB. For step 25, on the nine scores measured before the four small lipid and blood-pressure scores joined the list, `PRS_PREPARE` peaked at 11.5 GB and pgsc_calc's first task, run on one thread, needed about 14 to 18 GB on a GitHub runner without a genome, and 12.7 to 13.5 GiB on the three end-test genomes ([step 25, Memory](25-prs.md#memory)). Your own run records pgsc_calc's task peaks in `prs/pgsc_calc/pipeline_info/execution_trace.txt`; the `PRS` row of the pipeline's trace is only the launcher.
 
 **Minimum system RAM:** 16 GB. `run-all.sh` passes the machine's RAM as `--max_memory`, so on a 16 GB machine the 32 GB tasks run capped at about 15 GB. A task that needs more is killed (exit 137); its one retry asks for double, and the cap cuts that back to the same 15 GB. With the step scripts, lower DeepVariant's shards and memory with `THREADS` and `DV_MEM`, for example `THREADS=4 DV_MEM=12g ./scripts/03-deepvariant.sh <sample> [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.

docs/interpreting-results.md

@@ -657,6 +657,12 @@ One SMN1 copy suggests SMA carrier status; two do not rule it out (two copies on
 
 Expect many calls: a genome carries about 1,200 such insertions, nearly all common and harmless. Read the catalog first. **Seen** means an insertion of the same element was called inside the site's window; one copy of the FKTN or BBS1 insertion is carrier status for a recessive condition, and the TAF1 insertion on a male's single X is what causes X-linked dystonia-parkinsonism. **Not seen** is not proof of absence: short reads find an SVA's edges, not its body. **Not callable** means too few reads there (under 10x) to tell. Confirm anything that matters in a clinical laboratory. [Step 50](50-mei.md) explains the columns.
 
+## Ancient Populations (Step 44, opt-in)
+
+**Where to look:** `${SAMPLE}/ancient/${SAMPLE}_ancient.tsv`, the 20 ancient groups nearest to you on a PCA of present-day West Eurasians, with their dates, places and number of people (`n`), and `${SAMPLE}/ancient/${SAMPLE}_ancient_pca.svg`, the plot with you marked.
+
+Closeness on the PCA is genetic similarity, not descent: being near an ancient group does not mean you descend from those people. A group with `n` of 1 or 2 has a noisy position. The axes describe West Eurasia only, so a genome with most of its ancestry from elsewhere lands at the edge and its nearest groups say little. `sites_used` in `${SAMPLE}_ancient_info.tsv` tells how many panel sites your gVCF covered; a whole genome covers nearly all of them (1,133,951 of the 1,150,106 that lift to GRCh38, for a gVCF built from the public GIAB HG002 genotypes). [Step 44](44-ancient.md) explains the columns.
+
 ## Polygenic Risk Scores (Step 25)
 
 **Where to look:** `${SAMPLE}/prs/${SAMPLE}_prs_summary.tsv`, one row per condition, and the PRS section of the report.

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/               # <sample>_chip.tsv, <sample>_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 | `<sample>_traits.tsv` and `<sample>_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. `<sample>_genotype_rules.tsv` with WhatsHap's phase from a BAM, byte for byte, in case `genotype-rules-1-phase` |
 | archaic | `<sample>_archaic.tsv`, `<sample>_archaic_segments.tsv` and `<sample>_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 `<sample>_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/<sample>/`, the pipeline under `<outdir>/<sample>/`, mapped in the script).