From 7731e7c4f3dbe06d205d6f2136cd6e3f6309ee95 Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Wed, 23 Sep 2026 17:44:44 +0200 Subject: [PATCH 1/4] Carry every reference's variants to H37Rv from its own sequence The H37Rv coordinate of each variant was lifted once, from the first reference's sequence, for every reference of the run, and the map was applied by position alone. On the cohort mapped against two assemblies, 15,653 of the second reference's 16,296 variant sites (97%) got a wrong coordinate. LIFT_VARIANTS now runs once per reference, from its own sequence and for the positions on its own contigs, and the report looks each variant up by contig and position: 97% of the second reference's sites lift to the coordinate minimap2 gives them, and the rest drop. The liftover read only the first contig of each FASTA, 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. Every contig is now read and chained on its own, in whichever orientation it lies in the target, and the map names both contigs and the strand. The draft lifts as well as the complete assembly, a single-contig reference gets the same map as before, and the blind-spot mask reaches every contig. The canonical annotation ran the mapping-reference VCF through the H37Rv snpEff database in place: another reference's position was read as an H37Rv one, and since the database calls H37Rv's chromosome 'Chromosome' and a mapping reference is called something else, nothing was annotated. The unlifted position was then shown as the H37Rv coordinate, over the lift. lift_vcf.py now moves each SNP to its H37Rv position first, with H37Rv's base as REF and the allele complemented where the reference lies reversed, under the database's name for the chromosome (--canonical_chrom), and keeps the mapping coordinate in OPOS. A record lends only the coordinate it was lifted to, and the report keeps H37Rv's gene name beside its amino acid. --- bin/lift_vcf.py | 138 ++++++++++++ bin/pathotypr_liftover.README.md | 22 +- bin/pathotypr_liftover.py | 279 ++++++++++++++++++++----- bin/qc_report.py | 100 ++++++--- bin/qcreport/panels.py | 8 +- conf/organism.config | 6 +- conf/tuberculosis.config | 9 +- docs/configuration.md | 9 +- docs/outputs.md | 2 +- docs/pathotypr.md | 46 ++-- docs/qc-report.md | 4 +- docs/tutorials/other-organisms.md | 2 +- modules/annotation.nf | 28 ++- modules/reference.nf | 5 +- modules/report.nf | 38 ++-- nextflow.config | 24 ++- subworkflows/cohort_report.nf | 24 ++- tests/unit/test_lift_vcf.py | 98 +++++++++ tests/unit/test_pathotypr_liftover.py | 88 +++++++- tests/unit/test_qc_report_canonical.py | 108 ++++++++++ 20 files changed, 864 insertions(+), 174 deletions(-) create mode 100755 bin/lift_vcf.py create mode 100644 tests/unit/test_lift_vcf.py create mode 100644 tests/unit/test_qc_report_canonical.py 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..39d59d9 100644 --- a/bin/pathotypr_liftover.README.md +++ b/bin/pathotypr_liftover.README.md @@ -7,11 +7,19 @@ 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. + **`--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 +105,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..26199a0 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,26 @@ 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 +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 +425,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 +441,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 +483,75 @@ 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: + good, n_fwd, n_rev, n_drop = {}, 0, 0, n_outside + 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: + 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 + loc = T.local(best) if best is not None else None + if loc is None: n_drop += 1 - elif best_orient == 1: - good[p] = best; n_fwd += 1 + continue + good[(S.names[ci], p)] = (T.names[loc[0]], loc[1], "+" if best_orient == 1 else "-") + 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), + "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: t for (_, p), (_, t, _) in good.items()}, st + + +def _write_map(good, path, source_contigs): + """The contig-aware map: src_contig src_pos tgt_contig tgt_pos strand, in source order.""" + 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])): + 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 tc, tp, _ in good.values(): + by[tc].add(tp) + 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 +561,20 @@ 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["anchors"], st["inv_chains"], st["inv_anchors"], st["rd"])) + 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 +631,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 +669,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, @@ -541,7 +707,8 @@ def main(): 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..cad58f1 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,77 @@ def _index_dir(d): return dmg_by_sample +def read_lift_map(path): + """{(contig, position): canonical position} 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. 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) >= 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 + _t = _lift.get((_c, int(_p))) + if _t is None and legacy: + _t = _lift.get((None, int(_p))) + if _t is not None: + _v["pos_h37rv"] = "%s:%s" % (label, _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/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..d84164d 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,19 @@ 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. + +On the two references of a 185-sample cohort, checked against minimap2's alignment of +each to H37Rv, 97% of the variant sites lift to the same coordinate and the rest are +dropped; the handful that disagree sit in repeats, at a coordinate whose sequence +matches exactly. Cut into ten contigs with six of them reversed, the same reference +lifts just as well. ## Blind-spot masking (H37Rv problematic sites) @@ -96,7 +114,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..1232222 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -67,41 +67,43 @@ 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. Pure Python (bin/pathotypr_liftover.py lift + // --global-chain), about 15 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 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/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..155e2aa 100644 --- a/tests/unit/test_pathotypr_liftover.py +++ b/tests/unit/test_pathotypr_liftover.py @@ -629,10 +629,96 @@ 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" + + # --------------------------------------------------------------------- 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..f97b1e1 --- /dev/null +++ b/tests/unit/test_qc_report_canonical.py @@ -0,0 +1,108 @@ +"""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") From 06fe4a02c31f75149fa2de0639c6db7538784d7d Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Wed, 23 Sep 2026 17:44:57 +0200 Subject: [PATCH 2/4] Tag a trajectory with its catalogue entry in either notation The SNP-dynamics panel compared a variant's HGVS change (p.Ile66Met) with the resistance catalogue's short form (I66M), so no trajectory was ever tagged with its catalogue entry, on any reference. Both are now read in one form, a stop also as the catalogue's '!', and a stop or frameshift matches the gene's LoF entry. The H37Rv gene and numbering are tried first, since the catalogue is written in them and a reference of another lineage names and numbers its genes its own way. --- bin/report_assets/js/10_dynamics.js | 18 +++++++- tests/js/dynamics.test.mjs | 64 +++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+), 2 deletions(-) create mode 100644 tests/js/dynamics.test.mjs 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/tests/js/dynamics.test.mjs b/tests/js/dynamics.test.mjs new file mode 100644 index 0000000..4ca5b15 --- /dev/null +++ b/tests/js/dynamics.test.mjs @@ -0,0 +1,64 @@ +// Matching a trajectory's variant to the drug-resistance catalogue (10_dynamics.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); + }); +}); From 600fdcd82af5eed17f6650f20f4e8ed997ae09e2 Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Wed, 23 Sep 2026 17:45:09 +0200 Subject: [PATCH 3/4] Record the H37Rv coordinate fixes in the changelog --- CHANGELOG.md | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e966a2..4f62882 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -84,6 +84,32 @@ 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: 97% of the second + reference's sites lift to the coordinate minimap2 gives them, the rest are dropped. + - 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 From 0dd005b06f338626604f1c7a27172610f3d6b1fe Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Wed, 23 Sep 2026 18:32:17 +0200 Subject: [PATCH 4/4] Place what the lift used to drop by aligning between its anchors The lift interpolated only across an exactly colinear gap of at most ~2k between two shared unique k-mers, so a position beside an indel, or in a stretch without a shared unique k-mer, was dropped: 4 to 5% of the variant sites of the cohort's two references. With --align-gaps, which LIFT_VARIANTS now passes, the stretch between the two anchors is aligned instead, and a position is placed when a gap-free run of the alignment joins it to one of the anchors and its own context agrees. A position inside a stretch H37Rv lacks is written as such, and the report shows it as not in H37Rv. Checked against minimap2's alignment of each reference to H37Rv, every position of 27 resistance genes lifts to minimap2's coordinate, and so do 99.8% of the variant sites. None of those that differ fits H37Rv worse: they are ties in identical repeat units, bases at the junction of a structural difference, or places where minimap2 aligned another copy. About 0.2% are left without a coordinate, in PE_PGRS-type repeats where a position could sit on either side of an indel. --- CHANGELOG.md | 8 +- bin/pathotypr_liftover.README.md | 16 ++ bin/pathotypr_liftover.py | 279 ++++++++++++++++++++++++- bin/qc_report.py | 24 ++- bin/report_assets/js/01_prelude.js | 7 +- docs/pathotypr.md | 20 +- modules/report.nf | 8 +- tests/js/dynamics.test.mjs | 16 +- tests/unit/test_pathotypr_liftover.py | 89 ++++++++ tests/unit/test_qc_report_canonical.py | 14 ++ 10 files changed, 451 insertions(+), 30 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4f62882..bedc6f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -90,8 +90,12 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow - 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: 97% of the second - reference's sites lift to the coordinate minimap2 gives them, the rest are dropped. + 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 diff --git a/bin/pathotypr_liftover.README.md b/bin/pathotypr_liftover.README.md index 39d59d9..057178d 100644 --- a/bin/pathotypr_liftover.README.md +++ b/bin/pathotypr_liftover.README.md @@ -20,6 +20,22 @@ source it may name any (a BED in NC_000962.3 coordinates read on an ancestor FAS 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 diff --git a/bin/pathotypr_liftover.py b/bin/pathotypr_liftover.py index 26199a0..75efa53 100644 --- a/bin/pathotypr_liftover.py +++ b/bin/pathotypr_liftover.py @@ -395,6 +395,209 @@ def _chains(pairs, increasing, min_anchors, max_chains=1024): return out +_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 @@ -487,7 +690,29 @@ def bracket(ca, cb, p, orient): 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, n_outside + 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 @@ -498,11 +723,29 @@ def bracket(ca, cb, p, orient): 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: - n_drop += 1 + 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: @@ -513,22 +756,26 @@ def bracket(ca, cb, p, orient): 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), - "source_contigs": S.names, "target_contigs": T.names} + "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: t for (_, p), (_, t, _) in good.items()}, st + 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.""" + """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)) @@ -537,8 +784,9 @@ 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 tc, tp, _ in good.values(): - by[tc].add(tp) + 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))): @@ -569,8 +817,11 @@ def cmd_lift(args): 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, @@ -701,6 +952,16 @@ 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, diff --git a/bin/qc_report.py b/bin/qc_report.py index cad58f1..8d526c7 100644 --- a/bin/qc_report.py +++ b/bin/qc_report.py @@ -166,12 +166,14 @@ def _index_dir(d): def read_lift_map(path): - """{(contig, position): canonical position} from the liftover map, or {} without one. + """{(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. 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.""" + 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 @@ -181,6 +183,8 @@ def read_lift_map(path): 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: @@ -231,11 +235,13 @@ def load_variants(args): _c, _, _p = _key.rpartition(":") if not _p.isdigit(): continue - _t = _lift.get((_c, int(_p))) - if _t is None and legacy: - _t = _lift.get((None, int(_p))) - if _t is not None: - _v["pos_h37rv"] = "%s:%s" % (label, _t) + _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/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/docs/pathotypr.md b/docs/pathotypr.md index d84164d..899aae9 100644 --- a/docs/pathotypr.md +++ b/docs/pathotypr.md @@ -100,11 +100,23 @@ per reference), and every contig of it: each contig of a draft assembly is chain 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, 97% of the variant sites lift to the same coordinate and the rest are -dropped; the handful that disagree sit in repeats, at a coordinate whose sequence -matches exactly. Cut into ten contigs with six of them reversed, the same reference -lifts just as well. +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) diff --git a/modules/report.nf b/modules/report.nf index 1232222..1bcd5be 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -73,8 +73,10 @@ process LIFT_VARIANTS { // 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. Pure Python (bin/pathotypr_liftover.py lift - // --global-chain), about 15 s and 2.5 GB for a 4.4 Mb genome. + // 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 } @@ -95,7 +97,7 @@ process LIFT_VARIANTS { | sort -u > positions.txt if [ -s positions.txt ]; then python3 ${projectDir}/bin/pathotypr_liftover.py lift positions.txt ${ref_fa} ${params.canonical_ref} \\ - --out-map ${refId}.pos_liftover.tsv --kmer-size 21 --global-chain + --out-map ${refId}.pos_liftover.tsv --kmer-size 21 --global-chain --align-gaps else printf 'src_contig\\tsrc_pos\\ttgt_contig\\ttgt_pos\\tstrand\\n' > ${refId}.pos_liftover.tsv fi diff --git a/tests/js/dynamics.test.mjs b/tests/js/dynamics.test.mjs index 4ca5b15..51a2eb5 100644 --- a/tests/js/dynamics.test.mjs +++ b/tests/js/dynamics.test.mjs @@ -1,4 +1,5 @@ -// Matching a trajectory's variant to the drug-resistance catalogue (10_dynamics.js). +// 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 @@ -62,3 +63,16 @@ describe("dynDRmatch", () => { 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_pathotypr_liftover.py b/tests/unit/test_pathotypr_liftover.py index 155e2aa..bf0d727 100644 --- a/tests/unit/test_pathotypr_liftover.py +++ b/tests/unit/test_pathotypr_liftover.py @@ -719,6 +719,95 @@ def test_cmd_lift_global_chain_writes_a_bed_per_target_contig(tmp_path, capsys): 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 index f97b1e1..ccbd613 100644 --- a/tests/unit/test_qc_report_canonical.py +++ b/tests/unit/test_qc_report_canonical.py @@ -106,3 +106,17 @@ def test_an_in_place_record_pairs_by_contig_as_well_as_position(tmp_path): 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")