diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e966a2..bedc6f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -84,6 +84,36 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow ### Fixed +- **The H37Rv coordinate and numbering were wrong, or missing, on any reference but H37Rv.** On the + 185-sample cohort, mapped against two assemblies, the report gave 15,653 of the second reference's + 16,296 variant sites (97%) a wrong H37Rv coordinate. Four faults: + + - The variant positions of every reference were lifted once, from the first reference's sequence, + and the map was applied by position alone. `LIFT_VARIANTS` now runs once per reference, each from + its own sequence, and the report looks each variant up by its contig. Where a position cannot be + interpolated between two anchors (an indel between them, or a stretch without a shared unique + k-mer), the stretch between the anchors is aligned (`--align-gaps`). Every position of 27 + resistance genes now lifts to the coordinate minimap2 gives it, 99.8% of the report's variant sites + do, and none of the rest fits H37Rv worse than minimap2's coordinate; a position H37Rv lacks is + shown as *not in H37Rv*, and about 0.2%, in PE_PGRS-type repeats, are left without a coordinate. + - Only the first contig of a reference was read, so on a draft assembly a variant of any other + contig was lifted with the first contig's sequence: cut into ten contigs, one of the cohort's + references had half its positions misplaced. Each contig is now chained on its own, in whichever + orientation it lies, and the draft lifts as well as the complete assembly. The blind-spot mask + now reaches every contig as well. + - `--annotate_canonical`, on in the TB profile, annotated each VCF in place against the H37Rv snpEff + database: a position of another reference was read as an H37Rv one, and since the database names + H37Rv's chromosome `Chromosome`, which no mapping reference does, it annotated nothing. The + unlifted position was then shown as the H37Rv coordinate, over the lift. Each SNP is now moved to + its H37Rv position first (`bin/lift_vcf.py`: H37Rv's base as REF, the allele complemented where + the reference lies reversed, the mapping coordinate kept in `OPOS`), under the database's name for + the chromosome (`--canonical_chrom`, default `Chromosome`), and a record lends only the + coordinate it was lifted to. + - The SNP-dynamics panel compared a variant's HGVS change (`p.Ile66Met`) with the catalogue's short + form (`I66M`), so no trajectory was ever tagged with its catalogue entry. Both are now read in + one form, in H37Rv's gene name and numbering first, and a stop or frameshift matches the gene's + `LoF` entry. + - **Gene conversion was called from reads that could not carry it, and the cohort multiplied it.** On a 185-sample cohort, 1,520 of the 1,579 calls were rows the samples themselves had left `ambiguous`, promoted because another sample had called the same stretch outright. Those diff --git a/bin/lift_vcf.py b/bin/lift_vcf.py new file mode 100755 index 0000000..37442b9 --- /dev/null +++ b/bin/lift_vcf.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python3 +"""A sample's SNPs moved to the canonical reference's coordinates, to be annotated against its database. + +A variant VCF is in the coordinates of the reference the sample was mapped against. Annotated as it stands +against the canonical (H37Rv) snpEff database, each position is read as an H37Rv one: on any other reference +that names another gene and codon, and on every reference it names nothing at all, since the database calls +its chromosome 'Chromosome' and a mapping reference does not. This writes each SNP where the liftover map of +LIFT_VARIANTS puts it: + + - CHROM is the canonical contig, or --chrom when the canonical genome has one contig (the name the + database gives it), and POS the lifted position; + - REF is the canonical base there, and ALT the sample's allele, complemented where the mapping + reference lies reversed against the canonical genome; + - INFO/OPOS keeps the mapping coordinate, which is how the report pairs the record with its variant; + - a record whose allele IS the canonical base is left out: in canonical numbering nothing changed. + +Records the lift dropped, and anything but a SNP (the report reads SNPs only), are left out too. Every ANN, +which described the mapping reference, is stripped, and the records are sorted by their new position, as +tabix needs. +""" +from __future__ import annotations + +import argparse +import gzip +import sys + +COMP = {"A": "T", "C": "G", "G": "C", "T": "A"} +_OLD_ANNOTATION = ("ANN=", "LOF=", "NMD=", "EFF=") + + +def read_map(path): + """{(contig, position): (canonical contig, position, strand)} from a liftover map + (src_contig src_pos tgt_contig tgt_pos strand).""" + out = {} + try: + with open(path, encoding="utf-8", errors="replace") as fh: + for line in fh: + c = line.rstrip("\n").split("\t") + if len(c) >= 5 and c[1].isdigit() and c[3].isdigit(): + out[(c[0], int(c[1]))] = (c[2], int(c[3]), c[4]) + except OSError: + pass + return out + + +def read_fasta(path): + """[(contig, uppercase sequence)] in file order.""" + out, name, parts = [], None, [] + with open(path, encoding="utf-8", errors="replace") as fh: + for line in fh: + if line.startswith(">"): + if name is not None: + out.append((name, "".join(parts).upper())) + head = line[1:].strip() + name, parts = (head.split()[0] if head else ""), [] + elif name is not None: + parts.append(line.strip()) + if name is not None: + out.append((name, "".join(parts).upper())) + return out + + +def _open(path): + return gzip.open(path, "rt", encoding="utf-8", errors="replace") if str(path).endswith(".gz") \ + else open(path, encoding="utf-8", errors="replace") + + +def lift(lines, lmap, contigs, chrom=None): + """(header lines, record lines sorted by canonical position, counts) for the lines of a VCF.""" + seq = dict(contigs) + order = {n: i for i, (n, _) in enumerate(contigs)} + rename = chrom if (chrom and len(contigs) == 1) else None + head, recs, used, chrom_line = [], [], set(), None + n = {"snps": 0, "lifted": 0, "unlifted": 0, "canonical_base": 0, "not_snp": 0} + for line in lines: + if line.startswith("##"): + if not line.startswith(("##contig=", "##SnpEff", "##INFO=" % (rename or tc, len(seq[tc]))) + head.append('##INFO=') + head.append(chrom_line or "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO") + return head, [r for _, r in sorted(recs, key=lambda x: x[0])], n + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--vcf", required=True, help="the sample's VCF, in mapping-reference coordinates (may be .gz)") + ap.add_argument("--map", required=True, help="the liftover map (src_contig src_pos tgt_contig tgt_pos strand)") + ap.add_argument("--fasta", required=True, help="the canonical genome the map lifts onto (for its bases)") + ap.add_argument("--chrom", default=None, + help="the chromosome name of the canonical annotation database, written in place of the " + "canonical FASTA's contig name when that genome has one contig") + ap.add_argument("-o", "--output", required=True) + a = ap.parse_args(argv) + with _open(a.vcf) as fh: + head, recs, n = lift(fh, read_map(a.map), read_fasta(a.fasta), a.chrom) + with open(a.output, "w", encoding="utf-8") as out: + out.write("\n".join(head + recs) + "\n") + sys.stderr.write("[lift_vcf] %s: %d SNP(s), %d lifted, %d where the lift has no coordinate, %d carrying the " + "canonical base, %d other record(s) left out\n" + % (a.vcf, n["snps"], n["lifted"], n["unlifted"], n["canonical_base"], n["not_snp"])) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/bin/pathotypr_liftover.README.md b/bin/pathotypr_liftover.README.md index 4b91aa0..057178d 100644 --- a/bin/pathotypr_liftover.README.md +++ b/bin/pathotypr_liftover.README.md @@ -7,11 +7,35 @@ flanking-context k-mer and finds it in the target genome; the target position is ## `lift` — the recommended method (self-contained, synteny-anchored) ```bash -# lift the positions in POS (BED or 1-based list) from reference A onto reference B: -python3 pathotypr_liftover.py lift POS A.fasta B.fasta --out-map map.tsv --out-bed lifted.bed --contig B \ +# lift the positions in POS (BED, 'contigpos' or a 1-based list) from reference A onto reference B: +python3 pathotypr_liftover.py lift POS A.fasta B.fasta --out-map map.tsv --out-bed lifted.bed \ --kmer-size 21 --global-chain ``` +With `--global-chain` every contig of both FASTAs is read. Each source contig is chained on its own, so the +contigs of a draft assembly are placed wherever, and in whichever orientation, they lie in the target; the map +names both contigs and the strand (`src_contig src_pos tgt_contig tgt_pos strand`) and the BED is written +per target contig (`--contig` renames a single-contig target). A position names its contig; on a single-contig +source it may name any (a BED in NC_000962.3 coordinates read on an ancestor FASTA named otherwise), on a +multi-contig one a position without a contig, or on a contig the source lacks, is left out and counted. The +other modes read only the first contig of each FASTA. + +**`--align-gaps` (what `LIFT_VARIANTS` uses).** Interpolation only crosses an exactly colinear gap of ≤ ~2k, so a +position beside an indel, or in a stretch with no shared unique k-mer (SNP-dense, repeated), drops. With +`--align-gaps` the stretch between the two anchors is aligned instead (global, affine gaps with BWA-MEM's scores, +banded, indels left-aligned as VCF normalisation does; pure Python, `--max-align-gap` 20 kb): + +- a position is placed when a gap-free run of the alignment at ≥ `--align-min-identity` (0.9) joins it to one of + the two anchors, and its own 10 bases on each side match — or, on the side where they do not, a short indel + (≤ 20 bp) sits right there, or a long one beside which its context occurs only once in the stretch; +- a position inside a run the target lacks is written as absent (`.` as its target) when that run is the only + indel between the anchors, its flanks align at ≥ 90% and no anchor inside the gap points elsewhere; +- anything else — between two indels, beside a long indel in a repeat — is left out. + +Against minimap2 on two MTBC assemblies (L7 and A4) lifted to H37Rv, every position of 27 resistance genes lands +on minimap2's coordinate, 99.8% of random positions do, and of those that differ none fits H37Rv worse than +minimap2's; about 0.2% are left out, in PE_PGRS-type repeats. It adds ~10 s to a 4.4 Mb lift. + **`--global-chain` (recommended, what BAMpiro uses)** builds a whole-genome coordinate map: every k-mer unique in BOTH genomes is an anchor (~97 % of the MTBC genome — anchors roughly every base), chained into a collinear order (LIS), and any position is placed by **interpolating between its flanking anchors**. So a @@ -97,14 +121,16 @@ Flag style is mixed in pathotypr 0.1.0: `--tsv_pos` / `--ref_fasta` / `--tsv_gen but `--kmer-size` / `--fasta-genomes` use dashes. The classify main output is written to the `--output` name exactly (no extension); its columns are `genome k-mer k-merPOS SNPgenome SNPreference lineage`. -`map.tsv` = `src_pos tgt_pos`; `lifted.bed` = the positions collapsed into intervals in B's coords. +`map.tsv` = `src_pos tgt_pos` for `apply` and the per-position `lift`; `src_contig src_pos tgt_contig tgt_pos +strand` for `lift --global-chain`. `lifted.bed` = the positions collapsed into intervals in B's coords. ## The two BAMpiro use cases - **(A) SNP tables — variant position in H37Rv.** Lift the run's variant positions (mapping-reference - coords) **to H37Rv**: `A = mapping reference`, `B = H37Rv`. `snp_position` is the H37Rv coordinate; feed - it as the report's `pos_h37rv` (the same field the `--vcfs-h37rv` path fills). Alignment-free, replaces - the Picard/bcftools liftover. + coords) **to H37Rv**: `A = mapping reference`, `B = H37Rv`, once per reference of the run (`LIFT_VARIANTS`), + each with its own FASTA and the positions on its own contigs. The joined map feeds the report's + `pos_h37rv`, looked up by contig and position, and `bin/lift_vcf.py`, which moves each SNP to its H37Rv + position so the canonical snpEff annotation reads H37Rv's gene and codon. - **(B) Blind-spots mask on any reference.** Lift `assets/H37Rv_blindspots.bed` **from H37Rv to the run's reference**: `A = H37Rv`, `B = mapping reference`. The resulting BED (target coords) goes into the per-reference exclusion — so the mask works even when the reference is not H37Rv. diff --git a/bin/pathotypr_liftover.py b/bin/pathotypr_liftover.py index aac8087..75efa53 100644 --- a/bin/pathotypr_liftover.py +++ b/bin/pathotypr_liftover.py @@ -30,10 +30,19 @@ -- can be placed at the un-reflected coordinate (error bounded by the event size). Raise --min-identity (near- free on low-divergence references) to shrink it further; it is a hard k-mer-resolution limit, not tied to any self-similar context. + +`lift --global-chain` reads every contig of both FASTAs. Positions name their contig ('contigpos', or a BED), +and each source contig is chained on its own, so a draft assembly's contigs are placed wherever and in whichever +orientation they lie in the target. The map it writes names both contigs and the strand: +src_contig src_pos tgt_contig tgt_pos strand. A single-contig source takes every position whatever contig it +names (a BED in NC_000962.3 coordinates read on an H37Rv-colinear ancestor named otherwise). The other modes read +only the first contig of each FASTA. """ from __future__ import annotations import argparse +import bisect import sys +from collections import defaultdict def _read_first_contig(fasta): @@ -52,6 +61,93 @@ def _read_first_contig(fasta): return "".join(seq).upper() +def _read_contigs(fasta): + """[(name, uppercase sequence)] of every record of a FASTA, in file order. The name is the header up to + its first whitespace, as aligners and VCFs write it.""" + out, name, parts = [], None, [] + with open(fasta, encoding="utf-8", errors="replace") as fh: + for line in fh: + if line.startswith(">"): + if name is not None: + out.append((name, "".join(parts).upper())) + head = line[1:].strip() + name, parts = (head.split()[0] if head else ""), [] + elif name is not None: + parts.append(line.strip()) + if name is not None: + out.append((name, "".join(parts).upper())) + return out + + +class _Genome: + """A FASTA's contigs laid end to end, k bases of N apart, so one coordinate runs through the whole genome + and no k-mer spans two contigs (the rolling k-mer code resets on N). Global and contig positions are both + 1-based. A single contig is its own global coordinate, so a single-contig genome lifts exactly as before.""" + + def __init__(self, contigs, k): + self.names = [n for n, _ in contigs] + self.lengths = [len(s) for _, s in contigs] + self.index = {n: i for i, n in enumerate(self.names)} + self.starts, parts, off = [], [], 0 + for i, (_, s) in enumerate(contigs): + if i: + parts.append("N" * k) + off += k + self.starts.append(off) + parts.append(s) + off += len(s) + self.seq = "".join(parts) + + def local(self, g): + """(contig index, position) of global position g, or None in the gap between two contigs.""" + i = bisect.bisect_right(self.starts, g - 1) - 1 + if i < 0: + return None + p = g - self.starts[i] + return (i, p) if 1 <= p <= self.lengths[i] else None + + +def _named_positions(path): + """[(contig or None, 1-based position)] from a BED (contig start end ...; 0-based half-open, expanded), a + 'contigpos' list or a plain one-position-per-line list (contig None). Deduplicated, sorted. A line whose + coordinates are not digits (a header, a track line, a negative BED coordinate) is skipped, so a contig may be + called 'chromosome' without being taken for a header.""" + pos = set() + with open(path, encoding="utf-8", errors="replace") as fh: + for line in fh: + line = line.strip() + if not line or line.startswith("#"): + continue + c = line.replace(",", "\t").split("\t") + if len(c) >= 3 and c[1].isdigit() and c[2].isdigit(): + for p in range(int(c[1]) + 1, int(c[2]) + 1): + pos.add((c[0], p)) + elif len(c) >= 2 and c[1].isdigit() and not c[0].isdigit(): + pos.add((c[0], int(c[1]))) + elif c[0].isdigit(): + pos.add((None, int(c[0]))) + return sorted(pos, key=lambda x: (x[0] or "", x[1])) + + +def _resolve(named, genome): + """([(contig index, position)], positions on a contig the genome lacks, positions past a contig's end). + + A single-contig genome takes every position whatever contig it names: a BED written in NC_000962.3 + coordinates is read on an H37Rv-colinear ancestor whose record is named otherwise. A multi-contig genome + needs the name, and a position without one, or on a contig it does not have, is left out.""" + out, unknown, outside = [], 0, 0 + single = len(genome.names) == 1 + for name, p in named: + i = 0 if single else genome.index.get(name) + if i is None: + unknown += 1 + elif not 1 <= p <= genome.lengths[i]: + outside += 1 + else: + out.append((i, p)) + return out, unknown, outside + + def _positions(path): """Yield 1-based positions from either a BED (contig start end ...; 0-based half-open, expanded) or a plain one-position-per-line list. Deduplicated, sorted.""" @@ -299,22 +395,229 @@ def _chains(pairs, increasing, min_anchors, max_chains=1024): return out -def _lift_chain(args, positions): - """Whole-genome anchor-chain liftover. Anchors = k-mers unique in both genomes. The FORWARD chain (LIS) is - the collinear backbone; several REVERSE chains (successive longest decreasing runs of reverse-complement - anchors, one per inversion) cover inversions -- a single LIS holds only one, so multiple independent - inversions each need their own chain. A position is placed only when BRACKETED by two consecutive anchors - of some chain across a COLLINEAR gap (source span == target span +- indel_tol) that CONTAINS NO OTHER anchor - of any chain, AND whose interpolated coordinate passes a SEQUENCE-HOMOLOGY check (the anchors prove only the - gap ends correspond; the interior is verified by comparing the actual source/target context, so an inversion, - non-homologous filler or net-zero double-indel that a length-only test would accept is caught and dropped). - Of the verified chains that bracket it, the one with the TIGHTEST gap wins. Everything else DROPS: a - non-colinear gap, a gap straddling a rearrangement, an anchor desert, a non-homologous interior, or a - position beyond every chain's ends. No extrapolation. So a coordinate is correct or absent -- colinear/SNP - sites exact, inverted blocks mapped by the reflected coordinate, indel/RD interiors and rearrangement - boundaries dropped. RD deletions (target lost >= rd_min bp) are reported to --rd-out. A parity correction - makes the reverse chains exact for even k too.""" - import bisect +_ALN_MATCH, _ALN_MISMATCH, _ALN_OPEN, _ALN_EXTEND = 1, -4, -6, -1 # BWA-MEM's scores +_NEG = -(10 ** 9) + + +def _align(s, t, margin=16, max_cells=20_000_000): + """Global alignment of s against t with affine gaps (BWA-MEM's scores), banded around the path from corner + to corner. [(i, j)] of the alignment's columns, 0-based, None on the side of a gap; None when the band + would pass max_cells. On a tie the traceback takes the diagonal, so an indel in a repeat is placed as far + left as the scores allow, as VCF normalisation places it.""" + n, m = len(s), len(t) + lo, hi = min(0, m - n) - margin, max(0, m - n) + margin + if (n + 1) * (min(m, hi) - max(0, lo) + 1) > max_cells: + return None + o, x, mt, e_ = _ALN_OPEN, _ALN_MISMATCH, _ALN_MATCH, _ALN_EXTEND + # traceback byte per cell: bits 0-1 where H came from (0 diagonal, 1 E, 2 F); bit 2 E extended; bit 3 F extended + jb0 = min(m, hi) + prev_h = [0] + [o + e_ * j for j in range(1, jb0 + 1)] + prev_f = [_NEG] * (jb0 + 1) + tbs, starts = [bytearray([0]) + bytearray([1 | (4 if j > 1 else 0) for j in range(1, jb0 + 1)])], [0] + prev_ja = 0 + for i in range(1, n + 1): + ja, jb = max(0, i + lo), min(m, i + hi) + if ja > jb: + return None + w = jb - ja + 1 + h, f, tb = [_NEG] * w, [_NEG] * w, bytearray(w) + e = _NEG + si = s[i - 1] + np_ = len(prev_h) + for k in range(w): + j = ja + k + pk = j - prev_ja + if 0 <= pk < np_: + fo, fe = prev_h[pk] + o + e_, prev_f[pk] + e_ + fv, cf = (fe, 8) if fe >= fo else (fo, 0) + else: + fv, cf = _NEG, 0 + if k: + eo, ee = h[k - 1] + o + e_, e + e_ + e, ce = (ee, 4) if ee >= eo else (eo, 0) + else: + e, ce = _NEG, 0 + dk = pk - 1 + d = (prev_h[dk] + (mt if si == t[j - 1] and si != "N" else x)) if (j and 0 <= dk < np_) else _NEG + if d >= e and d >= fv: + h[k], src_ = d, 0 + elif e >= fv: + h[k], src_ = e, 1 + else: + h[k], src_ = fv, 2 + f[k] = fv + tb[k] = src_ | ce | cf + tbs.append(tb) + starts.append(ja) + prev_h, prev_f, prev_ja = h, f, ja + cols, i, j, state = [], n, m, 0 + while i > 0 or j > 0: + code = tbs[i][j - starts[i]] + if state == 0: + src_ = code & 3 + if src_ == 0: + cols.append((i - 1, j - 1)); i -= 1; j -= 1 + else: + state = src_ + elif state == 1: + cols.append((None, j - 1)); j -= 1 + state = 1 if code & 4 else 0 + else: + cols.append((i - 1, None)); i -= 1 + state = 2 if code & 8 else 0 + cols.reverse() + return cols + + +def _gap_positions(src, tgt, al, ar, bl, br, orient, k, min_identity, clean_flanks): + """{source position: target position, or None where the target has no counterpart} for the positions + strictly between two consecutive anchors of a chain, from an alignment of the stretch between them + (the anchors' own k-mers included, which match exactly). + + A position is placed when its column aligns two bases and it is joined to one of the two anchors by columns + without a gap, at `min_identity` or more: that anchor carries it, and any indel of the stretch lies on its + other side. Between two indels it could sit either side of either one, a unit of a repeat away, and is left + out. Against minimap2 on two MTBC assemblies, no placement under this rule was a worse fit than minimap2's: + those that differ are equally good, in repeats, or better. A position is + absent when it lies in a run the target lacks, that run is the only indel between the anchors, its flanks + align at 90% or more and `clean_flanks` (no anchor inside the gap points elsewhere). Anything else is left + out. Coordinates are global, 1-based.""" + half = k // 2 + s0, s1 = al - half, ar + half # source span, 1-based inclusive + if orient == 1: + t0, t1 = bl - half, br + half + t = tgt[t0 - 1:t1] + else: + t0, t1 = br - half, bl + half + t = "".join(_COMP.get(b, "N") for b in reversed(tgt[t0 - 1:t1])) + if s0 < 1 or t0 < 1 or s1 > len(src) or t1 > len(tgt): + return {} + s = src[s0 - 1:s1] + cols = _align(s, t) + if cols is None: + return {} + ncol = len(cols) + al_ok = [a is not None and b is not None for a, b in cols] + same = [ok and s[a] == t[b] and s[a] != "N" for ok, (a, b) in zip(al_ok, cols)] + pa, pm = [0], [0] + for ok, sm in zip(al_ok, same): + pa.append(pa[-1] + ok) + pm.append(pm[-1] + sm) + + # Where gaps fall in the alignment: a position joined to one of the anchors by columns without a gap is + # carried from that anchor, and whatever indel the stretch holds lies on its other side. A position with an + # indel between it and each anchor could sit either side of either one: in a repeat, a unit of it away. + pg = [0] + for ok in al_ok: + pg.append(pg[-1] + (not ok)) + + def carried(c): + left = pg[c] == 0 and pm[c] >= min_identity * pa[c] + right = (pg[ncol] - pg[c + 1] == 0 + and (pm[ncol] - pm[c + 1]) >= min_identity * (pa[ncol] - pa[c + 1])) + return (left or right) and beside(c, -1) and beside(c, 1) + + def beside(c, step, w=10, short=20): + """The w bases beside the position match along both sequences, or an indel sits right there: a short one + (at most `short` bases), or a long one where the position and the 2w bases on its other side occur only + once in the stretch. At the edge of a long indel in a repeat, the gap-free stretch can as well line up + with the next unit of the repeat, and then they occur again.""" + a, b = cols[c] + n = m_ = 0 + for d in range(1, w + 1): + i, j = a + step * d, b + step * d + if 0 <= i < len(s) and 0 <= j < len(t): + n += 1 + m_ += s[i] == t[j] and s[i] != "N" + if n == w and m_ >= min_identity * w: + return True + k_, run, seen = c + step, 0, 0 + while 0 <= k_ < ncol and seen < w: # the nearest gap run within w columns + if al_ok[k_]: + if run: + break + seen += 1 + else: + run += 1 + k_ += step + if not run: + return False # mismatches, and no indel to account for them + return run <= short or once(a, -step, 2 * w) + + def once(a, step, w): + """The position and up to w source bases on the given side (at least 10) match the target stretch at + min_identity in one place only.""" + lo_, hi_ = (max(0, a - w), a) if step < 0 else (a, min(len(s) - 1, a + w)) + if hi_ - lo_ < 10: + return False + probe = s[lo_:hi_ + 1] + need, n = min_identity * len(probe), 0 + for x in range(len(t) - len(probe) + 1): + if sum(1 for u, v in zip(probe, t[x:x + len(probe)]) if u == v) >= need: + n += 1 + if n > 1: + return False + return n == 1 + + def flank_ok(c, step): + """The 10 aligned columns next to a gap run, walking away from it, match at 90% or more.""" + n_al = n_m = 0 + while 0 <= c < ncol and n_al < 10: + if al_ok[c]: + n_al += 1 + n_m += same[c] + elif cols[c][0] is None: # a gap the source lacks right beside it: a replacement + return False + c += step + return n_al == 10 and n_m >= 9 + + def tpos(b): + return t0 + b if orient == 1 else t1 - b + + out, c = {}, 0 + while c < ncol: + a, b = cols[c] + if a is not None and b is not None: + p = s0 + a + if al < p < ar and carried(c): + out[p] = tpos(b) + c += 1 + elif a is not None: # a run of source bases the target lacks + run_end = c + while run_end < ncol and cols[run_end][1] is None and cols[run_end][0] is not None: + run_end += 1 + only_gap = pg[ncol] - (pg[run_end] - pg[c]) == 0 # the one indel between the two anchors + if clean_flanks and only_gap and flank_ok(c - 1, -1) and flank_ok(run_end, 1): + for cc in range(c, run_end): + p = s0 + cols[cc][0] + if al < p < ar: + out[p] = None + c = run_end + else: + c += 1 + return out + + +def _lift_chain_multi(args, named): + """Whole-genome anchor-chain liftover. Anchors = k-mers unique in both genomes. Each SOURCE CONTIG gets its + own chains: a FORWARD chain (LIS) as its collinear backbone and several REVERSE chains (successive longest + decreasing runs of reverse-complement anchors, one per inversion, or one for a whole contig that lies + reversed in the target) -- a single LIS holds only one, and a single genome-wide chain would keep the longest + run of contigs in the target's order and drop the rest. A position is placed only when BRACKETED by two + consecutive anchors of a chain of its contig across a COLLINEAR gap (source span == target span +- indel_tol) + that CONTAINS NO OTHER anchor of any chain and does not cross from one target contig to another, AND whose + interpolated coordinate passes a SEQUENCE-HOMOLOGY check (the anchors prove only the gap ends correspond; the + interior is verified by comparing the actual source/target context, so an inversion, non-homologous filler + or net-zero double-indel that a length-only test would accept is caught and dropped). Of the verified chains + that bracket it, the one with the TIGHTEST gap wins. Everything else DROPS: a non-colinear gap, a gap + straddling a rearrangement, an anchor desert, a non-homologous interior, or a position beyond every chain's + ends. No extrapolation. So a coordinate is correct or absent -- colinear/SNP sites exact, inverted blocks + mapped by the reflected coordinate, indel/RD interiors and rearrangement boundaries dropped. RD deletions + (target lost >= rd_min bp) are reported to --rd-out. A parity correction makes the reverse chains exact for + even k too. + + `named` is [(contig name or None, 1-based position)] on the source (see _resolve). Returns + ({(source contig, position): (target contig, position, '+' or '-')}, stats).""" k = args.kmer_size sample = max(1, args.sample) tol = max(0, args.indel_tol) @@ -325,8 +628,10 @@ def _lift_chain(args, positions): # sit there and pass a single-coordinate check), so drop it. Scales with sample (sparser anchors -> wider # legitimate colinear gaps). An explicit --max-gap overrides. mg = args.max_gap if args.max_gap is not None else 2 * k + 2 * sample + 2 - src = _read_first_contig(args.source_fasta) - tgt = _read_first_contig(args.target_fasta) + S = _Genome(_read_contigs(args.source_fasta), k) + T = _Genome(_read_contigs(args.target_fasta), k) + positions, n_unknown, n_outside = _resolve(named, S) + src, tgt = S.seq, T.seq vw = k // 2 # homology-verification half-window (~k bp) uA = _kmer_unique_map(src, k, sample) uB = _kmer_unique_map(tgt, k, sample) @@ -339,26 +644,37 @@ def _lift_chain(args, positions): if rc is not None: rev.append((a, rc + par)); in_r.add(a) amb = in_f & in_r # a source pos anchoring both strands -> drop - fa, fb = _lis(sorted((a, b) for a, b in fwd if a not in amb), True) - # A real inversion is a DENSE, CONTIGUOUS block of reverse anchors (~1 per `sample` bp of its span); a - # chain assembled from scattered spurious reverse anchors (a k-mer whose revcomp happens to be unique - # elsewhere) is SPARSE. Keep only dense chains, else a coincidental spurious pair could place a repeat- - # desert position at an unrelated coordinate. Density = anchors / (span / sample). - rev_chains = [(ca, cb) for ca, cb in _chains(sorted((a, b) for a, b in rev if a not in amb), False, max(2, args.min_chain)) - if len(ca) * sample >= args.min_density * (ca[-1] - ca[0] + 1)] - chains = ([(fa, fb, 1)] if fa else []) + [(ca, cb, -1) for ca, cb in rev_chains] + single = len(S.names) == 1 + by_f, by_r = defaultdict(list), defaultdict(list) + for pairs, by in ((fwd, by_f), (rev, by_r)): + for a, b in pairs: + if a not in amb: + by[0 if single else S.local(a)[0]].append((a, b)) + chains, rd = {}, [] # contig index -> [(ca, cb, orient)] + n_anchor = n_inv = n_inv_anchor = 0 + for ci in range(len(S.names)): + fa, fb = _lis(sorted(by_f.get(ci, ())), True) + # A real inversion is a DENSE, CONTIGUOUS block of reverse anchors (~1 per `sample` bp of its span); a + # chain assembled from scattered spurious reverse anchors (a k-mer whose revcomp happens to be unique + # elsewhere) is SPARSE. Keep only dense chains, else a coincidental spurious pair could place a repeat- + # desert position at an unrelated coordinate. Density = anchors / (span / sample). + rev_chains = [(ca, cb) for ca, cb in _chains(sorted(by_r.get(ci, ())), False, max(2, args.min_chain)) + if len(ca) * sample >= args.min_density * (ca[-1] - ca[0] + 1)] + chains[ci] = ([(fa, fb, 1)] if fa else []) + [(ca, cb, -1) for ca, cb in rev_chains] + n_anchor += len(fa) + n_inv += len(rev_chains) + n_inv_anchor += sum(len(c[0]) for c in rev_chains) + for i in range(1, len(fa)): # target lost >= rd_min bp + if (fa[i] - fa[i - 1]) - (fb[i] - fb[i - 1]) >= args.rd_min: + rd.append((ci, fa[i - 1], fa[i])) allpos = sorted((in_f | in_r) - amb) # forbid set: EVERY anchored source position - rd = [] # (a_left, a_right): target lost >= rd_min bp - for i in range(1, len(fa)): - if (fa[i] - fa[i - 1]) - (fb[i] - fb[i - 1]) >= args.rd_min: - rd.append((fa[i - 1], fa[i])) - def bracket(ca, cb, p, orient): """(coord, gap) if p sits between two consecutive anchors of this chain across a COLLINEAR gap - (|a_gap - orient*b_gap| <= tol, a_gap <= max_gap) that contains NO other anchored position; else - (None, None). orient is +1 forward, -1 reverse. No extrapolation past the chain's anchored span. The - no-other-anchor rule is what stops a gap from crossing an inversion/indel/rearrangement of ANY size.""" + (|a_gap - orient*b_gap| <= tol, a_gap <= max_gap) that contains NO other anchored position and whose ends + lie on one target contig; else (None, None). orient is +1 forward, -1 reverse. No extrapolation past the + chain's anchored span. The no-other-anchor rule is what stops a gap from crossing an inversion/indel/ + rearrangement of ANY size.""" j = bisect.bisect_left(ca, p) if j < len(ca) and ca[j] == p: return cb[j], 0 @@ -370,30 +686,120 @@ def bracket(ca, cb, p, orient): return None, None # gap too long, or an indel lives in it -> drop if bisect.bisect_right(allpos, al) < bisect.bisect_left(allpos, ar): return None, None # another anchor lies in the gap -> it spans a rearrangement + if (T.local(bl) or (None,))[0] != (T.local(br) or (-1,))[0]: + return None, None # the gap runs from one target contig into another return int(round(bl + (p - al) * (br - bl) / a_gap)), a_gap - good, n_fwd, n_rev, n_drop = {}, 0, 0, 0 - for p in positions: + align = getattr(args, "align_gaps", False) + max_align = getattr(args, "max_align_gap", 20000) + gaps = {} # (contig, chain, j) -> {position: target or None} + + def aligned(ci, n_chain, ca, cb, orient, g): + """(target position or None for absent, source gap) from an alignment of the stretch between the two + consecutive anchors of this chain around g, or (False, None) when it cannot say.""" + j = bisect.bisect_left(ca, g) + if j == 0 or j >= len(ca) or ca[j] == g: + return False, None + al, ar, bl, br = ca[j - 1], ca[j], cb[j - 1], cb[j] + if ar - al > max_align or abs(br - bl) > max_align: + return False, None + if (T.local(bl) or (None,))[0] != (T.local(br) or (-1,))[0]: + return False, None # from one target contig into another + key = (ci, n_chain, j) + if key not in gaps: + clean = bisect.bisect_right(allpos, al) >= bisect.bisect_left(allpos, ar) + gaps[key] = _gap_positions(src, tgt, al, ar, bl, br, orient, k, args.align_min_identity, clean) + got = gaps[key] + return (got[g], ar - al) if g in got else (False, None) + + good, n_fwd, n_rev, n_drop, n_aligned, n_absent = {}, 0, 0, n_outside, 0, 0 + for ci, p in positions: + g = S.starts[ci] + p best, best_gap, best_orient = None, None, 0 - for ca, cb, orient in chains: - c, g = bracket(ca, cb, p, orient) - # verify homology at the placed coordinate (skip the exact-anchor hit g==0: a shared unique k-mer + for ca, cb, orient in chains.get(ci, ()): + c, gp = bracket(ca, cb, g, orient) + # verify homology at the placed coordinate (skip the exact-anchor hit gp==0: a shared unique k-mer # is homologous by construction); a length-only-colinear but non-homologous interior fails here. - if c is not None and (g == 0 or _homologous(src, tgt, p, c, orient, vw, args.min_identity)): - if best is None or g < best_gap: - best, best_gap, best_orient = c, g, orient # tightest verified bracketing gap wins - if best is None: - n_drop += 1 - elif best_orient == 1: - good[p] = best; n_fwd += 1 + if c is not None and (gp == 0 or _homologous(src, tgt, g, c, orient, vw, args.min_identity)): + if best is None or gp < best_gap: + best, best_gap, best_orient = c, gp, orient # tightest verified bracketing gap wins + by_alignment = absent = False + if best is None and align: + # Not interpolable (an indel between the anchors, or no anchor for longer than 2k): the stretch + # between the anchors is aligned instead, and the tightest chain that places the position wins. + for n_chain, (ca, cb, orient) in enumerate(chains.get(ci, ())): + c, gp = aligned(ci, n_chain, ca, cb, orient, g) + if c is False: + continue + if c is None: + absent = True + elif best is None or gp < best_gap: + best, best_gap, best_orient = c, gp, orient + by_alignment = best is not None + loc = T.local(best) if best is not None else None + if loc is None: + if absent: + good[(S.names[ci], p)] = None # the target has no counterpart: inserted here + n_absent += 1 + else: + n_drop += 1 + continue + good[(S.names[ci], p)] = (T.names[loc[0]], loc[1], "+" if best_orient == 1 else "-") + n_aligned += by_alignment + if best_orient == 1: + n_fwd += 1 else: - good[p] = best; n_rev += 1 + n_rev += 1 if args.rd_out: with open(args.rd_out, "w", encoding="utf-8") as w: - for a0, a1 in rd: - w.write("%s\t%d\t%d\tRD_deletion\n" % (args.source_contig, a0, a1)) - return good, {"fwd": n_fwd, "rev": n_rev, "drop": n_drop, "anchors": len(fa), - "inv_chains": len(rev_chains), "inv_anchors": sum(len(c[0]) for c in rev_chains), "rd": len(rd)} + for ci, a0, a1 in rd: + name = args.source_contig if (single and args.source_contig) else S.names[ci] + w.write("%s\t%d\t%d\tRD_deletion\n" % (name, a0 - S.starts[ci], a1 - S.starts[ci])) + return good, {"fwd": n_fwd, "rev": n_rev, "drop": n_drop, "unknown_contig": n_unknown, "anchors": n_anchor, + "inv_chains": n_inv, "inv_anchors": n_inv_anchor, "rd": len(rd), "aligned": n_aligned, + "absent": n_absent, "source_contigs": S.names, "target_contigs": T.names} + + +def _lift_chain(args, positions): + """_lift_chain_multi for plain positions on a single-contig source: ({position: target position}, stats).""" + good, st = _lift_chain_multi(args, [(None, p) for p in positions]) + return {p: v[1] for (_, p), v in good.items() if v is not None}, st + + +def _write_map(good, path, source_contigs): + """The contig-aware map: src_contig src_pos tgt_contig tgt_pos strand, in source order. A position the + target has no counterpart of (a stretch inserted relative to it) is written with '.' as its target.""" + order = {n: i for i, n in enumerate(source_contigs)} + with open(path, "w", encoding="utf-8") as w: + w.write("src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\n") + for sc, sp in sorted(good, key=lambda x: (order.get(x[0], len(order)), x[1])): + if good[(sc, sp)] is None: # the target has no counterpart + w.write("%s\t%d\t.\t.\t.\n" % (sc, sp)) + continue + tc, tp, strand = good[(sc, sp)] + w.write("%s\t%d\t%s\t%d\t%s\n" % (sc, sp, tc, tp, strand)) + + +def _write_bed(good, path, target_contigs, name, contig=None): + """The lifted positions collapsed into BED intervals in the target's coordinates, contig by contig. `contig` + renames the target's contig when the target has only one (what callers passing --contig have always meant).""" + by = defaultdict(set) + for v in good.values(): + if v is not None: + by[v[0]].add(v[1]) + order = {n: i for i, n in enumerate(target_contigs)} + with open(path, "w", encoding="utf-8") as w: + for tc in sorted(by, key=lambda c: order.get(c, len(order))): + label = contig if (contig and len(target_contigs) == 1) else tc + ps = sorted(by[tc]) + start = prev = ps[0] + for q in ps[1:]: + if q == prev + 1: + prev = q + else: + w.write("%s\t%d\t%d\t%s\n" % (label, start - 1, prev, name)) + start = prev = q + w.write("%s\t%d\t%d\t%s\n" % (label, start - 1, prev, name)) def cmd_lift(args): @@ -403,14 +809,23 @@ def cmd_lift(args): With --global-chain, use the whole-genome anchor-chain instead (places SNP sites and other positions by interpolation between flanking anchors, not by their own k-mer).""" if args.global_chain: - positions = _positions(args.positions) - good, st = _lift_chain(args, positions) - _write_outputs(good, args.out_map, args.out_bed, args.contig, args.name) + good, st = _lift_chain_multi(args, _named_positions(args.positions)) + if args.out_map: + _write_map(good, args.out_map, st["source_contigs"]) + if args.out_bed: + _write_bed(good, args.out_bed, st["target_contigs"], args.name, args.contig) sys.stderr.write("[liftover] lift(anchor-chain): %d fwd + %d inverted = %d lifted ; %d dropped " "(indel shadow / RD / desert / boundary) ; %d collinear anchors, %d inversion " "chain(s)/%d anchors, %d RD deletion(s)\n" - % (st["fwd"], st["rev"], len(good), st["drop"], + % (st["fwd"], st["rev"], st["fwd"] + st["rev"], st["drop"], st["anchors"], st["inv_chains"], st["inv_anchors"], st["rd"])) + if getattr(args, "align_gaps", False): + sys.stderr.write("[liftover] %d of them placed by aligning the stretch between two anchors ; %d the " + "target has no counterpart of (inserted relative to it)\n" % (st["aligned"], st["absent"])) + if st["unknown_contig"]: + sys.stderr.write("[liftover] WARN %d position(s) name a contig %s does not have, or none while it has " + "%d; left out\n" % (st["unknown_contig"], args.source_fasta, + len(st["source_contigs"]))) return import bisect k = args.kmer_size @@ -467,7 +882,7 @@ def predict(p): # expected target coord from the ne good[p] = occ_sorted[0]; n_synteny += 1 else: n_drop += 1 - _write_outputs(good, args.out_map, args.out_bed, args.contig, args.name) + _write_outputs(good, args.out_map, args.out_bed, args.contig or "target", args.name) sys.stderr.write("[liftover] lift: %d anchored + %d recovered-by-synteny = %d ; %d dropped (of %d)\n" % (n_anchor, n_synteny, len(good), n_drop, len(pkmer))) @@ -505,7 +920,9 @@ def main(): l.add_argument("target_fasta", help="the reference to lift the positions onto") l.add_argument("--out-map", default=None, help="write a src_postgt_pos table") l.add_argument("--out-bed", default=None, help="collapse lifted positions into a BED in target coords") - l.add_argument("--contig", default="target", help="contig name for --out-bed") + l.add_argument("--contig", default=None, + help="contig name for --out-bed (default: the target's own; with --global-chain it renames the " + "target's contig only when the target has a single one)") l.add_argument("--name", default="blindspot", help="BED feature name") l.add_argument("--kmer-size", type=int, default=21) l.add_argument("--max-shift", type=int, default=100, @@ -535,13 +952,24 @@ def main(): "sampling, non-homologous filler, net-zero double-indel, diverged decoy) and it drops. " "Higher = catches smaller rearrangements (a ~3 bp inversion perturbs ~2 bp/half, caught " "at 0.9 but not 0.8) at ~no coverage cost on low-divergence MTBC references") + l.add_argument("--align-gaps", action="store_true", + help="--global-chain: where a position cannot be interpolated (an indel between its anchors, or no " + "anchor for longer than --max-gap), align the stretch between the two anchors and place it " + "from the alignment, or report it as absent from the target where the target lacks that " + "stretch (written with '.' as its target)") + l.add_argument("--align-min-identity", type=float, default=0.9, + help="--align-gaps: min identity of the gap-free columns joining a position to the anchor that " + "carries it") + l.add_argument("--max-align-gap", type=int, default=20000, + help="--align-gaps: longest stretch between two anchors (bp) that is aligned") l.add_argument("--sample", type=int, default=1, help="--global-chain: keep ~1/N of k-mers as anchors (FracMinHash) -> ~N x less memory") l.add_argument("--rd-out", default=None, help="--global-chain: write detected RD deletions (target lost sequence) as a BED (source coords)") l.add_argument("--rd-min", type=int, default=50, help="--global-chain: minimum size (bp) of a target deletion to call it an RD block") - l.add_argument("--source-contig", default="source", help="contig name for the --rd-out BED") + l.add_argument("--source-contig", default=None, + help="contig name for the --rd-out BED (default: the source's own; used for a single-contig source)") l.set_defaults(func=cmd_lift) args = ap.parse_args() diff --git a/bin/qc_report.py b/bin/qc_report.py index 5cc7f2a..8d526c7 100644 --- a/bin/qc_report.py +++ b/bin/qc_report.py @@ -85,13 +85,14 @@ def build_parser(): help="Optional per-sample annotated VCFs -> per-SNP allele frequencies for the dynamics panel.") ap.add_argument("--vcfs-h37rv", nargs="*", default=[], help="Optional per-sample VCFs annotated against a canonical reference -> the canonical " - "amino-acid position shown alongside the used-reference one (matched by sample + contig:pos).") + "amino-acid change and gene shown alongside the used-reference ones. Records lifted to the " + "canonical coordinates carry the mapping coordinate in INFO/OPOS and are paired by it.") ap.add_argument("--aa2-label", default="H37Rv", help="Label for the canonical-reference amino-acid numbering shown by --vcfs-h37rv (default H37Rv).") ap.add_argument("--pos-liftover", default=None, - help="TSV 'mapping_poscanonical_pos' (e.g. pathotypr_liftover.py apply --out-map) -> the " - "reference-of-interest COORDINATE per variant, alignment-free. Fills pos_h37rv for the SNP " - "tables; an alternative to --vcfs-h37rv when the references do not share coordinates.") + help="TSV 'src_contig src_pos tgt_contig tgt_pos strand' (pathotypr_liftover.py lift " + "--global-chain --out-map, one per reference joined) -> the canonical COORDINATE of each " + "variant, looked up by contig and position. Fills pos_h37rv for the SNP tables.") ap.add_argument("--depth-profiles", nargs="*", default=[], help="Per-sample depth tables from depth_profile.py (windows, stretches without reads, " "genes) -> deletions and SNPs per callable kb along the genome (optional).") @@ -164,44 +165,83 @@ def _index_dir(d): return dmg_by_sample +def read_lift_map(path): + """{(contig, position): canonical position, or None where the canonical genome has no counterpart} from + the liftover map, or {} without one. + + LIFT_VARIANTS writes src_contig src_pos tgt_contig tgt_pos strand, one map per reference joined into + one, and '.' as the target of a position inserted relative to the canonical genome. A map in the older + two-column form (src_pos tgt_pos) names no contig and is keyed on (None, pos): it can only be applied + where the cohort has a single contig, since the same position of another contig, or of another reference, + is another site.""" + out = {} + if not path or not os.path.exists(path) or os.path.basename(path).startswith("NO_FILE"): + return out + try: + with open(path, encoding="utf-8", errors="replace") as fh: + for line in fh: + c = line.rstrip("\n").split("\t") + if len(c) >= 4 and c[1].strip().isdigit() and c[3].strip().isdigit(): + out[(c[0].strip(), int(c[1]))] = int(c[3]) + elif len(c) >= 4 and c[1].strip().isdigit() and c[3].strip() == ".": + out[(c[0].strip(), int(c[1]))] = None + elif len(c) >= 2 and c[0].strip().isdigit() and c[1].strip().isdigit(): + out[(None, int(c[0]))] = int(c[1]) + except OSError as e: + sys.stderr.write("[qc_report] WARN pos-liftover (%s): %s\n" % (path, e)) + return out + + def load_variants(args): - """{sample: {'chrom:pos': variant}} from --vcfs, with the canonical-reference amino acid - (--vcfs-h37rv) and/or the canonical coordinate (--pos-liftover) attached where available.""" + """{sample: {'chrom:pos': variant}} from --vcfs, with the canonical-reference (H37Rv) amino acid, gene + and coordinate attached where they are known. + + A canonical record is paired with a variant by the mapping coordinate it carries (OPOS, or Picard's + OriginalContig/OriginalStart), which ANNOTATE_CANONICAL writes as it lifts each record: its position is + then the canonical coordinate. A canonical record without one sits at the mapping coordinate itself (a + VCF annotated in place, right only for a reference that shares the canonical coordinates): it lends its + amino acid to the variant at the same contig and position, but not a coordinate, which it never had. + The liftover map fills the coordinate of every variant still without one, looked up by contig.""" _variants = parse_vcfs(args.vcfs) # parsed once, feeds both the dynamics panel and the SNP matrix - if args.vcfs_h37rv: # attach the reference-of-interest (H37Rv) amino-acid change AND coordinate per variant + label = args.aa2_label + if args.vcfs_h37rv: _h37 = parse_vcfs(args.vcfs_h37rv) for _s, _pm in _variants.items(): - _hs = _h37.get(_s, {}) - # Pair each used-ref variant with its canonical-reference record. Prefer the ORIGINAL used-ref - # coordinate a liftover stamped (OriginalContig/Start or OPOS) so the two references need NOT - # share coordinates; else fall back to the bare position (references that DO share H37Rv coords). - _by_opos, _by_pos = {}, {} - for _ck, _cv in _hs.items(): + _by_opos, _in_place = {}, {} + for _ck, _cv in _h37.get(_s, {}).items(): if _cv.get('opos'): _by_opos[_cv['opos']] = (_ck, _cv) - _by_pos.setdefault(_ck.rpartition(':')[2], (_ck, _cv)) + else: + _in_place[_ck] = _cv for _key, _v in _pm.items(): - _hit = _by_opos.get(_key) or _by_pos.get(_key.rpartition(':')[2]) + _hit = _by_opos.get(_key) + _cv = _hit[1] if _hit else _in_place.get(_key) + if not _cv: + continue + if _cv.get('aa'): + _v['aa_h37rv'] = _cv['aa'] + if _cv.get('gene'): + _v['gene_h37rv'] = _cv['gene'] if _hit: - _ck, _cv = _hit - if _cv.get('aa'): - _v['aa_h37rv'] = _cv['aa'] - _v['pos_h37rv'] = _ck # reference-of-interest coordinate (contig:pos), possibly != the mapping one - if args.pos_liftover: # alignment-free canonical coordinate per variant (pathotypr k-mer liftover map) - _lift = {} - try: - with open(args.pos_liftover, encoding="utf-8", errors="replace") as _fh: - for _line in _fh: - _c = _line.rstrip("\n").split("\t") - if len(_c) >= 2 and _c[0].strip().isdigit() and _c[1].strip().isdigit(): - _lift[_c[0].strip()] = _c[1].strip() - except OSError as _e: - sys.stderr.write("[qc_report] WARN pos-liftover (%s): %s\n" % (args.pos_liftover, _e)) - for _s, _pm in _variants.items(): + _v['pos_h37rv'] = "%s:%s" % (label, _hit[0].rpartition(':')[2]) + _lift = read_lift_map(args.pos_liftover) + if _lift: + contigs = {k.rpartition(':')[0] for _pm in _variants.values() for k in _pm} + legacy = len(contigs) <= 1 # a map without contigs is only unambiguous on a single contig + for _pm in _variants.values(): for _key, _v in _pm.items(): - _p = _key.rpartition(":")[2] - if not _v.get("pos_h37rv") and _p in _lift: - _v["pos_h37rv"] = "%s:%s" % (args.aa2_label, _lift[_p]) + if _v.get("pos_h37rv"): + continue + _c, _, _p = _key.rpartition(":") + if not _p.isdigit(): + continue + _k = (_c, int(_p)) + if _k not in _lift and legacy: + _k = (None, int(_p)) + if _k in _lift: + # None: inserted relative to the canonical genome, which has no coordinate to give + _t = _lift[_k] + _v["pos_h37rv"] = "%s:%s" % (label, "absent" if _t is None else _t) return _variants diff --git a/bin/qcreport/panels.py b/bin/qcreport/panels.py index 4359d8b..64aaef5 100644 --- a/bin/qcreport/panels.py +++ b/bin/qcreport/panels.py @@ -63,6 +63,7 @@ def build_dynamics(metadata, variants, sample_meta=None, emerge=0.25, fix=0.90, series.append({'pos': pos, 'gene': (meta or {}).get('gene', ''), 'eff': (meta or {}).get('eff', ''), 'imp': (meta or {}).get('imp', ''), 'alt': (meta or {}).get('alt', ''), 'aa': (meta or {}).get('aa', ''), 'aa_h37rv': (meta or {}).get('aa_h37rv', ''), + 'gene_h37rv': (meta or {}).get('gene_h37rv', ''), 'pos_h37rv': (meta or {}).get('pos_h37rv', ''), 'traj': [round(x, 4) for x in traj], 'dp': dps, 'flags': flags}) if not series: @@ -265,10 +266,11 @@ def build_snp_matrix(variants, reference='', max_sites=50000): except ValueError: continue st = sites.setdefault(key, {'contig': contig, 'pos': ipos, 'ref': v.get('ref', ''), - 'alt': set(), 'gene': '', 'eff': '', 'aa': '', 'aa_h37rv': '', 'pos_h37rv': '', 'cells': {}}) + 'alt': set(), 'gene': '', 'eff': '', 'aa': '', 'aa_h37rv': '', + 'gene_h37rv': '', 'pos_h37rv': '', 'cells': {}}) if v.get('alt'): st['alt'].add(v['alt']) - for k in ('gene', 'eff', 'aa', 'aa_h37rv', 'pos_h37rv'): + for k in ('gene', 'eff', 'aa', 'aa_h37rv', 'gene_h37rv', 'pos_h37rv'): if v.get(k) and not st[k]: st[k] = v[k] st['cells'][sidx[s]] = [v['af'], v.get('dp')] @@ -279,7 +281,7 @@ def build_snp_matrix(variants, reference='', max_sites=50000): ordered.sort(key=lambda x: (x['contig'], x['pos'])) rows = [{'contig': x['contig'], 'pos': x['pos'], 'ref': x['ref'], 'alt': ','.join(sorted(x['alt'])), 'gene': x['gene'], 'eff': x['eff'], 'aa': x['aa'], 'aa_h37rv': x['aa_h37rv'], - 'pos_h37rv': x['pos_h37rv'], 'n': len(x['cells']), 'cells': x['cells']} + 'gene_h37rv': x['gene_h37rv'], 'pos_h37rv': x['pos_h37rv'], 'n': len(x['cells']), 'cells': x['cells']} for x in ordered] return {'samples': samples, 'reference': reference, 'rows': rows, 'total_sites': len(sites), 'truncated': truncated} diff --git a/bin/report_assets/js/01_prelude.js b/bin/report_assets/js/01_prelude.js index a1278fa..f5dcd6c 100644 --- a/bin/report_assets/js/01_prelude.js +++ b/bin/report_assets/js/01_prelude.js @@ -33,8 +33,11 @@ var AA2LBL=R.aa2_label||'H37Rv'; // label for the canonical-reference amino-ac function aaDual(aa,aaH){ if(!aa) return ''; return esc(aa)+((aaH&&aaH!==aa)?(' ['+esc(AA2LBL)+' '+esc(aaH)+']'):''); } var REFNAME=R.ref_name||''; // the reference the samples were mapped against, shown in the SNP-table headers // mapping-reference position, plus the reference-of-interest (H37Rv) coordinate in brackets when the pipeline -// provides one that differs (no shared-coordinate assumption; pos2 is a 'contig:pos' string from the lifted VCF) -function refPos(pos,pos2){ var m=esc(''+pos); var p2=pos2?(''+pos2).split(':').pop():''; if(p2&&p2!==(''+pos)){ m+=' ['+esc(AA2LBL)+' '+esc(p2)+']'; } return m; } +// provides one that differs (no shared-coordinate assumption; pos2 is a 'contig:pos' string from the lifted VCF). +// 'absent' is a stretch of the mapping reference that H37Rv does not have: no coordinate, and said so. +function refPos(pos,pos2){ var m=esc(''+pos); var p2=pos2?(''+pos2).split(':').pop():''; + if(p2==='absent') return m+' [not in '+esc(AA2LBL)+']'; + if(p2&&p2!==(''+pos)){ m+=' ['+esc(AA2LBL)+' '+esc(p2)+']'; } return m; } var thr=assign({},R.thresholds); var athr=assign({},R.anc_thresholds||{}); // ancient (aDNA) threshold view function actv(s){return (s.anc&&R.n_ancient)?athr:thr;} // active threshold set for a sample diff --git a/bin/report_assets/js/10_dynamics.js b/bin/report_assets/js/10_dynamics.js index 9ddf8a1..a6d99ae 100644 --- a/bin/report_assets/js/10_dynamics.js +++ b/bin/report_assets/js/10_dynamics.js @@ -69,7 +69,21 @@ function dynSelCls(traj){ var vals=[],i; for(i=0;i the trajectory axis -> excluded) and with >1 value across series. var groupMeta={}; D.groups.forEach(function(g){ groupMeta[g.group]=g.meta||{}; }); diff --git a/conf/organism.config b/conf/organism.config index 863e438..14a4672 100644 --- a/conf/organism.config +++ b/conf/organism.config @@ -42,9 +42,13 @@ params { // --- Dual amino-acid numbering (optional) ------------------------------ // OFF by default. If your organism has a canonical reference you want a - // second amino-acid numbering against, build/point a SnpEff DB and enable: + // second amino-acid numbering against, build/point a SnpEff DB and enable. + // Each SNP is lifted onto canonical_ref before it is annotated, so point it + // at that genome's FASTA, and name its chromosome as the DB does: // annotate_canonical = true // canonical_snpeff_db = "Escherichia_coli_str_k_12_substr_mg1655" // a SnpEff genome id + // canonical_ref = "/path/to/K-12.fasta" // the DB's genome, colinear with it + // canonical_chrom = "Chromosome" // its chromosome's name in the DB ('' keeps the FASTA's names) // canonical_label = "K-12" // label shown in the report // --- MTBC-only features: leave OFF for a non-TB run -------------------- diff --git a/conf/tuberculosis.config b/conf/tuberculosis.config index 59cab80..414dee6 100644 --- a/conf/tuberculosis.config +++ b/conf/tuberculosis.config @@ -35,15 +35,18 @@ params { // --- Dual H37Rv amino-acid numbering ----------------------------------- // Re-annotate variants against the canonical H37Rv SnpEff DB so the report - // shows each protein change in BOTH the used-reference and H37Rv numbering, - // with Mycobrowser (Rv...) gene links. + // shows each protein change and gene in BOTH the used-reference and H37Rv + // numbering, with Mycobrowser (Rv...) gene links. Each SNP is lifted to its + // H37Rv position first, so this holds on any MTBC reference. annotate_canonical = true canonical_snpeff_db = "Mycobacterium_tuberculosis_h37rv" // bundled in the image + canonical_chrom = "Chromosome" // the name that DB gives H37Rv's chromosome canonical_label = "H37Rv" // --- H37Rv COORDINATE per variant (k-mer liftover) --------------------- // Fill each variant's H37Rv coordinate for the report's SNP tables, - // alignment-free — correct even when the mapping reference is not H37Rv. + // alignment-free — correct even when the mapping reference is not H37Rv, + // for every reference of the run and every contig of each. variant_liftover = true // canonical_ref = "/opt/pathotypr/reference.fasta" // MTBC-ancestor (H37Rv-colinear); the lift target diff --git a/docs/configuration.md b/docs/configuration.md index 42c8b94..2014028 100644 --- a/docs/configuration.md +++ b/docs/configuration.md @@ -62,7 +62,7 @@ table below. | Flag | Default | Turns on | | :--- | :---: | :--- | | `--run_pathotypr` | :octicons-x-circle-fill-16:{ .red } off | Alignment-free lineage + WHO drug-resistance typing ([pathotypr](pathotypr.md)) | -| `--annotate_canonical` | :octicons-x-circle-fill-16:{ .red } off | Dual **amino-acid** numbering (re-annotate vs H37Rv) + Mycobrowser gene links | +| `--annotate_canonical` | :octicons-x-circle-fill-16:{ .red } off | Dual **amino-acid** numbering (each SNP lifted to H37Rv, then annotated there) + Mycobrowser gene links | | `--variant_liftover` | :octicons-x-circle-fill-16:{ .red } off | Dual **coordinate** - the H37Rv position per variant (alignment-free k-mer liftover) | | `--mask_blindspots` | :octicons-x-circle-fill-16:{ .red } off | Mask the H37Rv Illumina blind-spots ([Zenodo 3701840](https://zenodo.org/records/3701840)) | | `--output_cram` | :octicons-x-circle-fill-16:{ .red } off | Publish alignments as CRAM instead of BAM | @@ -157,14 +157,15 @@ Two documented, ready-to-use **organism** configs live in | | `--pathotypr_markers` | `/opt/pathotypr/lineage_markers.tsv` | Zenodo lineage markers (bundled). | | | `--pathotypr_dr_markers` | `/opt/pathotypr/dr_markers.tsv` | Zenodo WHO drug-resistance markers (bundled). | | | `--pathotypr_min_alt` | `95` | Min alt-allele % for a DR call. Lower it (10-25) to catch heteroresistant / minority alleles. | -| **Dual AA numbering** | `--annotate_canonical` | `false` | Re-annotate variants against a canonical snpEff DB for dual (H37Rv) amino-acid numbering. | +| **Dual AA numbering** | `--annotate_canonical` | `false` | Annotate each SNP against a canonical snpEff DB at its lifted canonical position, for dual (H37Rv) amino-acid numbering and gene names on any reference. | | | `--canonical_snpeff_db` | `Mycobacterium_tuberculosis_h37rv` | Canonical snpEff genome for the second annotation. | +| | `--canonical_chrom` | `Chromosome` | The name that database gives its chromosome, written for a single-contig `--canonical_ref`; empty keeps the FASTA's names. | | | `--canonical_label` | `H37Rv` | How that numbering is labelled in the report. | -| **Canonical coords & masking** | `--variant_liftover` | `false` | Fill the H37Rv **coordinate** per variant via the alignment-free k-mer liftover (report SNP tables); works for any reference. | +| **Canonical coords & masking** | `--variant_liftover` | `false` | Fill the H37Rv **coordinate** per variant via the alignment-free k-mer liftover (report SNP tables); every reference of the run from its own sequence, every contig of it. | | | `--canonical_ref` | `/opt/pathotypr/reference.fasta` | FASTA the canonical (H37Rv-colinear) coordinates are defined on - the liftover / blind-spot **target**. | | | `--mask_blindspots` | `false` | Add the H37Rv Illumina blind-spots ([Zenodo 3701840](https://zenodo.org/records/3701840)) to the exclusion mask. | | | `--blindspot_bed` | `assets/H37Rv_blindspots.bed` | The blind-spots BED (H37Rv / NC_000962.3 coordinates). | -| | `--blindspot_liftover` | `true` | Lift the blind-spots onto the run reference (k-mer liftover) so the mask is correct on **any** reference. | +| | `--blindspot_liftover` | `true` | Lift the blind-spots onto the run reference, every contig of it (k-mer liftover), so the mask is correct on **any** reference. | | **QC Report** | `--make_qc_report` | `true` | Build the interactive HTML QC report + `qc_flags.tsv`. | | | `--make_snp_matrix` | `true` | Also emit the master SNP-matrix TSV. | | | `--depth_window` | `1000` | Window (bp) of each sample's depth profile, behind the deletions and SNPs per callable kb. | diff --git a/docs/outputs.md b/docs/outputs.md index 8db6398..dfb2387 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -36,7 +36,7 @@ results_bampiro/ │ ├── MP00091.LENS.ann.vcf.gz # 🎯 MAIN OUTPUT: Annotated Variants (SNPs/Indels) ├── MP00091.LENS.ann.vcf.gz.tbi # Index for the main VCF - ├── MP00091.LENS.canonical.ann.vcf.gz # 🧬 Variants re-annotated vs the canonical ref (only if --annotate_canonical) + ├── MP00091.LENS.canonical.ann.vcf.gz # 🧬 SNPs at their canonical (H37Rv) position, annotated there; INFO/OPOS = the mapping coordinate (only if --annotate_canonical) │ ├── MP00091.LENS.all.pos.vcf.gz # 🦴 BACKBONE: VCF containing ALL positions (WT + Variants) │ # (Ideal for phylogenetic supermatrices) diff --git a/docs/pathotypr.md b/docs/pathotypr.md index 8af7e14..899aae9 100644 --- a/docs/pathotypr.md +++ b/docs/pathotypr.md @@ -61,16 +61,24 @@ the bundled `rf_model.pathotypr` is currently unused and no flag exposes it.) Se ## Dual amino-acid numbering (H37Rv / Mycobrowser) -With `--annotate_canonical true`, each sample's variants are re-annotated with a -canonical snpEff database (`--canonical_snpeff_db`, default H37Rv) *in addition to* -your mapping reference. The report then shows every amino-acid change in **both** -numberings (the used reference + `--canonical_label`, default `H37Rv`) and links each -gene to its Mycobrowser locus (`Rv…`). The **amino-acid** number is exact when your -mapping reference shares (or is lifted to) H37Rv coordinates; the canonical -**coordinate** shown beside it is provided independently and alignment-free by the -k-mer liftover below, so the SNP tables carry the H37Rv position even for a non-H37Rv -reference. Turn `--annotate_canonical` off for non-MTB organisms, or point -`--canonical_snpeff_db` at another snpEff genome. +With `--annotate_canonical true`, each sample's SNPs are also annotated against a +canonical snpEff database (`--canonical_snpeff_db`, default H37Rv), *in addition to* +your mapping reference. Each SNP is first moved to its H37Rv position with the k-mer +liftover below (`bin/lift_vcf.py`): its REF becomes H37Rv's base, its allele is +complemented where your reference lies reversed against H37Rv, and a SNP whose allele is +H37Rv's own base is left out, since in H37Rv numbering nothing changed. So the gene and +the amino-acid change are H37Rv's on **any** MTBC reference, whatever it calls and +numbers its genes. The report then shows every amino-acid change in **both** numberings +(the used reference + `--canonical_label`, default `H37Rv`), links each gene to its +Mycobrowser locus (`Rv…`), and matches trajectories to the resistance catalogue in +H37Rv's gene names and numbering. + +The database names H37Rv's chromosome `Chromosome` (`--canonical_chrom`), which is the +name the lifted records are written under. For another organism, point +`--canonical_snpeff_db` at its snpEff genome, `--canonical_ref` at that genome's FASTA, +and `--canonical_chrom` at the name the database gives its chromosome (empty keeps the +FASTA's own names, for a multi-chromosome genome whose names already match). Turn +`--annotate_canonical` off for non-MTB organisms that have no canonical genome. ## Reference-agnostic coordinates (k-mer liftover) @@ -84,9 +92,31 @@ interpolating between its flanking anchors — then the placement is **verified the actual sequence**. The method is **"correct or absent"**: SNP sites, indels and inversions are followed, and any position it cannot pin (a repeat / low-complexity desert, or a rearrangement below the k-mer resolution) is *dropped* rather than -mis-mapped. It replaces the Picard/bcftools liftover for the report coordinate (the -amino-acid numbering still comes from the snpEff re-annotation above). See -`bin/pathotypr_liftover.README.md` for the algorithm, options and limits. +mis-mapped. `--annotate_canonical` runs it too. See `bin/pathotypr_liftover.README.md` +for the algorithm, options and limits. + +Every reference of the run is lifted from its own sequence (`LIFT_VARIANTS`, one task +per reference), and every contig of it: each contig of a draft assembly is chained on +its own, so it is placed wherever and in whichever orientation it lies in H37Rv. The +maps name the contig, and the report looks each variant up by its contig and position. + +Where a position cannot be interpolated — an indel between its two anchors, or no +shared unique k-mer for longer than ~2k (a SNP-dense or repeated stretch) — the pipeline +aligns the stretch between the two anchors (`--align-gaps`). The anchors pin both ends, +so the alignment cannot drift to another copy of a repeat. A position is placed when a +gap-free, ≥90%-identical run of the alignment joins it to one of the anchors and its own +context agrees; one inside a stretch H37Rv does not have (an IS copy, a region H37Rv +lost) is reported as **not in H37Rv**, which the report shows beside the variant; and +one that could sit on either side of an indel in a repeat is left without a coordinate. + +On the two references of a 185-sample cohort, checked against minimap2's alignment of +each to H37Rv: every position of 27 resistance genes lifts to minimap2's coordinate, +99.8% of 20,000 random positions and of the report's variant sites do, and not one of +those that differ fits H37Rv worse than minimap2's coordinate — they are exact ties in +identical repeat units, bases at the junction of a structural difference, or places +where minimap2 aligned another copy. About 0.2% stay without a coordinate, in PE_PGRS-type +repeats and rearranged stretches that minimap2 mostly cannot place either. Cut into ten +contigs with six of them reversed, the same reference lifts just as well. ## Blind-spot masking (H37Rv problematic sites) @@ -96,7 +126,7 @@ low-mappability positions that are unreliable to call (Zenodo record `assets/H37Rv_blindspots.bed`, NC_000962.3 coordinates) — to each reference's exclusion mask, so those sites are dropped from variant calling, consensus and the report. Because the BED is in H37Rv coordinates, `--blindspot_liftover true` (default) -lifts it onto the run reference with the **same k-mer liftover**, treating +lifts it onto the run reference, every contig of it, with the **same k-mer liftover**, treating `--canonical_ref` as the FASTA the blind-spots are defined on (default: the H37Rv-colinear MTBC ancestor bundled at `/opt/pathotypr/reference.fasta`; point it at a true H37Rv FASTA for an exact lift), so the mask is correct even when the reference is diff --git a/docs/qc-report.md b/docs/qc-report.md index 4a7aa00..4f31d34 100644 --- a/docs/qc-report.md +++ b/docs/qc-report.md @@ -112,8 +112,8 @@ these from the run's own outputs; you only get what you have. | `--mask-bed ` | Mask-region toggle in the genome-landscape panel | | `--consensus ` | Consensus-completeness panel | | `--lineage-colors ` | Canonical lineage palette (`lineage#hex`, e.g. *mycolorsTB*) | -| `--vcfs-h37rv ` | Canonical (H37Rv) **amino-acid** numbering shown beside the used-reference one | -| `--pos-liftover ` | Canonical (H37Rv) **coordinate** per variant, alignment-free (`pathotypr_liftover.py`) | +| `--vcfs-h37rv ` | Canonical (H37Rv) **amino-acid** numbering and gene shown beside the used-reference ones (records lifted to H37Rv, paired by `INFO/OPOS`) | +| `--pos-liftover ` | Canonical (H37Rv) **coordinate** per variant, alignment-free (`pathotypr_liftover.py lift --global-chain`; `src_contig src_pos tgt_contig tgt_pos strand`, looked up by contig) | | `--aa2-label ` | Label for the canonical numbering (default `H37Rv`) | ## Optional metadata diff --git a/docs/tutorials/other-organisms.md b/docs/tutorials/other-organisms.md index c4c80be..cf098a1 100644 --- a/docs/tutorials/other-organisms.md +++ b/docs/tutorials/other-organisms.md @@ -51,7 +51,7 @@ Everything else is already generic. For the full table see [Configuration](../co The `generic` profile already handles this, but it is worth understanding: - **Lineage / DR typing** (`--run_pathotypr`) is MTBC-specific and **off by default** - leave it off. See [Lineage & Drug-Resistance Typing](../pathotypr.md). -- **`--annotate_canonical`** stays **off**: H37Rv amino-acid numbering is meaningless for a non-MTBC organism. If you *do* want a canonical numbering for your species, set `--annotate_canonical true` with `--canonical_snpeff_db` (and `--canonical_label`) pointing at another SnpEff genome. +- **`--annotate_canonical`** stays **off**: H37Rv amino-acid numbering is meaningless for a non-MTBC organism. If you *do* want a canonical numbering for your species, set `--annotate_canonical true` with `--canonical_snpeff_db` (and `--canonical_label`) pointing at another SnpEff genome, `--canonical_ref` at that genome's FASTA (each SNP is lifted onto it before it is annotated) and `--canonical_chrom` at the name the database gives its chromosome (empty keeps the FASTA's names). !!! note "Only the H37Rv SnpEff DB is bundled" diff --git a/modules/annotation.nf b/modules/annotation.nf index 43d2576..5463db1 100644 --- a/modules/annotation.nf +++ b/modules/annotation.nf @@ -9,11 +9,13 @@ include { getSavePath; getSampleDir } from './utils' ==================================================================== */ process ANNOTATE_CANONICAL { - // Re-annotate a sample's variants against a CANONICAL reference snpEff DB (H37Rv by default; set - // params.canonical_snpeff_db for another organism) so the report can show the amino-acid change in - // both the used-reference and the canonical numbering. Existing ANN is stripped first so only the - // canonical annotation remains (qc_report.py reads the first ANN). Positions must line up with the - // canonical reference (true when the mapping reference shares its coordinates). + // Annotate a sample's SNPs against a CANONICAL reference snpEff DB (H37Rv by default; set + // params.canonical_snpeff_db for another organism) so the report can show the amino-acid change and the + // gene in both the used-reference and the canonical numbering. The VCF is first moved to the canonical + // coordinates with the liftover map (bin/lift_vcf.py): annotated in place, a position of any other + // reference was read as an H37Rv one, and no position was read at all, since the database names its + // chromosome params.canonical_chrom ('Chromosome') and no mapping reference does. Each lifted record keeps + // its mapping coordinate in INFO/OPOS, which is how the report pairs it with its variant. tag "AnnCanonical: ${sampleId}" publishDir path: { "${params.outdir}/${getSampleDir(sampleId, params)}" }, mode: params.publish_mode, saveAs: { filename -> getSavePath(filename, params) } cpus 2 @@ -21,21 +23,29 @@ process ANNOTATE_CANONICAL { input: tuple val(sampleId), val(refId), path(vcf_in) + path(lift_map) // LIFT_VARIANTS' maps of every reference, joined val canonical_db output: tuple val(sampleId), path("${sampleId}.${refId}.canonical.ann.vcf.gz"), emit: out script: + def chrom = params.canonical_chrom ? "--chrom ${params.canonical_chrom}" : "" """ set -euo pipefail - # Drop any existing ANN, then annotate against the canonical genome DB. refId is in the name so a - # sample mapped to more than one reference does not produce a file-name collision downstream. - ( bcftools annotate -x INFO/ANN "${vcf_in}" -Ov 2>/dev/null || zcat -f "${vcf_in}" ) > stripped.vcf - snpEff ann -v ${canonical_db} stripped.vcf 2> ${sampleId}.${refId}.canonical.snpeff.log \\ + # refId is in the name so a sample mapped to more than one reference does not produce a file-name + # collision downstream. + python3 ${projectDir}/bin/lift_vcf.py --vcf "${vcf_in}" --map "${lift_map}" --fasta "${params.canonical_ref}" \\ + ${chrom} -o lifted.vcf + snpEff ann -v ${canonical_db} lifted.vcf 2> ${sampleId}.${refId}.canonical.snpeff.log \\ | awk 'BEGIN{FS="\\t"; OFS="\\t"} /^#/{print;next} /^\\[/{next} NF>=8{print}' \\ | bgzip -c > ${sampleId}.${refId}.canonical.ann.vcf.gz tabix -f -p vcf ${sampleId}.${refId}.canonical.ann.vcf.gz 2>/dev/null || : > ${sampleId}.${refId}.canonical.ann.vcf.gz.tbi + # A chromosome name the database does not have annotates nothing, and says so only in the log. + if grep -q 'ERROR_CHROMOSOME_NOT_FOUND' ${sampleId}.${refId}.canonical.snpeff.log 2>/dev/null \\ + || zcat -f ${sampleId}.${refId}.canonical.ann.vcf.gz | grep -v '^#' | grep -q 'ERROR_CHROMOSOME_NOT_FOUND'; then + echo "WARN ${sampleId}: ${canonical_db} has no chromosome named as the lifted records; set --canonical_chrom" >&2 + fi """ stub: diff --git a/modules/reference.nf b/modules/reference.nf index e05188b..fe08b39 100644 --- a/modules/reference.nf +++ b/modules/reference.nf @@ -108,9 +108,10 @@ process PREPARE_REFERENCE { if [[ "!{params.blindspot_liftover}" == "true" && -s "!{params.canonical_ref}" ]]; then # Whole-genome anchor-chain liftover of the H37Rv blind-spots onto THIS reference, so the mask is # correct without assuming shared coordinates: every position is placed by interpolation between - # flanking unique anchors (handles SNP sites + indels; anchor-desert positions drop, never mis-map). + # flanking unique anchors (handles SNP sites + indels; anchor-desert positions drop, never mis-map), + # on whichever contig of the reference it lies, and the BED names that contig. python3 !{projectDir}/bin/pathotypr_liftover.py lift "!{params.blindspot_bed}" "!{params.canonical_ref}" reference.fa \ - --out-bed bs_lifted.bed --contig "$CONTIG" --kmer-size 21 --global-chain --sample 10 + --out-bed bs_lifted.bed --kmer-size 21 --global-chain --sample 10 awk 'BEGIN{OFS="\t"} $2 ~ /^[0-9]+$/ && $3 ~ /^[0-9]+$/ {print $1, $2+1, $3, "blindspot", ""}' bs_lifted.bed >> "$out" else # reference already shares H37Rv coordinates: append the blind-spots directly (contig rewrite + 1-based) diff --git a/modules/report.nf b/modules/report.nf index 0511725..1bcd5be 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -67,41 +67,45 @@ process COLLECT_DR { } process LIFT_VARIANTS { - // Whole-genome anchor-chain liftover of the run's variant positions onto the canonical (H37Rv) reference - // -> a mapping_posh37rv_pos map for the report's SNP tables (--pos-liftover). Alignment-free, no - // shared-coordinate assumption; every position (incl. SNP sites) is placed by interpolation between - // flanking unique anchors -> never mis-mapped. Pure Python (bin/pathotypr_liftover.py lift --global-chain). - tag "LiftVariants: ${reference}" - cpus 4 + // Whole-genome anchor-chain liftover of ONE reference's variant positions onto the canonical (H37Rv- + // colinear) reference -> contig-aware map (src_contig src_pos tgt_contig tgt_pos strand) for the report's + // SNP tables (--pos-liftover) and for ANNOTATE_CANONICAL. One task per reference: every reference is lifted + // from its own sequence, and only the positions on its own contigs. Lifting a second reference's positions + // with the first one's sequence put 96% of them at a wrong coordinate. Alignment-free, no shared-coordinate + // assumption; every position (incl. SNP sites) is placed by interpolation between flanking unique anchors, + // on each of the reference's contigs -> correct or absent. Where it cannot interpolate (an indel between two + // anchors, or no anchor for longer than 2k) the stretch between the anchors is aligned (--align-gaps), and a + // position the canonical genome lacks is written as such. Pure Python (bin/pathotypr_liftover.py lift + // --global-chain), about 25 s and 2.5 GB for a 4.4 Mb genome. + tag "LiftVariants: ${refId}" + cpus 1 memory { 6.GB * task.attempt } input: - path(vcfs) // the report's per-sample annotated VCFs (mapping-reference coordinates; may be .gz) - path(ref_fa) // the mapping reference FASTA (k-mer context source) - val(reference) - val(basename) + tuple val(refId), path(ref_fa) // the mapping reference FASTA (k-mer context source) + path(vcfs) // the report's per-sample VCFs, of every reference (may be .gz) output: - path("${basename}_pos_liftover.tsv"), emit: map + path("${refId}.pos_liftover.tsv"), emit: map script: """ set -euo pipefail - # union of variant positions across the report VCFs (mapping-reference coordinates) - zcat -f ${vcfs} | awk '!/^#/ && \$2 ~ /^[0-9]+\$/ {print \$2}' | sort -un > positions.txt + # This reference's contigs, and the variant positions on them. A VCF of another reference names none. + grep '^>' ${ref_fa} | sed 's/^>//; s/[[:space:]].*//' | sort -u > contigs.txt + zcat -f ${vcfs} | awk 'BEGIN{FS=OFS="\\t"} NR==FNR {c[\$1]=1; next} !/^#/ && (\$1 in c) && \$2 ~ /^[0-9]+\$/ {print \$1, \$2}' contigs.txt - \\ + | sort -u > positions.txt if [ -s positions.txt ]; then - # synteny-anchored k-mer liftover of the variant positions (mapping-ref coords) onto the canonical - # reference; ambiguous (repeat) positions are resolved by the surrounding unique anchors, never mis-mapped. python3 ${projectDir}/bin/pathotypr_liftover.py lift positions.txt ${ref_fa} ${params.canonical_ref} \\ - --out-map ${basename}_pos_liftover.tsv --kmer-size 21 --global-chain + --out-map ${refId}.pos_liftover.tsv --kmer-size 21 --global-chain --align-gaps else - printf 'src_pos\\ttgt_pos\\n' > ${basename}_pos_liftover.tsv + printf 'src_contig\\tsrc_pos\\ttgt_contig\\ttgt_pos\\tstrand\\n' > ${refId}.pos_liftover.tsv fi """ stub: """ - printf 'src_pos\\ttgt_pos\\n' > ${basename}_pos_liftover.tsv + printf 'src_contig\\tsrc_pos\\ttgt_contig\\ttgt_pos\\tstrand\\n' > ${refId}.pos_liftover.tsv """ } diff --git a/nextflow.config b/nextflow.config index e474f3a..fa32e53 100644 --- a/nextflow.config +++ b/nextflow.config @@ -176,26 +176,30 @@ params { // contig name. Use ONLY with an H37Rv-coordinate reference (H37Rv / the MTBC ancestor share coords). mask_blindspots = false blindspot_bed = "${projectDir}/assets/H37Rv_blindspots.bed" - // Lift the H37Rv blind-spots onto the run's reference with pathotypr classify (alignment-free k-mer - // liftover, bin/pathotypr_liftover.py) so the mask is correct even when the reference does NOT share - // H37Rv coordinates. Needs the pathotypr container. If false, the BED is appended directly (only correct - // for an H37Rv-coordinate reference). canonical_ref = the FASTA the blind-spot coords are defined on. + // Lift the H37Rv blind-spots onto the run's reference (alignment-free k-mer liftover, + // bin/pathotypr_liftover.py), onto every contig of it, so the mask is correct even when the reference does + // NOT share H37Rv coordinates. If false, the BED is appended directly (only correct for an H37Rv-coordinate + // reference). canonical_ref = the FASTA the blind-spot coords are defined on. blindspot_liftover = true canonical_ref = "/opt/pathotypr/reference.fasta" // MTBC-ancestor fasta (H37Rv-colinear); override with an H37Rv fasta if preferred - // Report: fill each variant's reference-of-interest (H37Rv) COORDINATE with a pathotypr k-mer liftover - // (alignment-free) instead of the Picard --vcfs-h37rv path. Opt-in; needs the pathotypr container. + // Report: fill each variant's reference-of-interest (H37Rv) COORDINATE with the alignment-free k-mer + // liftover (LIFT_VARIANTS): every reference of the run is lifted from its own sequence, contig by contig, + // onto canonical_ref. Opt-in; annotate_canonical runs it too. variant_liftover = false make_consensus = true annotate_legacy_vcfs = true annotate_main_vcf = true bcftools_sort_mem = "2G" // Also annotate each sample's variants against a second, CANONICAL reference snpEff DB so the report - // shows the amino-acid change in both the used-reference and the canonical numbering (qc_report.py - // --vcfs-h37rv). Defaults to H37Rv (MTB); change canonical_snpeff_db / canonical_label for another - // organism, or set annotate_canonical=false to turn the dual numbering off. Only meaningful when the - // mapping reference shares the canonical reference's coordinates. + // shows the amino-acid change and the gene in both the used-reference and the canonical numbering + // (qc_report.py --vcfs-h37rv). Each SNP is first lifted to its canonical position (LIFT_VARIANTS, onto + // canonical_ref), so any mapping reference works, not only one that shares the canonical coordinates. + // canonical_ref must be colinear with the snpEff genome, and canonical_chrom is the name that genome + // gives its chromosome. Defaults to H37Rv (MTB); change canonical_snpeff_db / canonical_ref / + // canonical_chrom / canonical_label for another organism. annotate_canonical = false // master on/off switch for the dual numbering canonical_snpeff_db = "Mycobacterium_tuberculosis_h37rv" // any snpEff genome for the organism of interest + canonical_chrom = "Chromosome" // its chromosome's name in that DB (a single-contig canonical_ref is renamed to it; empty keeps canonical_ref's names) canonical_label = "H37Rv" // how that numbering is labelled in the report // --- Output publishing (disk footprint) --- diff --git a/subworkflows/cohort_report.nf b/subworkflows/cohort_report.nf index 7137d1c..def3c5e 100644 --- a/subworkflows/cohort_report.nf +++ b/subworkflows/cohort_report.nf @@ -103,11 +103,22 @@ workflow COHORT_REPORT { def dr_report = asBool(params.run_pathotypr) ? COLLECT_DR(patho_dr_results.map { sId, f -> f }.collect().ifEmpty([]), tsv_name).dr : file("${projectDir}/assets/NO_FILE_DR") + // Canonical (H37Rv) coordinates of the variant positions: one liftover per reference, each from its + // own sequence and for the positions on its own contigs, joined into one map keyed by contig. The + // report's SNP tables read it, and so does the canonical annotation below, which needs each record's + // canonical position. Lifting every reference with the first one's sequence put 96% of another + // reference's positions at a wrong coordinate. + def lift_needed = asBool(params.variant_liftover) || asBool(params.annotate_canonical) + def pos_liftover = lift_needed + ? LIFT_VARIANTS(ref_bundle.map { rId, fa, idx, excl -> tuple(rId, fa) }, report_vcfs).map + .collectFile(name: "${tsv_name}_pos_liftover.tsv", keepHeader: true) + .first() + : file("${projectDir}/assets/NO_FILE_LIFTOVER") // Canonical-reference-annotated VCFs for the dual amino-acid numbering (off unless annotate_canonical). - // vcf_for_stats is the per-sample (sId,rId,vcf) channel; its positions match report_vcfs, so the - // report's sample+position merge finds each variant's canonical amino-acid change. + // Each record is moved to its canonical position before it is annotated and keeps its mapping + // coordinate in OPOS, which is how the report pairs it with its variant. def report_vcfs_h37rv = asBool(params.annotate_canonical) - ? ANNOTATE_CANONICAL(vcf_for_stats, params.canonical_snpeff_db).out + ? ANNOTATE_CANONICAL(vcf_for_stats, pos_liftover, params.canonical_snpeff_db).out .map { sId, vcf -> vcf }.collect().ifEmpty([]) : file("${projectDir}/assets/NO_FILE_H37RV") // Kraken2 per-sample reports (deduped) -> Taxonomic composition panel; empty when Kraken is off. @@ -116,13 +127,6 @@ workflow COHORT_REPORT { // Empty when the stage is off, and QC_REPORT only passes --gene-conversion for a non-empty // file, so the panel hides itself rather than rendering an empty table. def report_gconv = gconv_cohort.ifEmpty([]) - // Alignment-free canonical COORDINATE per variant via pathotypr (alternative to the --vcfs-h37rv path): - // lift the run's variant positions (mapping-reference coords) onto H37Rv and hand the map to the report. - def report_ref_fa = ref_bundle.filter { rId, fa, idx, excl -> rId == report_ref } - .map { rId, fa, idx, excl -> fa }.first() - def pos_liftover = asBool(params.variant_liftover) - ? LIFT_VARIANTS(report_vcfs, report_ref_fa, report_ref, tsv_name).map - : file("${projectDir}/assets/NO_FILE_LIFTOVER") QC_REPORT(summ.summary, summ.gene_burden, cons_files, report_gff, report_mask, report_meta, report_vcfs, report_vcfs_h37rv, pos_liftover, dr_report, report_gconv, report_kraken, depth_profiles, snp_dist, snp_matrix, provenance, tsv_name) diff --git a/tests/js/dynamics.test.mjs b/tests/js/dynamics.test.mjs new file mode 100644 index 0000000..51a2eb5 --- /dev/null +++ b/tests/js/dynamics.test.mjs @@ -0,0 +1,78 @@ +// Matching a trajectory's variant to the drug-resistance catalogue (10_dynamics.js), and the H37Rv coordinate +// shown beside a variant's own (01_prelude.js). +// +// The catalogue writes a protein change as I66M (or LoF for any loss of function of the gene) and names +// genes as H37Rv does. A variant carries HGVS (p.Ile66Met), in the numbering and gene names of the +// reference it was mapped against, and in H37Rv's where the pipeline lifted it there. + +import { strict as assert } from "node:assert"; +import { describe, it } from "node:test"; +import { emptyReport, loadFunctions } from "./harness.mjs"; + +const calls = [ + { gene: "atpE", mutation: "I66M", drug: "BDQ", gn: 1 }, + { gene: "katG", mutation: "c.609C>T", drug: "INH", gn: 3 }, + { gene: "pncA", mutation: "LoF", drug: "PZA", gn: 1 }, + { gene: "katG", mutation: "W191!", drug: "INH", gn: 1 }, +]; + +function load() { + return loadFunctions(["dynDRindex", "dynDRmatch", "aaShort", "aaLoF"], { ...emptyReport(), dr: { calls } }); +} + +describe("aaShort, HGVS to the catalogue's short form", () => { + it("turns three-letter codes into one", () => { + const fn = load(); + assert.equal(fn.aaShort("p.Ile66Met"), "I66M"); + assert.equal(fn.aaShort("p.Trp191*"), "W191*"); + assert.equal(fn.aaShort("p.Gln53Ter"), "Q53*"); + assert.equal(fn.aaShort("p.Lys12fs"), "", "a frameshift has no short form"); + assert.equal(fn.aaShort(""), ""); + }); + + it("knows a loss of function", () => { + const fn = load(); + assert.ok(fn.aaLoF("p.Lys12fs") && fn.aaLoF("p.Trp191*") && fn.aaLoF("p.Gln53Ter")); + assert.ok(!fn.aaLoF("p.Ile66Met")); + }); +}); + +describe("dynDRmatch", () => { + it("matches a missense change written in HGVS", () => { + const fn = load(); + const hit = fn.dynDRmatch(fn.dynDRindex(), { gene: "atpE", aa: "p.Ile66Met" }); + assert.equal(hit && hit.drug, "BDQ"); + }); + + it("matches in the H37Rv gene and numbering where the reference names and numbers its own way", () => { + const fn = load(); + const v = { gene: "E1ASM0057_01312", aa: "p.Ile67Met", gene_h37rv: "atpE", aa_h37rv: "p.Ile66Met" }; + assert.equal(fn.dynDRmatch(fn.dynDRindex(), v).drug, "BDQ"); + }); + + it("matches a stop written with ! and a loss of function by gene", () => { + const fn = load(); + const idx = fn.dynDRindex(); + assert.equal(fn.dynDRmatch(idx, { gene: "katG", aa: "p.Trp191*" }).mutation, "W191!"); + assert.equal(fn.dynDRmatch(idx, { gene: "pncA", aa: "p.Gln10fs" }).drug, "PZA"); + }); + + it("does not match another change of the same gene", () => { + const fn = load(); + assert.equal(fn.dynDRmatch(fn.dynDRindex(), { gene: "atpE", aa: "p.Ala63Val" }), null); + assert.equal(fn.dynDRmatch(fn.dynDRindex(), { gene: "atpE", aa: "" }), null); + }); +}); + +describe("refPos, the H37Rv coordinate beside the mapping one", () => { + it("shows the lifted coordinate where it differs, and nothing where it is the same", () => { + const { refPos } = loadFunctions(["refPos"], emptyReport()); + assert.match(refPos(762828, "H37Rv:761155"), /762828 .*\[H37Rv 761155\]/); + assert.equal(refPos(761155, "H37Rv:761155"), "761155"); + }); + + it("says a stretch H37Rv lacks is not in H37Rv, rather than giving it no coordinate at all", () => { + const { refPos } = loadFunctions(["refPos"], emptyReport()); + assert.match(refPos(1800, "H37Rv:absent"), /\[not in H37Rv\]/); + }); +}); diff --git a/tests/unit/test_lift_vcf.py b/tests/unit/test_lift_vcf.py new file mode 100644 index 0000000..da427c4 --- /dev/null +++ b/tests/unit/test_lift_vcf.py @@ -0,0 +1,98 @@ +"""Unit tests for bin/lift_vcf.py: a sample's SNPs moved to the canonical coordinates before annotation.""" + +from __future__ import annotations + +import argparse +import gzip + +from conftest import load_script + +lv = load_script("lift_vcf") +qc = load_script("qc_report") + +HEAD = ("##fileformat=VCFv4.2\n##contig=\n" + "##INFO=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n") + + +def _rec(chrom, pos, ref, alt, info="DP=40;ANN=T|missense_variant|MODERATE|E1_00001|x"): + return f"{chrom}\t{pos}\t.\t{ref}\t{alt}\t50\tPASS\t{info}\tGT:AD:DP\t1:2,38:40\n" + + +CANON = [("MTB_anc", "ACGTACGTAC" * 10)] # canonical base at p is "ACGTACGTAC"[(p - 1) % 10] + + +def _lift(records, lmap, contigs=CANON, chrom="Chromosome"): + return lv.lift((HEAD + "".join(records)).splitlines(True), lmap, contigs, chrom) + + +def test_a_snp_moves_to_its_canonical_position_with_the_canonical_base(): + head, recs, n = _lift([_rec("E1", 500, "C", "T")], {("E1", 500): ("MTB_anc", 42, "+")}) + + c = recs[0].split("\t") + assert c[:5] == ["Chromosome", "42", ".", "C", "T"], "CHROM renamed for the database, REF the canonical base" + assert c[7] == "OPOS=E1:500;DP=40", "the mapping coordinate kept, the old annotation stripped" + assert "##contig=" in head + assert not any(h.startswith("##INFO=" in head + + +def test_the_lifted_record_is_paired_with_its_variant_by_the_report(tmp_path): + """End to end: lift_vcf writes OPOS, the report pairs by it and takes the canonical coordinate.""" + vcf = tmp_path / "s.vcf.gz" + with gzip.open(vcf, "wt") as fh: + fh.write(HEAD + _rec("E1", 500, "C", "T")) + lmap = tmp_path / "lift.tsv" + lmap.write_text("src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\nE1\t500\tMTB_anc\t42\t+\n") + fa = tmp_path / "canon.fa" + fa.write_text(">MTB_anc\n" + CANON[0][1] + "\n") + lifted = tmp_path / "lifted.vcf" + + assert lv.main(["--vcf", str(vcf), "--map", str(lmap), "--fasta", str(fa), "--chrom", "Chromosome", + "-o", str(lifted)]) == 0 + # what snpEff adds against the canonical database + annotated = tmp_path / "canonical.vcf" + annotated.write_text(lifted.read_text().replace( + "OPOS=E1:500;DP=40", "OPOS=E1:500;DP=40;ANN=T|missense_variant|MODERATE|rpoB|g|t|t|pc|1/1|c.1A>T|p.Ser450Leu")) + original = tmp_path / "s.vcf" + original.write_text(HEAD + _rec("E1", 500, "C", "T")) + + v = qc.load_variants(argparse.Namespace(vcfs=[str(original)], vcfs_h37rv=[str(annotated)], + pos_liftover=None, aa2_label="H37Rv"))["S1"]["E1:500"] + + assert (v["gene_h37rv"], v["aa_h37rv"], v["pos_h37rv"]) == ("rpoB", "p.Ser450Leu", "H37Rv:42") diff --git a/tests/unit/test_pathotypr_liftover.py b/tests/unit/test_pathotypr_liftover.py index ccbf84d..bf0d727 100644 --- a/tests/unit/test_pathotypr_liftover.py +++ b/tests/unit/test_pathotypr_liftover.py @@ -629,10 +629,185 @@ def test_cmd_lift_global_chain_dispatches_to_lift_chain(tmp_path, capsys): min_identity=0.9, rd_min=50, rd_out=None, source_contig="A") lift.cmd_lift(args) - assert (tmp_path / "lift.map").read_text() == "src_pos\ttgt_pos\n100\t100\n400\t436\n" + assert (tmp_path / "lift.map").read_text() == ( + "src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\nA\t100\tB\t100\t+\nA\t400\tB\t436\t+\n") assert "lift(anchor-chain): 2 fwd + 0 inverted = 2 lifted" in capsys.readouterr().err +# ------------------------------------------------- several contigs (a draft assembly, a split target) + + +def _write_multi_fasta(path, records, width=60): + lines = [] + for name, seq in records: + lines.append(">%s" % name) + lines += [seq[i:i + width] for i in range(0, len(seq), width)] + path.write_text("\n".join(lines) + "\n", encoding="utf-8") + return str(path) + + +def test_named_positions_reads_contig_lists_beds_and_plain_positions(tmp_path): + p = tmp_path / "pos.txt" + p.write_text("chromosome\t12\nctg2\t5\t7\tx\n30\ntrack name=x\nchrom\tstart\tend\n") + assert lift._named_positions(str(p)) == [ + (None, 30), ("chromosome", 12), ("ctg2", 6), ("ctg2", 7)] + + +def test_a_draft_is_lifted_contig_by_contig_in_any_order_and_orientation(tmp_path): + """A draft assembly's contigs lie anywhere in the target, and some reverse-complemented: each is + chained on its own, so every one lands, with the strand it lands on.""" + genome = _random_seq(3000, seed=21) + c1, c2, c3 = genome[:1000], genome[1000:2000], genome[2000:] + src = _write_multi_fasta(tmp_path / "draft.fa", [("c3", c3), ("c1", c1), ("c2rc", _rc(c2))]) + tgt = _write_fasta(tmp_path / "ref.fa", "chr", genome) + + good, stats = lift._lift_chain_multi(_chain_args(src, tgt), [("c1", 500), ("c3", 100), ("c2rc", 300)]) + + assert good[("c1", 500)] == ("chr", 500, "+") + assert good[("c3", 100)] == ("chr", 2100, "+") + assert good[("c2rc", 300)] == ("chr", 1000 + (1000 - 300 + 1), "-") + assert stats["drop"] == 0 and stats["unknown_contig"] == 0 + + +def test_a_position_is_placed_on_the_target_contig_it_belongs_to(tmp_path): + genome = _random_seq(2000, seed=22) + src = _write_fasta(tmp_path / "src.fa", "chr", genome) + tgt = _write_multi_fasta(tmp_path / "tgt.fa", [("t1", genome[:1000]), ("t2", genome[1000:])]) + + good, _ = lift._lift_chain_multi(_chain_args(src, tgt), [(None, 400), (None, 1400)]) + + assert good[("chr", 400)] == ("t1", 400, "+") + assert good[("chr", 1400)] == ("t2", 400, "+") + + +def test_a_multi_contig_source_needs_the_contig_named(tmp_path): + """A plain position, or one on a contig the source lacks, is left out rather than read on the first + contig: that used to lift a variant of contig 2 with the sequence of contig 1.""" + genome = _random_seq(2000, seed=23) + src = _write_multi_fasta(tmp_path / "src.fa", [("a", genome[:1000]), ("b", genome[1000:])]) + tgt = _write_fasta(tmp_path / "tgt.fa", "chr", genome) + + good, stats = lift._lift_chain_multi(_chain_args(src, tgt), [(None, 500), ("zz", 500), ("b", 500)]) + + assert good == {("b", 500): ("chr", 1500, "+")} + assert stats["unknown_contig"] == 2 + + +def test_a_single_contig_source_takes_positions_whatever_contig_they_name(tmp_path): + """The H37Rv blind-spot BED names NC_000962.3 and is read on an ancestor FASTA named otherwise.""" + seq = _random_seq(600, seed=11) + src = _write_fasta(tmp_path / "anc.fa", "MTB_anc", seq) + tgt = _write_fasta(tmp_path / "ref.fa", "E1", seq) + + good, _ = lift._lift_chain_multi(_chain_args(src, tgt), [("NC_000962.3", 300)]) + + assert good == {("MTB_anc", 300): ("E1", 300, "+")} + + +def test_cmd_lift_global_chain_writes_a_bed_per_target_contig(tmp_path, capsys): + genome = _random_seq(2000, seed=22) + src = _write_fasta(tmp_path / "src.fa", "chr", genome) + tgt = _write_multi_fasta(tmp_path / "tgt.fa", [("t1", genome[:1000]), ("t2", genome[1000:])]) + positions = tmp_path / "pos.bed" + positions.write_text("chr\t399\t402\nchr\t1399\t1401\n") + + args = _lift_args(tmp_path, str(positions), src, tgt, global_chain=True, contig=None, + sample=1, indel_tol=0, max_gap=None, min_chain=10, min_density=0.2, + min_identity=0.9, rd_min=50, rd_out=None, source_contig=None) + lift.cmd_lift(args) + + assert (tmp_path / "lift.bed").read_text() == "t1\t399\t402\tblindspot\nt2\t399\t401\tblindspot\n" + + +# ------------------------------------------- --align-gaps: what interpolation cannot place, aligned + + +def test_align_places_a_gap_as_far_left_as_the_scores_allow(): + cols = lift._align("ACGTAAAATGCA", "ACGTAAATGCA") + top = "".join("-" if i is None else "ACGTAAAATGCA"[i] for i, _ in cols) + bottom = "".join("-" if j is None else "ACGTAAATGCA"[j] for _, j in cols) + assert (top, bottom) == ("ACGTAAAATGCA", "ACGT-AAATGCA") + + +def test_align_finds_an_insertion(): + cols = lift._align("ACGTCCGGTTAACCGGTTAAGGCA", "ACGTGGCA") + assert sum(1 for _, j in cols if j is None) == 16 + assert sum(1 for i, _ in cols if i is None) == 0 + + +def test_align_gives_up_past_its_cell_budget(): + assert lift._align("A" * 400, "C" * 5, max_cells=100) is None + + +def _aligned_args(src, tgt, **kw): + return _chain_args(src, tgt, **{"align_gaps": True, "align_min_identity": 0.9, "max_align_gap": 20000, **kw}) + + +def test_positions_beside_an_indel_are_placed_by_aligning_between_the_anchors(tmp_path): + """A 5 bp deletion in the target: interpolation drops the positions around it, the alignment places + them, each at its own coordinate.""" + seq = _random_seq(800, seed=31) + src = _write_fasta(tmp_path / "src.fa", "A", seq) + tgt = _write_fasta(tmp_path / "tgt.fa", "B", seq[:400] + seq[405:]) + positions = [390, 395, 398, 410, 415] + + strict, _ = lift._lift_chain(_chain_args(src, tgt), positions) + good, stats = lift._lift_chain_multi(_aligned_args(src, tgt), [(None, p) for p in positions]) + + assert len(strict) < len(positions), "interpolation alone cannot place them all" + assert good[("A", 390)] == ("B", 390, "+") and good[("A", 395)] == ("B", 395, "+") + assert good[("A", 410)] == ("B", 405, "+") and good[("A", 415)] == ("B", 410, "+") + assert stats["aligned"] >= 1 + + +def test_a_stretch_the_target_lacks_is_reported_absent_and_its_flanks_placed(tmp_path): + """A 300 bp insertion in the source, like an IS copy the target does not have at that site.""" + seq, ins = _random_seq(800, seed=32), _random_seq(300, seed=33) + src = _write_fasta(tmp_path / "src.fa", "A", seq[:400] + ins + seq[400:]) + tgt = _write_fasta(tmp_path / "tgt.fa", "B", seq) + + good, stats = lift._lift_chain_multi(_aligned_args(src, tgt), [(None, p) for p in (395, 550, 705)]) + + assert good[("A", 395)] == ("B", 395, "+") + assert good[("A", 550)] is None, "inside the insertion: the target has no counterpart" + assert good[("A", 705)] == ("B", 405, "+") + assert stats["absent"] == 1 + + +def test_a_stretch_without_anchors_but_the_same_length_is_placed(tmp_path): + """A SNP every 8 bases for 200 bp leaves no shared 21-mer; the alignment is gap-free and places it.""" + seq = list(_random_seq(900, seed=34)) + tgt_seq = seq[:] + for i in range(350, 550, 8): + tgt_seq[i] = {"A": "C", "C": "G", "G": "T", "T": "A"}[tgt_seq[i]] + src = _write_fasta(tmp_path / "src.fa", "A", "".join(seq)) + tgt = _write_fasta(tmp_path / "tgt.fa", "B", "".join(tgt_seq)) + + good, _ = lift._lift_chain_multi(_aligned_args(src, tgt, align_min_identity=0.8), + [(None, 400), (None, 450), (None, 500)]) + + assert good == {("A", p): ("B", p, "+") for p in (400, 450, 500)} + + +def test_an_unrelated_stretch_is_neither_placed_nor_called_absent(tmp_path): + """200 bp of one sequence where the target has 200 bp of another: there is no counterpart to give, + and it is not missing from the target either.""" + seq = _random_seq(900, seed=35) + src = _write_fasta(tmp_path / "src.fa", "A", seq[:350] + _random_seq(200, seed=36) + seq[550:]) + tgt = _write_fasta(tmp_path / "tgt.fa", "B", seq[:350] + _random_seq(200, seed=37) + seq[550:]) + + good, stats = lift._lift_chain_multi(_aligned_args(src, tgt), [(None, 450)]) + + assert good == {} + assert stats["drop"] == 1 + + +def test_the_map_writes_an_absent_position_with_no_target(tmp_path): + out = tmp_path / "map.tsv" + lift._write_map({("A", 5): ("B", 9, "+"), ("A", 7): None}, str(out), ["A"]) + assert out.read_text() == "src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\nA\t5\tB\t9\t+\nA\t7\t.\t.\t.\n" + + # --------------------------------------------------------------------- cmd_markers diff --git a/tests/unit/test_qc_report_canonical.py b/tests/unit/test_qc_report_canonical.py new file mode 100644 index 0000000..ccbd613 --- /dev/null +++ b/tests/unit/test_qc_report_canonical.py @@ -0,0 +1,122 @@ +"""The canonical (H37Rv) coordinate, amino acid and gene the report attaches to each variant (bin/qc_report.py). + +A cohort can be mapped against several references, none of them H37Rv. Each variant's canonical coordinate +comes from the liftover map of ITS reference, looked up by contig and position: the same position on another +contig or reference is another site. A canonical record lends a coordinate only when it was lifted to the +canonical reference, which it says by carrying the mapping coordinate in OPOS. +""" + +from __future__ import annotations + +import argparse + +from conftest import load_script + +qc = load_script("qc_report") + + +def _vcf(path, sample, records): + """records: (chrom, pos, ref, alt, info).""" + head = ("##fileformat=VCFv4.2\n" + f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample}\n") + body = "".join(f"{c}\t{p}\t.\t{r}\t{a}\t50\tPASS\t{i}\tGT:AD:DP\t1:2,38:40\n" for c, p, r, a, i in records) + path.write_text(head + body) + return str(path) + + +def _ann(gene, aa): + return f"ANN=T|missense_variant|MODERATE|{gene}|g|transcript|t|protein_coding|1/1|c.1A>T|{aa}" + + +def _args(vcfs, vcfs_h37rv=(), pos_liftover=None): + return argparse.Namespace(vcfs=list(vcfs), vcfs_h37rv=list(vcfs_h37rv), pos_liftover=pos_liftover, + aa2_label="H37Rv") + + +def _map(path, rows): + path.write_text("src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\n" + + "".join(f"{c}\t{p}\tMTB_anc\t{t}\t{s}\n" for c, p, t, s in rows)) + return str(path) + + +def test_read_lift_map_reads_the_contig_aware_map_and_the_older_one(tmp_path): + new = _map(tmp_path / "new.tsv", [("E1", 100, 761155, "+")]) + old = tmp_path / "old.tsv" + old.write_text("src_pos\ttgt_pos\n100\t761155\n") + + assert qc.read_lift_map(new) == {("E1", 100): 761155} + assert qc.read_lift_map(str(old)) == {(None, 100): 761155} + assert qc.read_lift_map(str(tmp_path / "NO_FILE_LIFTOVER")) == {} + + +def test_each_reference_gets_its_own_canonical_coordinate(tmp_path): + """Two references, the same position: each variant takes the coordinate its own reference lifts to. + One map for the first reference, applied by position alone, gave both the first one's.""" + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", ".")]) + b = _vcf(tmp_path / "b.vcf", "B", [("E2", 100, "C", "T", "."), ("E2", 200, "G", "A", ".")]) + lift = _map(tmp_path / "lift.tsv", [("E1", 100, 761155, "+"), ("E2", 100, 900001, "-")]) + + v = qc.load_variants(_args([a, b], pos_liftover=lift)) + + assert v["A"]["E1:100"]["pos_h37rv"] == "H37Rv:761155" + assert v["B"]["E2:100"]["pos_h37rv"] == "H37Rv:900001" + assert not v["B"]["E2:200"].get("pos_h37rv"), "a position the lift dropped has no canonical coordinate" + + +def test_a_map_without_contigs_is_applied_only_to_a_single_contig_cohort(tmp_path): + old = tmp_path / "old.tsv" + old.write_text("src_pos\ttgt_pos\n100\t761155\n") + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", ".")]) + b = _vcf(tmp_path / "b.vcf", "B", [("E2", 100, "C", "T", ".")]) + + one = qc.load_variants(_args([a], pos_liftover=str(old))) + two = qc.load_variants(_args([a, b], pos_liftover=str(old))) + + assert one["A"]["E1:100"]["pos_h37rv"] == "H37Rv:761155" + assert not two["A"]["E1:100"].get("pos_h37rv") and not two["B"]["E2:100"].get("pos_h37rv") + + +def test_a_lifted_canonical_record_lends_its_coordinate_amino_acid_and_gene(tmp_path): + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", _ann("E1_00667", "p.Ser451Leu"))]) + h = _vcf(tmp_path / "h.vcf", "A", [("Chromosome", 761155, "C", "T", + "OPOS=E1:100;" + _ann("rpoB", "p.Ser450Leu"))]) + + v = qc.load_variants(_args([a], vcfs_h37rv=[h]))["A"]["E1:100"] + + assert (v["aa_h37rv"], v["gene_h37rv"], v["pos_h37rv"]) == ("p.Ser450Leu", "rpoB", "H37Rv:761155") + + +def test_a_canonical_record_annotated_in_place_lends_no_coordinate(tmp_path): + """A VCF annotated against the canonical database without being lifted sits at the mapping coordinate. + Its position is not a canonical one, so it must not stand in for the lift, as it used to.""" + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", ".")]) + h = _vcf(tmp_path / "h.vcf", "A", [("E1", 100, "C", "T", _ann("rpoB", "p.Ser450Leu"))]) + lift = _map(tmp_path / "lift.tsv", [("E1", 100, 761155, "+")]) + + v = qc.load_variants(_args([a], vcfs_h37rv=[h], pos_liftover=lift))["A"]["E1:100"] + + assert v["pos_h37rv"] == "H37Rv:761155", "the lift, not the in-place record's own position" + assert v["aa_h37rv"] == "p.Ser450Leu" + + +def test_an_in_place_record_pairs_by_contig_as_well_as_position(tmp_path): + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", ".")]) + h = _vcf(tmp_path / "h.vcf", "A", [("E2", 100, "C", "T", _ann("rpoB", "p.Ser450Leu"))]) + + v = qc.load_variants(_args([a], vcfs_h37rv=[h]))["A"]["E1:100"] + + assert not v.get("aa_h37rv") + + +def test_a_position_the_canonical_genome_lacks_is_said_to_be_absent(tmp_path): + """The lift writes '.' for a position inserted relative to H37Rv: it has no coordinate, which the report + says instead of leaving it blank like a position the lift could not place.""" + a = _vcf(tmp_path / "a.vcf", "A", [("E1", 100, "C", "T", "."), ("E1", 200, "G", "A", ".")]) + lift = tmp_path / "lift.tsv" + lift.write_text("src_contig\tsrc_pos\ttgt_contig\ttgt_pos\tstrand\nE1\t100\t.\t.\t.\n") + + v = qc.load_variants(_args([a], pos_liftover=str(lift)))["A"] + + assert qc.read_lift_map(str(lift)) == {("E1", 100): None} + assert v["E1:100"]["pos_h37rv"] == "H37Rv:absent" + assert not v["E1:200"].get("pos_h37rv")