Give every variant its H37Rv coordinate and numbering, on any reference

CHANGELOG.md

@@ -84,6 +84,36 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Fixed
 
+- **The H37Rv coordinate and numbering were wrong, or missing, on any reference but H37Rv.** On the
+  185-sample cohort, mapped against two assemblies, the report gave 15,653 of the second reference's
+  16,296 variant sites (97%) a wrong H37Rv coordinate. Four faults:
+
+  - The variant positions of every reference were lifted once, from the first reference's sequence,
+    and the map was applied by position alone. `LIFT_VARIANTS` now runs once per reference, each from
+    its own sequence, and the report looks each variant up by its contig. Where a position cannot be
+    interpolated between two anchors (an indel between them, or a stretch without a shared unique
+    k-mer), the stretch between the anchors is aligned (`--align-gaps`). Every position of 27
+    resistance genes now lifts to the coordinate minimap2 gives it, 99.8% of the report's variant sites
+    do, and none of the rest fits H37Rv worse than minimap2's coordinate; a position H37Rv lacks is
+    shown as *not in H37Rv*, and about 0.2%, in PE_PGRS-type repeats, are left without a coordinate.
+  - Only the first contig of a reference was read, so on a draft assembly a variant of any other
+    contig was lifted with the first contig's sequence: cut into ten contigs, one of the cohort's
+    references had half its positions misplaced. Each contig is now chained on its own, in whichever
+    orientation it lies, and the draft lifts as well as the complete assembly. The blind-spot mask
+    now reaches every contig as well.
+  - `--annotate_canonical`, on in the TB profile, annotated each VCF in place against the H37Rv snpEff
+    database: a position of another reference was read as an H37Rv one, and since the database names
+    H37Rv's chromosome `Chromosome`, which no mapping reference does, it annotated nothing. The
+    unlifted position was then shown as the H37Rv coordinate, over the lift. Each SNP is now moved to
+    its H37Rv position first (`bin/lift_vcf.py`: H37Rv's base as REF, the allele complemented where
+    the reference lies reversed, the mapping coordinate kept in `OPOS`), under the database's name for
+    the chromosome (`--canonical_chrom`, default `Chromosome`), and a record lends only the
+    coordinate it was lifted to.
+  - The SNP-dynamics panel compared a variant's HGVS change (`p.Ile66Met`) with the catalogue's short
+    form (`I66M`), so no trajectory was ever tagged with its catalogue entry. Both are now read in
+    one form, in H37Rv's gene name and numbering first, and a stop or frameshift matches the gene's
+    `LoF` entry.
+
 - **Gene conversion was called from reads that could not carry it, and the cohort multiplied
   it.** On a 185-sample cohort, 1,520 of the 1,579 calls were rows the samples themselves had
   left `ambiguous`, promoted because another sample had called the same stretch outright. Those

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=<ID=ANN,", "##INFO=<ID=LOF,",
+                                    "##INFO=<ID=NMD,", "##INFO=<ID=EFF,", "##INFO=<ID=OPOS,")):
+                head.append(line.rstrip("\n"))
+            continue
+        if line.startswith("#"):
+            chrom_line = line.rstrip("\n")
+            continue
+        c = line.rstrip("\n").split("\t")
+        if len(c) < 8 or not c[1].isdigit():
+            continue
+        ref, alts = c[3].upper(), c[4].upper().split(",")
+        if ref not in COMP or any(a not in COMP for a in alts):
+            n["not_snp"] += 1
+            continue
+        n["snps"] += 1
+        hit = lmap.get((c[0], int(c[1])))
+        if hit is None or hit[0] not in seq or not 1 <= hit[1] <= len(seq[hit[0]]):
+            n["unlifted"] += 1
+            continue
+        tc, tp, strand = hit
+        cref = seq[tc][tp - 1]
+        calts = alts if strand == "+" else [COMP[a] for a in alts]
+        if cref not in COMP or calts[0] == cref:
+            n["canonical_base"] += 1
+            continue
+        info = [f for f in c[7].split(";") if f and f != "." and not f.startswith(_OLD_ANNOTATION)
+                and not f.startswith("OPOS=")]
+        c[7] = ";".join(["OPOS=%s:%s" % (c[0], c[1])] + info)
+        c[0], c[1], c[3], c[4] = rename or tc, str(tp), cref, ",".join(calts)
+        used.add(tc)
+        recs.append(((order[tc], tp), "\t".join(c)))
+        n["lifted"] += 1
+    for tc in sorted(used, key=order.get):
+        head.append("##contig=<ID=%s,length=%d>" % (rename or tc, len(seq[tc])))
+    head.append('##INFO=<ID=OPOS,Number=1,Type=String,Description="The contig:position of the reference the '
+                'sample was mapped against, which this record was lifted from">')
+    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())

bin/pathotypr_liftover.README.md

@@ -7,11 +7,35 @@ flanking-context k-mer and finds it in the target genome; the target position is
 ## `lift` — the recommended method (self-contained, synteny-anchored)
 
 ```bash
-# lift the positions in POS (BED or 1-based list) from reference A onto reference B:
-python3 pathotypr_liftover.py lift POS A.fasta B.fasta --out-map map.tsv --out-bed lifted.bed --contig B \
+# lift the positions in POS (BED, 'contig<TAB>pos' or a 1-based list) from reference A onto reference B:
+python3 pathotypr_liftover.py lift POS A.fasta B.fasta --out-map map.tsv --out-bed lifted.bed \
         --kmer-size 21 --global-chain
 ```
 
+With `--global-chain` every contig of both FASTAs is read. Each source contig is chained on its own, so the
+contigs of a draft assembly are placed wherever, and in whichever orientation, they lie in the target; the map
+names both contigs and the strand (`src_contig  src_pos  tgt_contig  tgt_pos  strand`) and the BED is written
+per target contig (`--contig` renames a single-contig target). A position names its contig; on a single-contig
+source it may name any (a BED in NC_000962.3 coordinates read on an ancestor FASTA named otherwise), on a
+multi-contig one a position without a contig, or on a contig the source lacks, is left out and counted. The
+other modes read only the first contig of each FASTA.
+
+**`--align-gaps` (what `LIFT_VARIANTS` uses).** Interpolation only crosses an exactly colinear gap of ≤ ~2k, so a
+position beside an indel, or in a stretch with no shared unique k-mer (SNP-dense, repeated), drops. With
+`--align-gaps` the stretch between the two anchors is aligned instead (global, affine gaps with BWA-MEM's scores,
+banded, indels left-aligned as VCF normalisation does; pure Python, `--max-align-gap` 20 kb):
+
+- a position is placed when a gap-free run of the alignment at ≥ `--align-min-identity` (0.9) joins it to one of
+  the two anchors, and its own 10 bases on each side match — or, on the side where they do not, a short indel
+  (≤ 20 bp) sits right there, or a long one beside which its context occurs only once in the stretch;
+- a position inside a run the target lacks is written as absent (`.` as its target) when that run is the only
+  indel between the anchors, its flanks align at ≥ 90% and no anchor inside the gap points elsewhere;
+- anything else — between two indels, beside a long indel in a repeat — is left out.
+
+Against minimap2 on two MTBC assemblies (L7 and A4) lifted to H37Rv, every position of 27 resistance genes lands
+on minimap2's coordinate, 99.8% of random positions do, and of those that differ none fits H37Rv worse than
+minimap2's; about 0.2% are left out, in PE_PGRS-type repeats. It adds ~10 s to a 4.4 Mb lift.
+
 **`--global-chain` (recommended, what BAMpiro uses)** builds a whole-genome coordinate map: every k-mer
 unique in BOTH genomes is an anchor (~97 % of the MTBC genome — anchors roughly every base), chained into a
 collinear order (LIS), and any position is placed by **interpolating between its flanking anchors**. So a
@@ -97,14 +121,16 @@ Flag style is mixed in pathotypr 0.1.0: `--tsv_pos` / `--ref_fasta` / `--tsv_gen
 but `--kmer-size` / `--fasta-genomes` use dashes. The classify main output is written to the `--output`
 name exactly (no extension); its columns are `genome  k-mer  k-merPOS  SNPgenome  SNPreference  lineage`.
 
-`map.tsv` = `src_pos <TAB> tgt_pos`; `lifted.bed` = the positions collapsed into intervals in B's coords.
+`map.tsv` = `src_pos <TAB> 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.