fix(archaic): bring the Neanderthal share on HG002 in line with the published figure for the same method

bin/archaic_report.py

@@ -23,7 +23,13 @@
              genotypes of the Altai Neanderthal, the Vindija 33.19 Neanderthal
              and the Altai Denisovan at those sites and wherever one of the
              three carries a non-reference allele, as a sites-only VCF that
-             Picard lifts to GRCh38.
+             Picard lifts to GRCh38. Each archaic genome is read only inside
+             its own mask (--mask NAME=BED, the MPI-EVA FilterBed of that
+             genome): a site outside it is unknown (9) for that genome. A
+             site inside an excluded region (--exclude BED: segmental
+             duplications, the 5 bp around a modern indel) is left out
+             altogether. Both are the filtering Chen et al. apply before
+             IBDmix.
   genotypes  (step 43) IBDmix's genotype file of each chromosome: the three
              archaic columns, SAMPLE (yours, from the VCF) and the stand-ins;
              and spans.tsv, the span and site count of each chromosome.
@@ -106,6 +112,51 @@ def read_archaic(path, idx):
             yield int(c[0]), 1 + idx, (c[1], c[2], c[3])
 
 
+def read_mask(path, chrom):
+    """The intervals of one chromosome of a mask BED, sorted and merged.
+    A BED is 0-based and half-open, so [start, end) holds the 1-based
+    positions start + 1 to end. The chromosome is matched with and without
+    the chr prefix, because the archaic masks name it 22 and others chr22."""
+    num = str(chrom)[3:] if str(chrom).startswith("chr") else str(chrom)
+    want = {num, f"chr{num}"}
+    iv = []
+    with opener(path) as f:
+        for line in f:
+            if line.startswith(("#", "track", "browser")):
+                continue
+            c = line.split()
+            if len(c) < 3 or c[0] not in want:
+                continue
+            s, e = int(c[1]), int(c[2])
+            if e <= s:
+                continue          # an empty or reversed interval holds no position
+            iv.append((s, e))
+    iv.sort()
+    out = []
+    for s, e in iv:
+        if out and s <= out[-1][1]:
+            out[-1][1] = max(out[-1][1], e)
+        else:
+            out.append([s, e])
+    return out
+
+
+class Mask:
+    """Whether a 1-based position is inside a mask. Positions must come in
+    ascending order, as they do in the table build: the walk over the
+    intervals never goes back."""
+
+    def __init__(self, intervals):
+        self.iv = intervals
+        self.i = 0
+
+    def holds(self, pos):
+        p = pos - 1
+        while self.i < len(self.iv) and self.iv[self.i][1] <= p:
+            self.i += 1
+        return self.i < len(self.iv) and self.iv[self.i][0] <= p
+
+
 def archaic_code(rec, ref, alt):
     """The ALT count of one archaic record at a row with this REF and ALT, or 9."""
     if rec is None:
@@ -160,6 +211,25 @@ def cmd_table(a):
         streams.append(read_archaic(path, i))
     if names != list(ARCHAICS):
         die(f"--archaic must name {', '.join(ARCHAICS)} in that order, got {', '.join(names)}", 2)
+    masks = {}
+    for spec in a.mask:
+        name, _, path = spec.partition("=")
+        if name not in ARCHAICS or not path:
+            die(f"--mask takes NAME=FILE with NAME one of {', '.join(ARCHAICS)}, got {spec!r}", 2)
+        if name in masks:
+            die(f"--mask names {name} twice", 2)
+        iv = read_mask(path, a.chrom)
+        if not iv:
+            die(f"{path} holds no interval on chromosome {a.chrom}: a mask that matches nothing "
+                "would make every site unknown")
+        masks[name] = Mask(iv)
+    excluded = []
+    for path in a.exclude:
+        iv = read_mask(path, a.chrom)
+        if not iv:
+            die(f"{path} holds no interval on chromosome {a.chrom}: an empty exclusion list is a "
+                "download that went wrong, not a chromosome without segmental duplications")
+        excluded.append(Mask(iv))
     chrom = a.chrom if a.chrom.startswith("chr") else f"chr{a.chrom}"
     out = open(a.out, "w")
     out.write("##fileformat=VCFv4.2\n")
@@ -176,7 +246,17 @@ def flush():
         nonlocal rows
         if cur is None:
             return
-        r = table_row(cur, modern, arch)
+        # Every Mask is walked for every position, before any early return,
+        # because the walk only goes forward.
+        drop = [m.holds(cur) for m in excluded]
+        if any(drop):
+            return
+        # A genome is read only inside its own mask: outside it the record is
+        # dropped, so the row gets code 9 for that genome and the site cannot
+        # give it the row's ALT either.
+        inside = [None if (n in masks and not masks[n].holds(cur)) else r
+                  for n, r in zip(ARCHAICS, arch)]
+        r = table_row(cur, modern, inside)
         if r is None:
             return
         ref, alt, codes, afs = r
@@ -554,6 +634,10 @@ def main(argv=None):
     t.add_argument("--chrom", required=True)
     t.add_argument("--modern", required=True)
     t.add_argument("--archaic", action="append", default=[], help="NAME=FILE, once each for Altai, Vindija, Denisova")
+    t.add_argument("--mask", action="append", default=[],
+                   help="NAME=BED: read that archaic genome only inside this mask (its MPI-EVA FilterBed)")
+    t.add_argument("--exclude", action="append", default=[],
+                   help="BED: leave out every site inside these regions (segmental duplications, indel neighbourhoods)")
     t.add_argument("--out", required=True)
     g = sub.add_parser("genotypes", help="IBDmix's genotype files for one sample (step 43)")
     g.add_argument("--table-dir", required=True)

docs/43-archaic.md

@@ -37,15 +37,30 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags.
 
 ## Data
 
-`scripts/setup.sh --archaic <genome_dir>` builds `reference/archaic/snpad2016-1kgp3-v1/`, one `chrN.tsv.gz` per autosome on GRCh38: position, REF, ALT, the ALT count of each archaic genome (9 when it was not called there) and the 1000 Genomes ALT frequency of all and of each continental group. The rows are the 1000 Genomes biallelic SNVs where an archaic genome was called, plus every site where one of the three carries a non-reference allele. Those are the sites IBDmix's own `generate_gt` keeps from the same files.
+`scripts/setup.sh --archaic <genome_dir>` builds `reference/archaic/snpad2016-1kgp3-v2/`, one `chrN.tsv.gz` per autosome on GRCh38: position, REF, ALT, the ALT count of each archaic genome (9 when it was not called there) and the 1000 Genomes ALT frequency of all and of each continental group. The rows are the 1000 Genomes biallelic SNVs where an archaic genome was called, plus every site where one of the three carries a non-reference allele. Those are the sites IBDmix's own `generate_gt` keeps from the same files.
 
 | Source | What | Pin |
 |---|---|---|
 | MPI-EVA, https://cdna.eva.mpg.de/neandertal/Vindija/VCF/ | the snpAD genotypes of Altai, Vindija33.19 and Denisova, GRCh37, `chrN_mq25_mapab100.vcf.gz` | each file against the folder's `md5sums.txt`, whose sha256 `versions.env` pins |
+| MPI-EVA, https://cdna.eva.mpg.de/neandertal/Vindija/FilterBed/ | the mask of each archaic genome, GRCh37, `<genome>/chrN_mask.bed.gz` | the sha256 of each of the 66 files, in `versions.env` |
+| UCSC | the segmental duplications of GRCh37, `genomicSuperDups.txt.gz` | sha256 |
 | 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` |
+| 1000 Genomes | the GRCh38 strict mask, `20160622.allChr.mask.bed` (128 MB) | sha256 |
 | 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 150 GB to hash it.
+MPI-EVA publishes md5 lists, not sha256 sums, for 66 genotype 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 FilterBed files have no checksum file beside them and are about 1.5 MB each, so each one's sha256 is pinned directly.
+
+### The masks
+
+Chen et al. do not run IBDmix on every site. Each archaic genome comes with its own mask from the same project: the positions that pass its coverage, mapability 35, MQ25, tandem-repeat and indel filters. On top of that they exclude regions nobody is read at. The step does the same:
+
+- a genome's call is read only inside that genome's FilterBed; outside it the site is unknown (9) for that genome, so IBDmix scores nothing there;
+- a table row outside the 1000 Genomes GRCh38 strict mask is dropped;
+- a table row inside a segmental duplication (UCSC's GRCh37 track), or inside the window of a 1000 Genomes indel, is dropped. The window is 5 bp on each side of the indel, reaching past the longer of its REF and ALT, as the published script writes it.
+
+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.
 
 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.
 
@@ -81,7 +96,9 @@ 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 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.
+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.
+
+With the filters the step is comparable with the published figures again, and the e2e case holds it there: HG002's chr22 Neanderthal total must land inside 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. The case prints where HG002 falls inside that spread; it does not gate on the 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. The whole-genome figure is measured by a `public-evidence.yml` run, which installs the 150 GB of archaic data once per release.
 
 ## Notes
 
@@ -90,7 +107,9 @@ The same run gave HG002 a Neanderthal share of 4.13% with the 1000 Genomes ALL f
 - A segment counts when either of your two copies carries it. The share is of the analysed span, not of the whole genome.
 - The table has autosomes only: chrX and chrY are not assessed.
 - Sites the VCF has no record for are read as the reference allele. A VCF that holds only variants (no reference calls) cannot tell a reference call from a gap in coverage.
-- On the e2e fixture the case builds the table for GRCh37 chr22:25-35 Mb. It checks HG002's GIAB v4.2.1 genotypes there, a copy of HG002 with one chromosome over chr22:28-29 Mb taken from Vindija (a Neanderthal segment must cover it), and a person built from the table who carries only the majority allele (no segment).
+- The share is of the analysed span, and the span still runs from the first to the last table site of each chromosome whatever the masks remove. A segment runs from its first site to the site after its last, so it can cross masked sequence and count those bases; a wrong genotype at a site that was kept can raise the share as well. The number is an estimate in both directions, not a floor. Compare it with the published per-person figures, which come from masked data too.
+- The CpG sites Chen et al. exclude are kept here. A C to T change at a CpG site is the commonest ancient-DNA error, so a few of the archaic genotypes at those sites are wrong in both directions, and the sites are over-represented among the matches a method like this finds. Their own filter needs data this pipeline does not install, so the share carries that error; it is one of the reasons a number from this step is an estimate.
+- On the e2e fixture the case builds the table for GRCh37 chr22 and checks HG002's GIAB v4.2.1 genotypes against the published chr22 spread, a copy of HG002 with one chromosome over chr22:28-29 Mb taken from Vindija (a Neanderthal segment must cover it), a person built from the table who carries only the majority allele (a quarter of HG002's Neanderthal sequence at most), that no table row falls outside the strict mask, and that the filters took sites out.
 
 ## References
 

scripts/ci/e2e-case-times.tsv

@@ -48,7 +48,7 @@ ancient-similarity-1	shard	100	-
 annotsv-delly-depth-1-chain	shard	24	-
 annotsv-report-section-1-chain	shard	4	-
 annotsv-sv-inputs-1-manta	shard	1	-
-archaic-ancestry-1	shard	1021	-
+archaic-ancestry-1	shard	1900	-
 bash-step-08-hla	base	224	-
 bash-step-09-str	shard	1	-
 bash-step-10-telomere	shard	10	-