diff --git a/bin/archaic_report.py b/bin/archaic_report.py index e37c5f72..2ad55061 100755 --- a/bin/archaic_report.py +++ b/bin/archaic_report.py @@ -23,16 +23,20 @@ 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..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 @@ -40,7 +44,11 @@ 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..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=\n") for n in ARCHAICS: - out.write(f'##INFO=\n') + out.write(f'##INFO=\n') for p in POPS: out.write(f'##INFO=\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,7 +451,17 @@ 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 @@ -438,6 +469,13 @@ def cmd_genotypes(a): 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") diff --git a/docs/43-archaic.md b/docs/43-archaic.md index fd02411a..c4f45847 100644 --- a/docs/43-archaic.md +++ b/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 ` 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 ` 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 diff --git a/modules/local/archaic/main.nf b/modules/local/archaic/main.nf index c956c006..ec31864d 100644 --- a/modules/local/archaic/main.nf +++ b/modules/local/archaic/main.nf @@ -6,9 +6,10 @@ scripts/43-archaic.sh: ARCHAIC_GENOTYPES bin/archaic_report.py genotypes: IBDmix's genotype file of each chromosome (the archaic genomes, the sample's VCF and - the 1000 Genomes stand-ins of --archaic_pop) + the 1000 Genomes stand-ins of --archaic_pop), and each + archaic genome's mask ARCHAIC_IBDMIX IBDmix against each archaic genome, LOD 4, the options of - Chen et al. 2020 + Chen et al. 2020, with that genome's mask ARCHAIC_REPORT bin/archaic_report.py report: the segments of 50 kb or more, the share, the ten longest with their genes @@ -70,7 +71,7 @@ process ARCHAIC_IBDMIX { script: // One chromosome after another, the three archaic genomes of each; the - // same options as scripts/43-archaic.sh. + // same options and per-genome mask (-r) as scripts/43-archaic.sh. """ mkdir ibdmix cp ${gt}/spans.tsv ibdmix/spans.tsv @@ -78,7 +79,7 @@ process ARCHAIC_IBDMIX { c=\$(basename "\$f" .gt.gz) gzip -dc "\$f" > "\$c.gt" for a in Altai Vindija Denisova; do - ibdmix -g "\$c.gt" -s ${gt}/samples.txt -n "\$a" -o "ibdmix/\$c.\$a.txt" \\ + ibdmix -g "\$c.gt" -s ${gt}/samples.txt -n "\$a" -r "${gt}/\$c.\$a.mask.bed" -o "ibdmix/\$c.\$a.txt" \\ -d 4 -m 1 -a 0.01 -e 0.002 -c 2 done rm -f "\$c.gt" diff --git a/scripts/43-archaic.sh b/scripts/43-archaic.sh index 1a599972..c6bcb62a 100755 --- a/scripts/43-archaic.sh +++ b/scripts/43-archaic.sh @@ -65,7 +65,9 @@ run_in --cpus 1 --memory 4g -v "${PGP_ROOT}/bin:/pgp-bin:ro" "$PYTHON_IMAGE" \ --pop "$POP" --out-dir "${W}/gt" # One IBDmix run per chromosome and archaic genome, THREADS at a time, with the -# options of Chen et al. 2020 (their j.cell.2020.01.012-workflow/config.yaml). +# options of Chen et al. 2020 (their j.cell.2020.01.012-workflow/config.yaml), +# and, as there, the genome's mask (-r): the positions where its call is +# masked score nothing except a homozygous mismatch. echo "[2/3] IBDmix against each archaic genome..." # shellcheck disable=SC2016 # $1 and $2 belong to the inner bash run_in --cpus "$THREADS" --memory 8g "$IBDMIX_IMAGE" bash -euo pipefail -c ' @@ -73,7 +75,7 @@ run_in --cpus "$THREADS" --memory 8g "$IBDMIX_IMAGE" bash -euo pipefail -c ' c=$(basename "$1" .gt.gz) gzip -dc "$1" > "$1.txt" for a in Altai Vindija Denisova; do - ibdmix -g "$1.txt" -s "$2/gt/samples.txt" -n "$a" -o "$2/ibdmix/$c.$a.txt.part" \ + ibdmix -g "$1.txt" -s "$2/gt/samples.txt" -n "$a" -r "$2/gt/$c.$a.mask.bed" -o "$2/ibdmix/$c.$a.txt.part" \ -d 4 -m 1 -a 0.01 -e 0.002 -c 2 > "$2/ibdmix/$c.$a.log" 2>&1 \ || { tail -n 20 "$2/ibdmix/$c.$a.log"; exit 1; } mv "$2/ibdmix/$c.$a.txt.part" "$2/ibdmix/$c.$a.txt" diff --git a/tests/e2e/archaic-ancestry-1.sh b/tests/e2e/archaic-ancestry-1.sh index f85e2605..2b616f0c 100755 --- a/tests/e2e/archaic-ancestry-1.sh +++ b/tests/e2e/archaic-ancestry-1.sh @@ -4,9 +4,10 @@ # - setup.sh --archaic builds the table for GRCh37 chr22 only # (ARCHAIC_REGIONS): the three archaic chr22 files from MPI-EVA, checked # against their pinned md5 lists, and the 1000 Genomes sites, lifted to -# GRCh38, with each archaic genome read only inside its own FilterBed, -# every row inside the 1000 Genomes GRCh38 strict mask, and no row in a -# segmental duplication or within 5 bp of a 1000 Genomes indel; +# GRCh38, every row inside the 1000 Genomes GRCh38 strict mask, and each +# archaic call outside its genome's FilterBed, in a segmental duplication +# or within 5 bp of a 1000 Genomes indel marked masked (code 10 to 12), +# which step 43 hands IBDmix as its mask; # - HG002 (GIAB's v4.2.1 genotypes on chr22): the step runs, the shares are # numbers, the analysed span is most of the chromosome, and the # Neanderthal total sits inside the range Chen et al.'s own released @@ -61,12 +62,15 @@ check "only chr22 was built" bash -c '[ "$(cd "$1" && ls chr*.tsv.gz)" = chr22.t TABLE="${DATA}/chr22.tsv.gz" ROWS=$(gzip -dc "$TABLE" 2>/dev/null | awk 'END {print NR - 1}') check_ge "table rows on chr22" "$ROWS" 100000 -# The masks and the excluded regions must take sites out. Measured on this -# chromosome: 672,659 rows with neither (data version snpad2016-1kgp3-v1, run -# 41621716972), 648,493 with the masks alone (run 41847280046). A build that -# keeps 645,000 or more has lost the exclusions, and one that keeps 650,000 or -# more has lost both. -check "the masks and the exclusions dropped sites (fewer than 645,000 rows)" test "${ROWS:-0}" -lt 645000 +# The FilterBeds and the excluded regions mark calls masked (10 plus the ALT +# count) instead of dropping them, as Chen et al. keep those sites in IBDmix's +# input and mask them there. A build that lost the masks has no such code. +for i in 5 6 7; do + check_ge "rows where archaic column $((i - 4)) is masked" \ + "$(gzip -dc "$TABLE" | awk -F'\t' -v i="$i" 'NR > 1 && $i >= 10 && $i <= 12 {n++} END {print n + 0}')" 10000 +done +check "every archaic code is 0, 1, 2, 9 or 10 to 12" \ + bash -c 'gzip -dc "$1" | awk -F"\t" "NR > 1 { for (i = 5; i <= 7; i++) if (\$i !~ /^(0|1|2|9|10|11|12)\$/) bad = 1 } END {exit bad}"' _ "$TABLE" check "every row is on GRCh38 chr22" \ bash -c 'gzip -dc "$1" | awk -F"\t" "NR > 1 && \$1 != \"chr22\" {bad = 1} END {exit bad}"' _ "$TABLE" check_ge "rows where Vindija carries the ALT inside the plant window" \ diff --git a/tests/test_archaic_report.py b/tests/test_archaic_report.py index 898529a5..087409cd 100644 --- a/tests/test_archaic_report.py +++ b/tests/test_archaic_report.py @@ -3,14 +3,16 @@ Checks: 1. table: a 1000 Genomes site keeps its frequencies and gets each archaic - genome's ALT count; an archaic record with another ALT or another REF is + genome's ALT count, 10 plus it outside that genome's mask or inside an + excluded region; an archaic record with another ALT or another REF is unknown (9); a site only an archaic genome varies at gets frequency 0; a site no archaic genome was called at, or only called as reference with no modern variation, is left out; a position given twice by 1000 Genomes is left out; - 2. genotypes: the stand-in columns carry the population's ALT count; your - column is your ALT count, 0 where your VCF has no record, 9 for a - non-PASS record; a site where your VCF has another ALT is left out; a + 2. genotypes: a masked archaic call (10 to 12) goes to IBDmix as its ALT + count and its position into that genome's mask file; the stand-in + columns carry the population's ALT count; your column is your ALT + count, 0 where your VCF has no record, 9 for a non-PASS record; a site where your VCF has another ALT is left out; a chromosome your VCF has no record on is not analysed; in a split VCF an SNV keeps its call beside an insertion, a multi-allelic record or the 0/0 half of a split site, in either order, a usable record of the same @@ -25,7 +27,11 @@ 4. with IBDmix itself (IBDMIX_IMAGE, docker): a planted 600 kb stretch where the synthetic person carries the Neanderthal haplotype comes out as one Neanderthal segment over it, and a planted Denisovan one as Denisovan; - the same person without the planted stretches gets no segment at all. + the same person without the planted stretches gets no segment at all; + where the person contradicts both Neanderthals inside the planted + stretch and their calls there are masked, the segment stops at the + masked mismatches, while the same calls written as unknown (9, the + table before the fix) let it run through. Set ARCHAIC_TEST_NO_DOCKER=1 to skip part 4 on a machine without docker; CI never sets it. @@ -122,18 +128,28 @@ def table_info(out, extra=()): "--mask", f"Vindija={TMP}/vindija.bed", "--mask", f"Denisova={TMP}/denisova.bed.gz")) check("inside its mask a genome keeps its call", mi[100]["Altai"] == "2" and mi[100]["Denisova"] == "1", mi.get(100)) -check("outside its mask a genome is unknown (9), where the control read 0", - mi[100]["Vindija"] == "9" and nm[100]["Vindija"] == "0", (mi.get(100), nm.get(100))) -check("a row that only a masked-out reference call kept is dropped (200 goes, and the control kept it)", - sorted(m_rows) == [100, 500] and 200 in nm_rows, (sorted(m_rows), sorted(nm_rows))) -check("an archaic-only site inside its mask stays, with that genome's ALT", - m_rows[500][4] == "C" and mi[500]["Altai"] == "1" and mi[500]["Vindija"] == "9", (m_rows.get(500), mi.get(500))) +# Outside its mask a call is kept as masked (10 plus its ALT count), as Chen +# et al. keep the site in IBDmix's input and mask it there; the old table +# coded it 9 (unknown), which let IBDmix's segments run through it. +check("outside its mask a genome's call is masked (10 for 0/0), where the control read 0", + mi[100]["Vindija"] == "10" and nm[100]["Vindija"] == "0", (mi.get(100), nm.get(100))) +check("a row whose only call is masked stays, with that call masked (200: Altai 0/0 is 10)", + sorted(m_rows) == [100, 200, 500] and mi[200]["Altai"] == "10" and nm[200]["Altai"] == "0", + (sorted(m_rows), mi.get(200))) +check("an archaic-only site inside its mask keeps that genome's ALT, and the other genome's call is masked", + m_rows[500][4] == "C" and mi[500]["Altai"] == "1" and mi[500]["Vindija"] == "10", (m_rows.get(500), mi.get(500))) +check("unknown stays 9 inside or outside a mask", + mi.get(200, {}).get("Vindija") == "9" and mi.get(500, {}).get("Denisova") == "9", (mi.get(200), mi.get(500))) check("a mask file is read with and without the chr prefix", mi[100]["Denisova"] == "1", mi.get(100)) -# The excluded regions: a site inside one is left out whatever the genomes say. +# The excluded regions: every call at a site inside one is masked, whatever +# the genome's own mask says; the sites outside are untouched. write(f"{TMP}/excl.bed", "22\t90\t120\tsegdup\n") -ex_rows, _ = table_info(f"{TMP}/excl.vcf", ("--exclude", f"{TMP}/excl.bed")) -check("a site inside an excluded region is left out, and the rest stay", - sorted(ex_rows) == [200, 500] and 100 in nm_rows, (sorted(ex_rows), sorted(nm_rows))) +ex_rows, ex = table_info(f"{TMP}/excl.vcf", ("--exclude", f"{TMP}/excl.bed")) +check("every call at a site inside an excluded region is masked (100: 12, 10, 11), where the control read 2, 0, 1", + (ex[100]["Altai"], ex[100]["Vindija"], ex[100]["Denisova"]) == ("12", "10", "11") + and (nm[100]["Altai"], nm[100]["Vindija"], nm[100]["Denisova"]) == ("2", "0", "1"), (ex.get(100), nm.get(100))) +check("and the sites outside it keep their calls", sorted(ex_rows) == [100, 200, 500] + and (ex[200]["Altai"], ex[500]["Altai"]) == ("0", "1"), (sorted(ex_rows), ex.get(200), ex.get(500))) write(f"{TMP}/excl_other.bed", "7\t90\t120\tsegdup\n") r = run(*common, "--exclude", f"{TMP}/excl_other.bed", "--out", f"{TMP}/x.vcf") check("an exclusion list with no interval on the chromosome stops the build", @@ -232,6 +248,36 @@ def table_info(out, extra=()): (("A", "G", 0), ("A", "T", 0)), (("A", "G", 1), ("A", "T", 2)), (("A", "G", 9), ("A", "T", 0)), (("A", "G", 1), ("A", "G", 2))])) +# Masked calls: IBDmix gets the call itself and a mask per genome of the +# positions where it is masked, the runs of written rows closed by the next +# unmasked one. Codes per row: Altai, Vindija, Denisova. +mdir = f"{TMP}/table_masked" +os.makedirs(mdir) +mrows = {100: ("A\tG", "12\t0\t9"), 200: ("C\tT", "10\t11\t0"), 300: ("G\tA", "2\t0\t0"), + 400: ("T\tC", "10\t0\t0"), 500: ("A\tC", "0\t0\t0"), 600: ("C\tT", "11\t0\t0")} +with gzip.open(f"{mdir}/chr22.tsv.gz", "wt") as f: + f.write(hdr + "".join(f"chr22\t{p}\t{ra}\t{c}\t0.2\t0\t0\t0\t0\t0\n" for p, (ra, c) in mrows.items())) +write(f"{TMP}/m.vcf", "##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" + "chr22\t500\t.\tA\tT\t50\tPASS\t.\tGT\t0/1\n") # another ALT: row 500 is left out +r = run("genotypes", "--table-dir", mdir, "--vcf", f"{TMP}/m.vcf", "--pseudo", "10", "--out-dir", f"{TMP}/gtm") +check("genotypes on a table with masked calls exits 0", r.returncode == 0, r.stderr) +gm = {int(l.split("\t")[1]): l.split("\t")[4:7] for l in gzip.open(f"{TMP}/gtm/chr22.gt.gz", "rt").read().splitlines()[1:]} +check("IBDmix gets each masked call as its ALT count (12, 10, 11 are 2, 0, 1)", + gm.get(100) == ["2", "0", "9"] and gm.get(200) == ["0", "1", "0"] and gm.get(600) == ["1", "0", "0"], gm) +masks = {g: open(f"{TMP}/gtm/chr22.{g}.mask.bed").read() if os.path.exists(f"{TMP}/gtm/chr22.{g}.mask.bed") else None + for g in ar.ARCHAICS} +check("Altai's mask: 100-200, then 400-600 (the left-out row 500 does not close it)", + masks["Altai"] == "22\t99\t200\n22\t399\t600\n", masks["Altai"]) +check("Vindija's mask: 200 alone", masks["Vindija"] == "22\t199\t200\n", masks["Vindija"]) +check("Denisova's mask is written and empty", masks["Denisova"] == "", masks["Denisova"]) +check("an unmasked table writes empty masks", open(f"{TMP}/gt/chr22.Altai.mask.bed").read() == "") +bdir = f"{TMP}/table_bad" +os.makedirs(bdir) +with gzip.open(f"{bdir}/chr22.tsv.gz", "wt") as f: + f.write(hdr + "chr22\t100\tA\tG\t7\t0\t0\t0.2\t0\t0\t0\t0\t0\n") +r = run("genotypes", "--table-dir", bdir, "--vcf", f"{TMP}/m.vcf", "--out-dir", f"{TMP}/gtb") +check("a code the table never holds (7) stops the step", r.returncode == 1 and "code '7'" in r.stderr, r.stderr) + # --- 3. report --------------------------------------------------------------------- idir = f"{TMP}/ibd" os.makedirs(idir) @@ -305,12 +351,21 @@ def table_info(out, extra=()): DEN = (3200000, 3600000) -def synthetic(d, plant): +MID = (2200000, 2400000) + + +def synthetic(d, plant, mid=None): """A 4 Mb chr22 table and one person. A third of the sites carry an allele the Neanderthals share and modern people rarely do, a third one only the Denisovan carries, a third are common. The person draws both copies from the frequencies, except inside the planted stretches, where one copy - carries the archaic allele.""" + carries the archaic allele. + + mid: the middle of the planted Neanderthal stretch (MID) holds sites the + person contradicts (0/0 where both Neanderthals are 1/1), and the two + Neanderthals' calls there are outside their masks. "masked" writes them + as this table does (10 plus the call); "unknown" as the table before the + fix did (9).""" rnd = random.Random(11) os.makedirs(f"{d}/table") rows, calls = [], [] @@ -323,10 +378,15 @@ def synthetic(d, plant): else: f = rnd.uniform(0.1, 0.9) al = vi = de = sum(rnd.random() < f for _ in range(2)) - rows.append(f"chr22\t{pos}\tA\tG\t{al}\t{vi}\t{de}\t" + "\t".join([f"{f:.4f}"] * 6)) + ta, tv = al, vi # the codes the table writes g = [rnd.random() < f, rnd.random() < f] if plant and NEA[0] <= pos < NEA[1]: g[0] = al == 2 + if mid and MID[0] <= pos < MID[1]: + if kind == 0: + g = [False, False] + ta, tv = (al + ar.MASKED, vi + ar.MASKED) if mid == "masked" else (9, 9) + rows.append(f"chr22\t{pos}\tA\tG\t{ta}\t{tv}\t{de}\t" + "\t".join([f"{f:.4f}"] * 6)) if plant and DEN[0] <= pos < DEN[1]: g[0] = de == 2 if any(g): @@ -343,7 +403,7 @@ def ibdmix_run(d): os.makedirs(f"{d}/ibd") # The flags of scripts/43-archaic.sh. sh = ("set -e; gzip -dc /w/gt/chr22.gt.gz > /w/gt/chr22.gt; for a in Altai Vindija Denisova; do " - "ibdmix -g /w/gt/chr22.gt -s /w/gt/samples.txt -n $a -o /w/ibd/chr22.$a.txt " + "ibdmix -g /w/gt/chr22.gt -s /w/gt/samples.txt -n $a -r /w/gt/chr22.$a.mask.bed -o /w/ibd/chr22.$a.txt " "-d 4 -m 1 -a 0.01 -e 0.002 -c 2; done") r = subprocess.run(["docker", "run", "--rm", "--network", "none", "-u", f"{os.getuid()}:{os.getgid()}", "-v", f"{d}:/w", image("IBDMIX_IMAGE"), "sh", "-c", sh], capture_output=True, text=True) @@ -378,6 +438,32 @@ def ibdmix_run(d): else: check("without the planted stretches: no segment", body == [], body) + # The masked middle. Chen et al. keep a masked site in IBDmix's input and + # mask it there, so a homozygous mismatch still counts against a segment; + # the table before the fix wrote such calls as unknown (9), and a segment + # ran through them unopposed. "unknown" is that control: the planted + # segment crosses the contradicted middle. "masked" is this table: it + # must not. + cover_mid = {} + for mid in ("unknown", "masked"): + d = f"{TMP}/real_mid_{mid}" + os.makedirs(d) + os.chmod(d, 0o777) + synthetic(d, True, mid) + r = ibdmix_run(d) + check(f"IBDmix chain exits 0 (middle {mid})", r.returncode == 0, r.stderr + r.stdout) + if r.returncode: + continue + body = [l.split("\t") for l in open(f"{d}/out/S_archaic_segments.tsv").read().splitlines()[1:]] + print(f" segments (middle {mid}): {[(b[1], b[2], b[4]) for b in body]}") + cover_mid[mid] = sum(max(0, min(int(b[2]), MID[1]) - max(int(b[1]), MID[0])) + for b in body if b[4] == "Neanderthal") + span = MID[1] - MID[0] + check("control, calls written as unknown (9): the Neanderthal segment runs through the contradicted middle", + cover_mid.get("unknown", 0) >= 0.9 * span, cover_mid) + check("calls written as masked: IBDmix's mask keeps the mismatches, and no segment covers the middle", + "masked" in cover_mid and cover_mid["masked"] <= 0.1 * span, cover_mid) + shutil.rmtree(TMP, ignore_errors=True) print(f"{FAILS} check(s) failed" if FAILS else "all checks passed") sys.exit(1 if FAILS else 0)