Deliver indels, and read two SNPs of one codon as one change

.github/dockerfile/Dockerfile

@@ -49,6 +49,7 @@ ENV MAMBA_ROOT_PREFIX=/opt/conda
 # - kraken2: Decontamination
 # - multiqc: Reporting
 # - pathotypr: alignment-free MTBC lineage + WHO drug-resistance genotyping (k-mer, reference-agnostic)
+# - get_mnv: codon-level (MNV) amino-acid changes, phased on the reads (PathoGenOmics-Lab)
 # - python=3.12: Latest stable python requested
 RUN micromamba create -n bampiro-env -c conda-forge -c bioconda -y \
     'python' \
@@ -73,6 +74,9 @@ RUN micromamba create -n bampiro-env -c conda-forge -c bioconda -y \
     # Pinned, like its marker files below. The catalogue and the binary carry the same
     # version on purpose: v1.0.2 renumbered the data files so the two match again.
     'pathotypr=1.0.2' \
+    # Codon-level changes: two SNPs in one codon read whole, on the reads that carry them (MNVs).
+    # Pinned, like pathotypr: the codon table and its columns are what the run's outputs are read from.
+    'get_mnv=1.1.5' \
     'pandas' \
     # pandas already pulls numpy in, but bin/gconv_model.py and bin/build_min_unique_len.py import
     # it directly, and a direct import has no business resting on somebody else's dependency.

CHANGELOG.md

@@ -8,6 +8,33 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Added
 
+- **Indels, as a deliverable.** Indels are called with the definition a SNP has (the depth, strand
+  support and allele fraction of the hom and het rules), and each sample's
+  `<sample>.<ref>.indels.ann.vcf.gz` keeps every indel FreeBayes made, normalised, with `FILTER`
+  `PASS` where it meets the rule and `LowSupport` where not: a minority frameshift, the way
+  resistance through Rv0678 often starts, keeps its fraction instead of vanishing. Until now the
+  main VCF held SNPs only and a minority indel was in no filtered file. `<samplesheet>_indel_matrix.tsv`
+  holds the run's indels as the SNP matrix holds its SNPs: every indel a sample calls with a `PASS`,
+  its length, gene, effect and HGVS, and per sample the fraction, the depth and whether that call
+  passes; a sample without the indel reads AF 0 at the depth of the anchor base where it was read.
+  The report's matrix panel gains an Indels view, with frameshifts marked. The legacy split set
+  gains `var.het.indel.vcf`. The matrix is off with `--make_indel_matrix false`; the per-sample
+  indel VCFs are always written.
+
+- **Codon-level changes, read on the reads that carry them.** Two SNPs in one codon change one
+  amino acid, and annotated a base at a time they name two changes that are not there. A new step,
+  `GET_MNV`, runs [get_MNV](https://github.com/PathoGenOmics-Lab/get_MNV) 1.1.5 on each sample's
+  passing SNPs and indels with its alignment, and reads every codon whole on the reads that span it:
+  per-sample tables under `mnv/`, and `<samplesheet>_mnv.tsv` with the codons at least
+  `--mnv_min_reads` reads carry whole. The SNP matrix gains `codon_change` and
+  `codon_change_samples`, and the report marks those SNPs and their cells. On reads simulated from
+  H37Rv, two SNPs of rpoB codon 445 on the same reads read His445Glu (base by base: His445Asp and
+  His445Gln), a minority pair in codon 460 reads Glu460Leu where its first base alone says a stop,
+  and two SNPs of codon 430 on different molecules stay two changes: get_MNV writes a combined
+  change for them too, with no read carrying it, and the table leaves it out. get_MNV runs in its
+  bioconda image, pinned by digest (`--mnv_container`), until the pipeline image carries it; the
+  Dockerfile now installs it. Off with `--run_mnv false`.
+
 - **How low an allele frequency can be trusted.** A *Minority variants* panel counts each
   sample's calls below fixation, their spread of allele fractions and how many rest on three
   alternate reads or fewer, where an error and a real minority look the same. When the
@@ -84,6 +111,16 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Fixed
 
+- **The report left out every SNP FreeBayes wrote as an MNP.** Two changes of one codon on the same
+  reads come out of FreeBayes as one record (CAC>GAG), and the report read only single-base records,
+  so those SNPs were in the master SNP matrix and in no table of the report: not in the variant
+  matrix, the series, the minority variants or the dynamics. It now reads such a record as the
+  SNPs it is made of, each with the record's fraction and its annotation, the codon read whole, as
+  the SNP matrix does.
+
+- **The outputs page said the main VCF held indels.** `<sample>.<ref>.ann.vcf.gz` holds SNPs, one
+  record per base; the indels are in `<sample>.<ref>.indels.ann.vcf.gz` (see Added).
+
 - **The Garnatxa launcher asked for a branch that no longer exists.** `conf/garnatxa.sbatch` ran
   `feat/gene-conversion`, deleted once it was merged, which Nextflow can no longer pull. It now
   runs `main`, pulled again on every launch (`-latest`): without that, Nextflow runs the copy it

README.md

@@ -45,6 +45,7 @@ tuberculosis* but is **organism-agnostic** - point it at any reference genome +
 - Alignment-free MTBC lineage + WHO drug-resistance typing ([Pathotypr](docs/pathotypr.md))
 - A self-contained, interactive [HTML QC report](docs/qc-report.md) with 21 panels
 - Dual amino-acid numbering (used reference + H37Rv / Mycobrowser)
+- Indels and codon-level (MNV) amino-acid changes, phased on the reads
 - One digest-pinned container with every tool and marker panel built in
 
 ## Features
@@ -54,6 +55,7 @@ tuberculosis* but is **organism-agnostic** - point it at any reference genome +
 | 🧬 Any bacterial genome | Reference-agnostic mapping + variant calling; TB-tuned defaults, works on any species |
 | 🧹 Repeat & mappability masking | `nucmer` repeat exclusion plus a length-aware `genmap` read filter |
 | 🧪 Variants & backbone | `FreeBayes` (ploidy 1/2) + "all-sites" VCFs for phylogenetic supermatrices |
+| 🧩 Indels & codon-level changes | Every indel with the call rules a SNP has (PASS or LowSupport) + an indel matrix; SNPs of one codon read whole on the reads that carry them ([get_MNV](https://github.com/PathoGenOmics-Lab/get_MNV)) |
 | 🩺 Lineage & drug resistance | Alignment-free MTBC lineage + WHO DR typing ([Pathotypr](docs/pathotypr.md)), reference-agnostic |
 | 📊 Interactive QC report | Self-contained HTML dashboard, [21 panels](docs/qc-report.md) + per-sample `qc_flags.tsv` |
 | 🔤 Dual amino-acid numbering | Protein changes in both the used reference and H37Rv/Mycobrowser numbering |

assets/NO_FILE.README.md

@@ -13,6 +13,8 @@ stages one of these empty files instead, and the process recognises it by the
 | `NO_FILE_LIFTOVER` | variant coordinate map (`QC_REPORT`) | `--variant_liftover false` |
 | `NO_FILE_DISTANCES` | pairwise SNP distances (`QC_REPORT`) | `--make_snp_distances false`, or no consensus sequences |
 | `NO_FILE_MATRIX` | the master SNP matrix (`QC_REPORT`) | `--make_snp_matrix false` |
+| `NO_FILE_INDELS` | the master indel matrix (`QC_REPORT`) | `--make_indel_matrix false` |
+| `NO_FILE_MNV` | the codon-level table (`SNP_MATRIX`, `QC_REPORT`) | `--run_mnv false`, or no sample reached get_MNV |
 
 Two properties matter, and both are load-bearing:
 

bin/build_indel_matrix.py

@@ -0,0 +1,167 @@
+#!/usr/bin/env python3
+"""Master indel matrix.
+
+The indels of a run in one wide TSV, as the SNP matrix holds its SNPs:
+  rows    = every indel (contig, position, REF, ALT) that at least one sample calls with a PASS:
+            the hom or het rule of bin/vcf_filter_rules.py, the depth, strand support and allele
+            fraction a SNP needs
+  columns = the reference and the annotation (the indel's length, gene, effect, HGVS.c, HGVS.p),
+            how many samples call it with a PASS, then THREE columns per sample: allele fraction
+            (|AF), depth (|DP) and whether that sample's own call passes (|FT: PASS or LowSupport)
+
+It reads the per-sample indel VCFs of CALL_FREEBAYES (every normalised indel call, soft-filtered),
+annotated by snpEff when the run annotates. A sample that calls a row's indel below the rule keeps
+its fraction, marked LowSupport: a minority indel is where resistance through a frameshift often
+starts, and a blank would say the sample does not have it.
+
+A sample with no call at all reads as in the SNP matrix, from its all-positions VCF (--depth-vcfs):
+AF 0 with the depth at the indel's anchor base where it was read, AF blank and DP 0 where no read
+covers it, both blank where the file has no record there (a site of another reference). That depth
+counts the reads that show a base at the anchor, so for a long indel it is a floor: a read that
+ends inside the event is not counted.
+"""
+from __future__ import annotations
+
+import argparse
+import os
+import sys
+from concurrent.futures import ProcessPoolExecutor
+
+from build_snp_matrix import _af_dp, _init_worker, _open, _read_depths, _site_index
+
+PASSING = ("PASS", ".", "")        # an unfiltered VCF (no soft filter written) reads as passing
+
+
+def _ann(info):
+    """(gene, effect, HGVS.c, HGVS.p) from the first snpEff ANN entry."""
+    for f in info.split(';'):
+        if f.startswith('ANN='):
+            a = f[4:].split(',')[0].split('|')
+            get = lambda k: a[k] if len(a) > k else ''
+            return get(3), get(1), get(9), get(10)
+    return '', '', '', ''
+
+
+def read_indels(paths):
+    """(samples in first-seen order, {(contig, pos, ref, alt): site}) from per-sample VCFs."""
+    samples, sites = [], {}
+    for p in paths:
+        if not p or not os.path.exists(p):
+            continue
+        try:
+            with _open(p) as fh:
+                sample = None
+                for line in fh:
+                    if line.startswith('##'):
+                        continue
+                    if line.startswith('#CHROM'):
+                        cols = line.rstrip('\n').split('\t')
+                        sample = cols[9] if len(cols) > 9 else os.path.basename(p).split('.')[0]
+                        if sample not in samples:
+                            samples.append(sample)
+                        continue
+                    if sample is None:
+                        continue
+                    c = line.rstrip('\n').split('\t')
+                    if len(c) < 8:
+                        continue
+                    contig, ref, alt = c[0], c[3], c[4].split(',')[0]
+                    if alt in ('.', '', '*') or alt.startswith('<') or len(ref) == len(alt):
+                        continue        # SNPs and MNPs are the SNP matrix's
+                    try:
+                        pos = int(c[1])
+                    except ValueError:
+                        continue
+                    af, dp = _af_dp(c[8], c[9]) if len(c) >= 10 else (1.0, None)
+                    if af is None:
+                        continue
+                    st = sites.setdefault((contig, pos, ref, alt), {'gene': '', 'eff': '', 'hgvs_c': '',
+                                                                    'hgvs_p': '', 'cells': {}})
+                    gene, eff, hc, hp = _ann(c[7])
+                    for k, v in (('gene', gene), ('eff', eff), ('hgvs_c', hc), ('hgvs_p', hp)):
+                        if v and not st[k]:
+                            st[k] = v
+                    if sample not in st['cells']:
+                        st['cells'][sample] = (af, dp, c[6] if len(c) > 6 else '')
+        except OSError as e:
+            sys.stderr.write(f"[indel_matrix] WARN could not read {p}: {e}\n")
+    return samples, sites
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description="Build a master indel matrix from per-sample indel VCFs.")
+    ap.add_argument('--vcfs', nargs='+', required=True, help="per-sample indel VCFs (annotated or not)")
+    ap.add_argument('--reference', default='', help="reference id the samples were mapped against")
+    ap.add_argument('--depth-vcfs', nargs='*', default=[],
+                    help="per-sample all-positions VCFs: the depth at the anchor base of an indel a "
+                         "sample does not call")
+    ap.add_argument('--threads', type=int, default=1, help="all-positions VCFs read in parallel")
+    ap.add_argument('-o', '--output', required=True)
+    args = ap.parse_args(argv)
+
+    samples, sites = read_indels(args.vcfs)
+    keys = sorted(k for k, st in sites.items()
+                  if any(ft in PASSING for _, _, ft in st['cells'].values()))
+
+    # the depth at each row's anchor base, for the samples that do not call it
+    anchors = sorted({(contig, pos) for contig, pos, _, _ in keys})
+    at = {k: i for i, k in enumerate(anchors)}
+    depths = {}
+    depth_vcfs = [p for p in args.depth_vcfs if p and os.path.exists(p)]
+    if depth_vcfs and anchors:
+        index = _site_index(anchors)
+        if args.threads > 1 and len(depth_vcfs) > 1:
+            with ProcessPoolExecutor(max_workers=args.threads, initializer=_init_worker,
+                                     initargs=(index, len(anchors))) as pool:
+                results = list(pool.map(_read_depths, depth_vcfs))
+        else:
+            results = [_read_depths(p, index, len(anchors)) for p in depth_vcfs]
+        for res in results:
+            if res is None:
+                continue
+            sample, d, _, _ = res
+            if sample not in samples:
+                continue
+            have = depths.get(sample)
+            if have is None:
+                depths[sample] = d
+            else:          # one sample mapped against two references: disjoint contigs
+                for i, v in enumerate(d):
+                    if v >= 0 and have[i] < 0:
+                        have[i] = v
+
+    header = ['reference', 'contig', 'pos', 'ref_allele', 'alt_allele', 'length', 'gene', 'effect',
+              'hgvs_c', 'hgvs_p', 'n_pass']
+    for s in samples:
+        header += [f'{s}|AF', f'{s}|DP', f'{s}|FT']
+    with open(args.output, 'w') as out:
+        out.write('\t'.join(header) + '\n')
+        for contig, pos, ref, alt in keys:
+            st = sites[(contig, pos, ref, alt)]
+            n_pass = sum(1 for _, _, ft in st['cells'].values() if ft in PASSING)
+            row = [args.reference or contig, contig, str(pos), ref, alt, str(len(alt) - len(ref)),
+                   st['gene'], st['eff'], st['hgvs_c'], st['hgvs_p'], str(n_pass)]
+            i = at[(contig, pos)]
+            for s in samples:
+                if s in st['cells']:
+                    af, dp, ft = st['cells'][s]
+                    row += [f'{af:.4f}', '' if dp is None else str(dp), 'PASS' if ft in PASSING else ft]
+                    continue
+                d = depths.get(s)
+                dp = d[i] if d is not None else -1
+                if dp > 0:
+                    row += ['0', str(dp), '']      # read at the anchor, and no such indel called
+                elif dp == 0:
+                    row += ['', '0', '']           # nobody read it: absence says nothing
+                else:
+                    row += ['', '', '']            # no record: another reference, or no depths given
+            out.write('\t'.join(row) + '\n')
+    n_low = sum(1 for k in keys for _, _, ft in sites[k]['cells'].values() if ft not in PASSING)
+    sys.stderr.write(f"[indel_matrix] {len(keys)} indel(s) called with a PASS in at least one of "
+                     f"{len(samples)} sample(s) ({len(sites) - len(keys)} more only below the rule, left "
+                     f"out; {n_low} cell(s) kept below it) -> {args.output}\n")
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())

bin/build_snp_matrix.py

@@ -27,6 +27,14 @@
 
 Only SNPs (single-base REF and ALT) are kept. The '--reference' value labels the
 reference the samples were mapped against.
+
+A SNP that shares its codon with another one, on the same reads, is half of one amino-acid change:
+annotated alone it names a change that is not there (a stop, where the codon read whole gives a
+leucine). With the run's codon-level table (--mnv, collect_mnv.py over get_MNV), two columns after
+aa_change say so: codon_change, the change of every such codon the site takes part in, and
+codon_change_samples, the samples whose reads carry it with the fraction that does, as
+`His445Glu=S1:0.99,S2:1.00` (codons separated by ';'). Both are empty for a SNP that stands alone,
+and neither column is written without the table, so a run without get_MNV keeps the old layout.
 """
 from __future__ import annotations
 import argparse
@@ -200,6 +208,33 @@ def _ann(info):
     return '', '', ''
 
 
+def read_mnv(path):
+    """{(contig, pos): {aa_change: [(sample, mnv_frequency)]}} from a codon-level table (collect_mnv.py)."""
+    out = {}
+    if not path or not os.path.exists(path) or os.path.basename(path).startswith('NO_FILE'):
+        return out
+    import csv
+    with open(path, newline='', encoding='utf-8', errors='replace') as fh:
+        for r in csv.DictReader(fh, delimiter='\t'):
+            aa = (r.get('aa_change') or '').strip()
+            if not aa:
+                continue
+            for p in (r.get('positions') or '').split(','):
+                p = p.strip()
+                if p.isdigit():
+                    out.setdefault((r.get('contig', ''), int(p)), {}).setdefault(aa, []).append(
+                        (r.get('sample', ''), (r.get('mnv_frequency') or '').strip()))
+    return out
+
+
+def _codon_cells(mnv):
+    """(codon_change, codon_change_samples) of one site."""
+    if not mnv:
+        return '', ''
+    aas = sorted(mnv)
+    return ';'.join(aas), ';'.join(f"{aa}=" + ','.join(f"{s}:{f}" if f else s for s, f in mnv[aa]) for aa in aas)
+
+
 def main():
     ap = argparse.ArgumentParser(description="Build a master SNP matrix from per-sample annotated VCFs.")
     ap.add_argument('--vcfs', nargs='+', required=True, help="per-sample (annotated) VCFs")
@@ -210,8 +245,13 @@ def main():
                          "blank, and a blank is left only where nothing was read")
     ap.add_argument('--threads', type=int, default=1,
                     help="all-positions VCFs read in parallel")
+    ap.add_argument('--mnv', default=None,
+                    help="the run's codon-level table (collect_mnv.py): which SNPs are one amino-acid change "
+                         "together with another one on the same reads")
     ap.add_argument('-o', '--output', required=True)
     args = ap.parse_args()
+    with_codons = bool(args.mnv) and os.path.exists(args.mnv) and not os.path.basename(args.mnv).startswith('NO_FILE')
+    mnv = read_mnv(args.mnv) if with_codons else {}
 
     samples = []           # preserve first-seen order
     sites = {}             # (contig, pos) -> {ref, alt:set, gene, eff, aa, cells:{sample:(af,dp)}}
@@ -314,6 +354,8 @@ def main():
                 f"keyed on contig and position, cannot tell their sites apart\n")
 
     header = ['reference', 'contig', 'pos', 'ref_allele', 'alt_allele', 'gene', 'effect', 'aa_change']
+    if with_codons:
+        header += ['codon_change', 'codon_change_samples']
     for s in samples:
         header += [f'{s}|AF', f'{s}|DP']
     n_absent = n_unread = 0
@@ -323,6 +365,8 @@ def main():
             st = sites[(contig, pos)]
             row = [args.reference or contig, contig, str(pos), st['ref'],
                    ','.join(sorted(st['alt'])), st['gene'], st['eff'], st['aa']]
+            if with_codons:
+                row += _codon_cells(mnv.get((contig, pos)))
             for s in samples:
                 if s in st['cells']:
                     af, dp = st['cells'][s]
@@ -340,7 +384,9 @@ def main():
                 else:
                     row += ['', '']          # no record: another reference, or no depths given
             out.write('\t'.join(row) + '\n')
-    sys.stderr.write(f"[snp_matrix] {len(sites)} SNP site(s) x {len(samples)} sample(s) -> {args.output}\n")
+    in_codon = sum(1 for k in keys if k in mnv)
+    sys.stderr.write(f"[snp_matrix] {len(sites)} SNP site(s) x {len(samples)} sample(s)"
+                     + (f", {in_codon} of them part of a codon-level change" if in_codon else "") + f" -> {args.output}\n")
     if args.depth_vcfs and keys:
         missing = [s for s in samples if s not in depths]
         more = f" and {len(missing) - 5} more" if len(missing) > 5 else ""