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

bin/archaic_report.py

@@ -23,24 +23,32 @@
              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. 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.
+             Picard lifts to GRCh38. Each archaic genome is read inside its
+             own mask (--mask NAME=BED, the MPI-EVA FilterBed of that
+             genome): a call outside it is kept as masked, code 10 plus its
+             ALT count. Every call at a site inside an excluded region
+             (--exclude BED: segmental duplications, the 5 bp around a modern
+             indel) is masked the same way. That is how Chen et al. filter:
+             the site stays in IBDmix's input and IBDmix's --mask stops it
+             scoring, except that a homozygous mismatch there still counts
+             against a segment.
   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.
+             the mask of each archaic genome (chrN.<genome>.mask.bed, the
+             positions where its call is masked); and spans.tsv, the span and
+             site count of each chromosome.
   report     (step 43) IBDmix's segments for SAMPLE against each archaic genome,
              kept at LOD 4 and 50 kb or longer (the cutoffs of Chen et al.),
              merged into Neanderthal and Denisovan segments, with the share of
              the analysed span, the ten longest segments and the protein-coding
              genes in them.
 
 Codes: an archaic or modern genotype is the number of ALT alleles (0, 1, 2), or
-9 when it is unknown, the coding IBDmix's genotype file uses.
+9 when it is unknown, the coding IBDmix's genotype file uses. In the table an
+archaic call outside that genome's mask is 10, 11 or 12 (MASKED plus the ALT
+count); `genotypes` writes it to IBDmix as 0, 1 or 2 and lists the position in
+that genome's mask file (chrN.<genome>.mask.bed), which step 43 passes to
+IBDmix with -r.
 
 Used by scripts/setup.sh, scripts/43-archaic.sh and the ARCHAIC_* processes
 (bin/ is on the task PATH); tests/test_archaic_report.py checks it. Standard
@@ -57,6 +65,7 @@
 NEANDERTHALS = ("Altai", "Vindija")
 POPS = ("ALL", "AFR", "AMR", "EAS", "EUR", "SAS")
 TABLE_HEADER = ["chrom", "pos", "ref", "alt", "Altai", "Vindija", "Denisova"] + [f"AF_{p}" for p in POPS]
+MASKED = 10           # table code: MASKED + ALT count, a call outside the genome's mask
 PSEUDO = 100          # stand-in modern people: 200 alleles, about one 1000 Genomes population
 MIN_LOD = 4.0         # Chen et al. 2020: LOD 4
 MIN_LENGTH = 50000    # and 50 kb
@@ -235,7 +244,7 @@ def cmd_table(a):
     out.write("##fileformat=VCFv4.2\n")
     out.write(f"##contig=<ID={chrom}>\n")
     for n in ARCHAICS:
-        out.write(f'##INFO=<ID={n},Number=1,Type=Integer,Description="ALT alleles in the {n} genome (9: unknown)">\n')
+        out.write(f'##INFO=<ID={n},Number=1,Type=Integer,Description="ALT alleles in the {n} genome (9: unknown; 10 to 12: 0 to 2 outside its mask)">\n')
     for p in POPS:
         out.write(f'##INFO=<ID=AF_{p},Number=1,Type=Float,Description="1000 Genomes phase 3 ALT frequency ({p})">\n')
     out.write("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n")
@@ -248,18 +257,19 @@ def flush():
             return
         # 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)
+        excl = any([m.holds(cur) for m in excluded])
+        outside = [excl or (n in masks and not masks[n].holds(cur)) for n in ARCHAICS]
+        r = table_row(cur, modern, arch)
         if r is None:
             return
         ref, alt, codes, afs = r
+        # A genome's call outside its mask, or at an excluded site, is kept but
+        # marked masked (MASKED + code). Step 43 hands IBDmix the call and a
+        # mask of those positions, as Chen et al. do: IBDmix scores nothing
+        # there except a homozygous mismatch, which still counts against a
+        # segment. Dropping the call instead (code 9) let segments run through
+        # masked sequence nothing could contradict.
+        codes = [c + MASKED if (out and c != 9) else c for c, out in zip(codes, outside)]
         info = ";".join([f"{n}={c}" for n, c in zip(ARCHAICS, codes)] + [f"AF_{p}={v}" for p, v in zip(POPS, afs)])
         out.write(f"{chrom}\t{cur}\t.\t{ref}\t{alt}\t.\t.\t{info}\n")
         rows += 1
@@ -329,6 +339,7 @@ def read_vcf(path):
 
 
 CLASH = object()
+ARCHAIC_CODES = {str(v) for v in (0, 1, 2, 9, MASKED, MASKED + 1, MASKED + 2)}
 
 
 def pick_call(prev, call):
@@ -410,10 +421,20 @@ def cmd_genotypes(a):
         n = first = last = dropped = carried = 0
         num = chrom[3:]
         path = os.path.join(a.out_dir, f"{chrom}.gt.gz")
+        # Per archaic genome, the runs of written rows whose call is masked:
+        # [first, last] positions, closed at the next unmasked row.
+        runs = {g: [] for g in ARCHAICS}
+        open_run = {g: None for g in ARCHAICS}
         with gzip.open(path, "wt", compresslevel=1) as out:
             out.write(head + "\n")
             for c in read_table(tables[chrom]):
                 pos, ref, alt = int(c[1]), c[2], c[3]
+                arch = []
+                for g, v in zip(ARCHAICS, c[4:7]):
+                    if v not in ARCHAIC_CODES:
+                        die(f"{tables[chrom]}: {chrom}:{pos} has {g} code {v!r}; the table holds 0, 1, 2, 9 "
+                            f"or {MASKED} to {MASKED + 2}")
+                    arch.append(int(v))
                 call = calls.get(pos)
                 if call is None:
                     code = 0
@@ -430,14 +451,31 @@ def cmd_genotypes(a):
                 except ValueError:
                     ac = 0
                 ac = min(max(ac, 0), two_n)
-                out.write(f"{num}\t{pos}\t{ref}\t{alt}\t{c[4]}\t{c[5]}\t{c[6]}\t{code}\t{pseudo[ac]}\n")
+                for g, v in zip(ARCHAICS, arch):
+                    if v >= MASKED:
+                        if open_run[g] is None:
+                            open_run[g] = [pos, pos]
+                        else:
+                            open_run[g][1] = pos
+                    elif open_run[g] is not None:
+                        runs[g].append(open_run[g])
+                        open_run[g] = None
+                gts = "\t".join(str(v - MASKED if v >= MASKED else v) for v in arch)
+                out.write(f"{num}\t{pos}\t{ref}\t{alt}\t{gts}\t{code}\t{pseudo[ac]}\n")
                 if n == 0:
                     first = pos
                 last = pos
                 n += 1
         if n == 0:
             os.remove(path)
             continue
+        # IBDmix's mask (-r): chromosome as in the genotype file, and a line
+        # covers the positions start + 1 to end, as in a BED.
+        for g in ARCHAICS:
+            if open_run[g] is not None:
+                runs[g].append(open_run[g])
+            with open(os.path.join(a.out_dir, f"{chrom}.{g}.mask.bed"), "w") as f:
+                f.writelines(f"{num}\t{s0 - 1}\t{e0}\n" for s0, e0 in runs[g])
         spans.append((chrom, first, last, n, carried, dropped))
     with open(os.path.join(a.out_dir, "spans.tsv"), "w") as f:
         f.write("chrom\tfirst\tlast\tsites\tsites_with_your_alt\tsites_other_allele\n")

docs/43-archaic.md

@@ -37,7 +37,7 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags.
 
 ## Data
 
-`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.
+`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, 10 plus the count when the call is outside that genome's mask, see below) 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 |
 |---|---|---|
@@ -52,15 +52,18 @@ MPI-EVA publishes md5 lists, not sha256 sums, for 66 genotype files of about 2.5
 
 ### 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:
+Chen et al. do not score 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. Their workflow keeps the masked sites in IBDmix's input and passes the mask to IBDmix (`--mask`). IBDmix then scores nothing at a masked site except a homozygous mismatch (the archaic genome 0/0 where you are 1/1, or the reverse), which still counts against a segment. 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.
+- a genome's call outside that genome's FilterBed is kept and marked masked in the table (10, 11 or 12: 10 plus its ALT count);
+- every call at a site inside a segmental duplication (UCSC's GRCh37 track), or inside the window of a 1000 Genomes indel, is marked masked the same way. 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;
+- step 43 writes each masked call to IBDmix as its ALT count, and the masked positions to a mask file per chromosome and archaic genome, which it passes to IBDmix with `-r`;
+- a table row outside the 1000 Genomes GRCh38 strict mask is dropped.
+
+The first version of the masks dropped the masked calls instead (code 9, or no row). IBDmix then had nothing to score at those sites, the mismatches included, so a segment could run through masked sequence that contradicts it. That is one reason HG002's share stayed above the published figure after the masks; see Runtime below. The strict mask still drops rows rather than marking them: that filter runs in `setup.sh` after the liftover, and moving it to the same marking is the remaining difference from the published workflow in how masked sites are handled.
 
 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 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 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.
 
@@ -98,7 +101,7 @@ Measured on public HG002 by the `archaic-hg002` jobs of `public-evidence.yml`, o
 
 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.
+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 holds it inside the published spread: 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. On this table the e2e fork run gave 1,177 kb: closer to the published figures, but still above the 95th percentile and about twice the median, so the gap is narrowed, not closed. The strict mask still drops rows in `scripts/setup.sh` (marking them instead is the remaining step), and the whole-genome total must land in 49.9 to 82.9 Mb before the share counts as in line with the published one. 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
 
@@ -107,9 +110,9 @@ With the filters the step is comparable with the published figures again, and th
 - 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.
-- 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 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 homozygous mismatch at a masked site still ends it, as in the published workflow); 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; today it is closer to them than before, not yet in line.
 - 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.
+- 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 FilterBeds and the exclusions marked calls masked.
 
 ## References