diff --git a/.github/dockerfile/Dockerfile b/.github/dockerfile/Dockerfile
index 340703f..ca60f0e 100644
--- a/.github/dockerfile/Dockerfile
+++ b/.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.
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f0b50ba..9afafaf 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/README.md b/README.md
index d608736..dcfe657 100644
--- a/README.md
+++ b/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 |
diff --git a/assets/NO_FILE.README.md b/assets/NO_FILE.README.md
index b54f532..dabc397 100644
--- a/assets/NO_FILE.README.md
+++ b/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:
 
diff --git a/assets/NO_FILE_INDELS b/assets/NO_FILE_INDELS
new file mode 100644
index 0000000..e69de29
diff --git a/assets/NO_FILE_MNV b/assets/NO_FILE_MNV
new file mode 100644
index 0000000..e69de29
diff --git a/bin/build_indel_matrix.py b/bin/build_indel_matrix.py
new file mode 100644
index 0000000..5d50b6f
--- /dev/null
+++ b/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())
diff --git a/bin/build_snp_matrix.py b/bin/build_snp_matrix.py
index cc024cc..7b85261 100644
--- a/bin/build_snp_matrix.py
+++ b/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 ""
diff --git a/bin/collect_mnv.py b/bin/collect_mnv.py
new file mode 100644
index 0000000..3170359
--- /dev/null
+++ b/bin/collect_mnv.py
@@ -0,0 +1,109 @@
+#!/usr/bin/env python3
+"""The codon-level changes of a run: one row per sample and codon where reads carry an MNV.
+
+get_MNV groups a sample's SNVs that fall in one codon and reads the combined codon off the reads
+that span it (the BAM). Two SNVs in one codon change one amino acid, and annotating them one base
+at a time names two changes that are neither: CAC>GAG is His>Glu, where each base alone says
+His>Asp and His>Gln, and GAG>TTG is Glu>Leu where the first base alone says a stop.
+
+Its TSV writes the combined amino acid for every codon that holds more than one SNV, also when no
+read carries them together: two SNVs on different molecules of a mixed population (type SNP/MNV,
+0 MNV reads). A codon change is only real when molecules carry it, so a row is kept here when at
+least --min-reads reads carry every change of the codon at once.
+
+    collect_mnv.py --manifest manifest.tsv --min-reads 2 -o <samplesheet>_mnv.tsv
+
+The manifest has one line per sample: sample, reference, get_MNV TSV (the file names stay free).
+Columns: sample, reference, contig, gene, positions, ref_bases, alt_bases, ref_codon, mnv_codon,
+aa_change (the codon read whole), snp_aa_changes (each SNV alone: what a per-base annotation says),
+change_type, consequence_shift (MNV-masked: the SNVs alone name another consequence, such as a stop;
+MNV-gained; Concordant), mnv_reads, total_reads, mnv_frequency, phasing_support (the share of the
+reads spanning the codon with its leading change that carry the whole MNV).
+"""
+from __future__ import annotations
+
+import argparse
+import csv
+import os
+import sys
+
+COLUMNS = ["sample", "reference", "contig", "gene", "positions", "ref_bases", "alt_bases", "ref_codon",
+           "mnv_codon", "aa_change", "snp_aa_changes", "change_type", "consequence_shift", "mnv_reads",
+           "total_reads", "mnv_frequency", "phasing_support"]
+
+
+def _int(v):
+    try:
+        return int(float(v))
+    except (TypeError, ValueError):
+        return None
+
+
+def _clean(v):
+    v = (v or "").strip()
+    return "" if v in ("-", ".") else v
+
+
+def _list(v):
+    """'6139, 6141' -> '6139,6141'."""
+    return ",".join(x.strip() for x in _clean(v).split(",") if x.strip())
+
+
+def read_mnv(path, sample, reference, min_reads):
+    """The rows of one get_MNV TSV whose MNV at least min_reads reads carry."""
+    out = []
+    with open(path, newline="", encoding="utf-8", errors="replace") as fh:
+        for r in csv.DictReader(fh, delimiter="\t"):
+            if "MNV" not in (r.get("Variant Type") or ""):
+                continue
+            reads = _int(r.get("MNV Reads"))
+            if reads is None or reads < min_reads:
+                continue
+            out.append({
+                "sample": sample, "reference": reference, "contig": r.get("Chromosome", ""),
+                "gene": _clean(r.get("Gene")), "positions": _list(r.get("Positions")),
+                "ref_bases": _list(r.get("Reference Bases")), "alt_bases": _list(r.get("Base Changes")),
+                "ref_codon": _clean(r.get("Reference Codon")), "mnv_codon": _clean(r.get("MNV Codon")),
+                "aa_change": _clean(r.get("AA Changes")), "snp_aa_changes": _clean(r.get("SNP AA Changes")),
+                "change_type": _clean(r.get("Change Type")),
+                "consequence_shift": _clean(r.get("MNV Consequence Shift")),
+                "mnv_reads": str(reads), "total_reads": _clean(r.get("Total Reads")),
+                "mnv_frequency": _clean(r.get("MNV Frequencies")),
+                "phasing_support": _clean(r.get("MNV Phasing Support")),
+            })
+    return out
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--manifest", required=True, help="TSV: sample, reference, get_MNV TSV (one line per sample)")
+    ap.add_argument("--min-reads", type=int, default=2, help="reads that must carry every change of the codon")
+    ap.add_argument("-o", "--output", required=True)
+    args = ap.parse_args(argv)
+
+    rows, n_files, n_skipped = [], 0, 0
+    with open(args.manifest, encoding="utf-8") as fh:
+        for line in fh:
+            parts = line.rstrip("\n").split("\t")
+            if len(parts) < 3 or not parts[2]:
+                continue
+            sample, reference, path = parts[0], parts[1], parts[2]
+            if not os.path.exists(path) or os.path.getsize(path) == 0:
+                n_skipped += 1
+                continue
+            n_files += 1
+            rows += read_mnv(path, sample, reference, args.min_reads)
+    rows.sort(key=lambda r: (r["contig"], _int(r["positions"].split(",")[0]) or 0, r["sample"]))
+    with open(args.output, "w", newline="") as out:
+        w = csv.DictWriter(out, fieldnames=COLUMNS, delimiter="\t", lineterminator="\n")
+        w.writeheader()
+        w.writerows(rows)
+    masked = sum(1 for r in rows if r["consequence_shift"] == "MNV-masked")
+    sys.stderr.write(f"[mnv] {len(rows)} phased codon change(s) with >= {args.min_reads} reads in {n_files} "
+                     f"sample(s)" + (f", {masked} of them masked by a per-base annotation" if masked else "")
+                     + (f"; {n_skipped} empty or missing file(s)" if n_skipped else "") + f" -> {args.output}\n")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/bin/qc_report.py b/bin/qc_report.py
index 2510bfd..3bb7317 100644
--- a/bin/qc_report.py
+++ b/bin/qc_report.py
@@ -32,6 +32,7 @@
                               flag_sample, het_frac, is_ancient, lineage_counts_parsed, lineage_fracs, robust)
 from qcreport.panels import build_dynamics, build_epistasis, build_snp_matrix
 from qcreport.parsers import (NBINS, clean_str, consensus_stats, mapdamage_stats, mask_profile, parse_bed,
+                              parse_indel_matrix, parse_mnv_table,
                               parse_collection_dates, parse_dose, parse_dr, parse_gene_burden,
                               parse_gene_conversion, parse_gff,
                               parse_kraken, parse_lineage_colors, parse_metadata, parse_pnps, parse_profile,
@@ -122,6 +123,12 @@ def build_parser():
                     help="The master SNP matrix TSV (build_snp_matrix.py with --depth-vcfs) -> what each "
                          "series gained since its first time point, telling a site absent there from one "
                          "not read (optional).")
+    ap.add_argument("--indel-matrix", default=None,
+                    help="The master indel matrix TSV (build_indel_matrix.py) -> the indel view of the matrix "
+                         "panel (optional).")
+    ap.add_argument("--mnv", default=None,
+                    help="The run's codon-level table (collect_mnv.py over get_MNV) -> which SNPs of the matrix "
+                         "are one amino-acid change together (optional).")
     ap.add_argument("--min-dp", type=int, default=7,
                     help="Depth a site needs before a sample without a call there counts as lacking the "
                          "allele (--consensus_min_dp).")
@@ -436,6 +443,9 @@ def profile_of(sid):
     _dynamics = build_dynamics(series_meta, _variants, _sample_meta)   # feeds dynamics + epistasis
     # front-load 'dose' so it survives the correlation matrix's top-N view
     dist_keys = (["dose"] + DIST) if dose_map else DIST
+    # The SNP view of the matrix panel, with the SNPs that are one codon-level change marked, and the indel
+    # view in the same column order
+    _snp_mx = build_snp_matrix(_variants, provenance.get('reference', ''), mnv=parse_mnv_table(args.mnv))
     now = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M UTC")
     payload = {"generated": now, "version": clean_str(args.version) or "", "repo_url": REPO_URL,
                "counts": counts, "thresholds": thr, "dist": dist_keys,
@@ -456,7 +466,8 @@ def profile_of(sid):
                "n_ancient": n_ancient, "anc_thresholds": anc_thr, "provenance": provenance,
                "dynamics": _dynamics,
                "epistasis": build_epistasis(_dynamics),
-               "snp_matrix": build_snp_matrix(_variants, provenance.get('reference', '')),
+               "snp_matrix": _snp_mx,
+               "indel_matrix": parse_indel_matrix(args.indel_matrix, order=(_snp_mx or {}).get('samples')),
                "coverage": coverage,
                "relatedness": build_relatedness(distances, groups, args.cluster_snps),
                "series": _series,
diff --git a/bin/qcreport/panels.py b/bin/qcreport/panels.py
index 64aaef5..c26756f 100644
--- a/bin/qcreport/panels.py
+++ b/bin/qcreport/panels.py
@@ -248,11 +248,15 @@ def build_epistasis(dynamics, min_r=0.8, min_points=3, top=300, perm=2000, max_v
             'n_strong': sum(1 for p in out if p['tier'] == 'strong')}
 
 
-def build_snp_matrix(variants, reference='', max_sites=50000):
+def build_snp_matrix(variants, reference='', max_sites=50000, mnv=None):
     """Sparse SNP matrix for the report: samples + one row per SNP site with its ref/alt/annotation
     and only the cells (sample index -> [af, dp]) actually called. The payload is gzip-compressed in the
     HTML, so this cap is a high safety backstop (most-shared sites kept) rather than a size limit; the
-    pipeline also writes the complete matrix TSV separately."""
+    pipeline also writes the complete matrix TSV separately.
+
+    `mnv` (parse_mnv_table) marks the SNPs of a sample that are one codon-level change with another
+    one on the same reads: the row gets `mnv` (sample index -> [codon change, fraction of the reads,
+    what the SNPs say alone, consequence shift]) and `mnv_aa`, its distinct codon-level changes."""
     if not variants:
         return None
     samples = list(variants.keys())
@@ -283,5 +287,12 @@ def build_snp_matrix(variants, reference='', max_sites=50000):
              'gene': x['gene'], 'eff': x['eff'], 'aa': x['aa'], 'aa_h37rv': x['aa_h37rv'],
              'gene_h37rv': x['gene_h37rv'], 'pos_h37rv': x['pos_h37rv'], 'n': len(x['cells']), 'cells': x['cells']}
             for x in ordered]
+    if mnv:
+        for r in rows:
+            key = f"{r['contig']}:{r['pos']}"
+            hits = {i: mnv[s][key] for s, i in sidx.items() if i in r['cells'] and key in mnv.get(s, {})}
+            if hits:
+                r['mnv'] = hits
+                r['mnv_aa'] = sorted({h[0] for h in hits.values()})
     return {'samples': samples, 'reference': reference, 'rows': rows,
             'total_sites': len(sites), 'truncated': truncated}
diff --git a/bin/qcreport/parsers.py b/bin/qcreport/parsers.py
index 54e16a3..d9e921c 100644
--- a/bin/qcreport/parsers.py
+++ b/bin/qcreport/parsers.py
@@ -1066,7 +1066,12 @@ def _dyn_opos(info):
 
 def parse_vcfs(paths):
     """{sample: {'chrom:pos': {ref,alt,gene,eff,imp,af}}} from annotated per-sample VCFs (SNPs only).
-    Sample = the VCF #CHROM last column, so filenames are irrelevant."""
+    Sample = the VCF #CHROM last column, so filenames are irrelevant.
+
+    An MNP record (REF and ALT of one length, as FreeBayes writes two changes of one codon it sees on
+    the same reads) is every base it changes, each a SNP at its own position with the record's fraction
+    and annotation, which is the codon read whole: the master SNP matrix reads it the same way. Skipped,
+    those SNPs were missing from every table of the report."""
     out = {}
     for p in paths or []:
         if not p or not os.path.exists(p):
@@ -1088,16 +1093,122 @@ def parse_vcfs(paths):
                     if len(c) < 8:
                         continue
                     chrom, pos, ref, alt = c[0], c[1], c[3], c[4]
-                    if alt in ('.', '') or len(ref) != 1 or any(len(a) != 1 for a in alt.split(',')):
+                    if alt in ('.', ''):
+                        continue
+                    alt1 = alt.split(',')[0]
+                    if len(ref) == 1 and all(len(a) == 1 for a in alt.split(',')):
+                        changes = [(0, ref, alt1)]
+                    elif len(ref) == len(alt1) > 1:
+                        changes = [(k, ref[k], alt1[k]) for k in range(len(ref)) if ref[k].upper() != alt1[k].upper()]
+                    else:
+                        continue   # indels and complex records are not SNPs
+                    try:
+                        start = int(pos)
+                    except ValueError:
                         continue
                     af = _dyn_af(c[8], c[9]) if len(c) >= 10 else 1.0
                     if af is None:
                         continue
                     dp = _dyn_dp(c[8], c[9]) if len(c) >= 10 else None
                     gene, eff, imp, aa = _dyn_ann(c[7])
-                    out[sample][f'{chrom}:{pos}'] = {'ref': ref, 'alt': alt.split(',')[0], 'gene': gene,
-                                                     'eff': eff, 'imp': imp, 'aa': aa, 'aa_h37rv': '',
-                                                     'af': round(af, 4), 'dp': dp, 'opos': _dyn_opos(c[7])}
+                    opos = _dyn_opos(c[7])
+                    for k, r1, a1 in changes:
+                        key = f'{chrom}:{start + k}'
+                        if k and key in out[sample]:
+                            continue   # a record of its own already gave this sample's call here
+                        o = opos
+                        if o and k:
+                            oc, _, op = o.rpartition(':')
+                            o = f'{oc}:{int(op) + k}' if op.isdigit() else o
+                        out[sample][key] = {'ref': r1, 'alt': a1, 'gene': gene,
+                                            'eff': eff, 'imp': imp, 'aa': aa, 'aa_h37rv': '',
+                                            'af': round(af, 4), 'dp': dp, 'opos': o}
         except OSError:
             continue
     return out
+
+
+def _no_file(path):
+    return not path or not os.path.exists(path) or os.path.basename(path).startswith('NO_FILE')
+
+
+def parse_indel_matrix(path, order=None, max_rows=20000):
+    """The indel matrix TSV (build_indel_matrix.py) as the matrix panel holds a matrix: the samples, and
+    one row per indel with its annotation and only the cells a sample calls, [af, dp, filter] (PASS or
+    LowSupport: a call below the rule SNPs are held to, kept because a minority indel starts that way).
+
+    `order` puts the samples in the SNP view's order, so switching views keeps every column in place.
+    The cap keeps the most-called indels, as the SNP view does; the run's TSV has all of them."""
+    if _no_file(path):
+        return None
+    rows, samples = [], []
+    try:
+        with open(path, encoding='utf-8', errors='replace') as fh:
+            header = fh.readline().rstrip('\n').split('\t')
+            col = {h: i for i, h in enumerate(header)}
+            if not {'contig', 'pos', 'ref_allele', 'alt_allele'} <= set(col):
+                return None
+            found = [h[:-3] for h in header if h.endswith('|AF')]
+            samples = [s for s in (order or []) if s in found] + [s for s in found if s not in (order or [])]
+            for line in fh:
+                c = line.rstrip('\n').split('\t')
+                if len(c) < len(header):
+                    continue
+                cells = {}
+                for k, s in enumerate(samples):
+                    try:
+                        af = float(c[col[s + '|AF']])
+                    except (ValueError, KeyError):
+                        continue
+                    if af <= 0:
+                        continue          # read without the indel, or not read: the TSV tells which
+                    dp = c[col[s + '|DP']]
+                    ft = c[col[s + '|FT']] if (s + '|FT') in col else ''
+                    cells[k] = [round(af, 4), int(dp) if dp.isdigit() else None, ft or 'PASS']
+                if not cells:
+                    continue
+                get = lambda k: c[col[k]] if k in col else ''
+                try:
+                    pos, length = int(get('pos')), int(get('length') or 0)
+                except ValueError:
+                    continue
+                rows.append({'contig': get('contig'), 'pos': pos, 'ref': get('ref_allele'), 'alt': get('alt_allele'),
+                             'len': length, 'gene': get('gene'), 'eff': get('effect'), 'hgvs_c': get('hgvs_c'),
+                             'aa': get('hgvs_p'), 'n': len(cells), 'cells': cells})
+    except OSError:
+        return None
+    if not rows:
+        return None
+    total, truncated = len(rows), len(rows) > max_rows
+    if truncated:
+        rows = sorted(rows, key=lambda r: -r['n'])[:max_rows]
+        rows.sort(key=lambda r: (r['contig'], r['pos']))
+    return {'samples': samples, 'rows': rows, 'total_sites': total, 'truncated': truncated}
+
+
+def parse_mnv_table(path):
+    """{sample: {'contig:pos': [codon change, fraction of the reads carrying it, the changes the SNPs
+    name one at a time, consequence shift]}} from the run's codon-level table (collect_mnv.py): which SNP
+    of a sample is half of one amino-acid change with another one on the same reads."""
+    out = {}
+    if _no_file(path):
+        return out
+    import csv
+    try:
+        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
+                try:
+                    freq = round(float(r.get('mnv_frequency') or ''), 4)
+                except ValueError:
+                    freq = None
+                entry = [aa, freq, (r.get('snp_aa_changes') or '').strip(), (r.get('consequence_shift') or '').strip()]
+                for pos in (r.get('positions') or '').split(','):
+                    pos = pos.strip()
+                    if pos:
+                        out.setdefault(r.get('sample', ''), {})[f"{r.get('contig', '')}:{pos}"] = entry
+    except OSError:
+        return {}
+    return out
diff --git a/bin/qcreport/render.py b/bin/qcreport/render.py
index 8c5cb4c..174863c 100644
--- a/bin/qcreport/render.py
+++ b/bin/qcreport/render.py
@@ -177,7 +177,12 @@ def build_html(title, payload):
                  "SNPs, NOT proof of a functional interaction.",
     "snpmatrix": "Every SNP site (rows) by sample (columns); each cell is the allele frequency with its depth, "
                  "plus a reference column and samplesheet metadata as column-header levels. Filter by gene / "
-                 "position or by metadata, and download the full matrix as a TSV.",
+                 "position or by metadata, and download the full matrix as a TSV. A SNP that shares its codon "
+                 "with another one on the same reads carries a 'codon' badge: the amino-acid change is the codon "
+                 "read whole (get_MNV), not the one it names alone, and its cells are outlined. The Indels view "
+                 "holds every indel a sample calls with a PASS, with its length, a frameshift marked, and HGVS; "
+                 "a hatched cell is a call below the rule SNPs are held to, kept because a minority indel starts "
+                 "that way.",
     "drug": "Resistance-associated mutations detected by pathotypr against the WHO catalogue (H37Rv numbering), "
             "as a sample x drug matrix (worst grade per drug) and a per-mutation table with the WHO confidence "
             "grade. Alignment-free (k-mer), so it works regardless of the mapping reference. A genomic screen, "
diff --git a/bin/report_assets/css/04_panels.css b/bin/report_assets/css/04_panels.css
index 19f4866..34b57d2 100644
--- a/bin/report_assets/css/04_panels.css
+++ b/bin/report_assets/css/04_panels.css
@@ -156,6 +156,16 @@ table.snpmx .snpmx-h{writing-mode:vertical-rl;transform:rotate(180deg);font-size
 table.snpmx .snpmx-info{position:sticky;left:0;background:var(--panel);z-index:4;text-align:left;padding:5px 13px;border-right:1px solid var(--line);white-space:nowrap;font-size:12.5px}
 table.snpmx thead .snpmx-info{z-index:6;background:var(--soft)}
 .snpmx-aa{color:var(--accent);font-weight:700}
+/* the matrix's two views, and what an indel or a codon-level change adds to a row */
+.snpmx-mode{margin-right:6px} .snpmx-mode .snpmx-modec{font-size:11px;opacity:.75;margin-left:3px;font-variant-numeric:tabular-nums}
+.snpmx-len{font-size:11px;color:var(--txt2);border:1px solid var(--line);border-radius:6px;padding:0 5px;white-space:nowrap}
+.snpmx-len.fs{color:#b3261e;border-color:rgba(179,38,30,.45);background:rgba(179,38,30,.07);font-weight:600}
+html.dark .snpmx-len.fs{color:#ff8a80;border-color:rgba(255,138,128,.45);background:rgba(255,138,128,.08)}
+.snpmx-mnv{font-size:11px;font-weight:600;color:#7a3db8;border:1px dashed rgba(122,61,184,.55);border-radius:6px;padding:0 5px;white-space:nowrap;cursor:help}
+html.dark .snpmx-mnv{color:#c9a6f0;border-color:rgba(201,166,240,.55)}
+table.snpmx td.snpmx-cell.snpmx-low{opacity:.55;background-image:repeating-linear-gradient(135deg,transparent 0 4px,rgba(255,255,255,.55) 4px 6px)}
+table.snpmx td.snpmx-cell.snpmx-inmnv{box-shadow:inset 0 0 0 2px #7a3db8}
+html.dark table.snpmx td.snpmx-cell.snpmx-inmnv{box-shadow:inset 0 0 0 2px #c9a6f0}
 .snpmx-refn{color:var(--accent);font-weight:700}                        /* the mapping reference name in the matrix corner */
 .snpmx-refn2{color:var(--mut);font-weight:600;font-size:.92em}          /* the reference of interest (H37Rv) */
 table.snpmx td.snpmx-cell{min-width:var(--mxcw,32px);text-align:center;color:#0f2431;font-variant-numeric:tabular-nums;padding:3px 2px}
diff --git a/bin/report_assets/js/12_snpmatrix.js b/bin/report_assets/js/12_snpmatrix.js
index 120c017..67573a6 100644
--- a/bin/report_assets/js/12_snpmatrix.js
+++ b/bin/report_assets/js/12_snpmatrix.js
@@ -5,11 +5,18 @@ var SNPMX_PAL_DARK=['#2f4a6b','#6b4c2f','#2f5a40','#6b3a4a','#4a3a6b','#2f5a55',
 var snpmxFilter={};
 var snpmxAll=false;   // false = show the top MAXR most-shared sites; true = virtualized scroll over every site
 var snpmxZoom=32;     // SNP-matrix cell width in px (zoom slider); CSS var --mxcw, no re-render (row height untouched -> virtualization safe)
+var snpmxMode='snp';  // which variants the matrix shows: 'snp' or 'indel' (the run's indel matrix)
+// an indel's allele, short enough for a row label: long ones keep their ends and say how long they are
+function mxAllele(a){ a=String(a||''); return a.length>9?(esc(a.slice(0,4))+'&#8230;'+esc(a.slice(-2))):esc(a); }
 function renderSnpMatrix(){
   var host=el('snpmx_body'), sec=el('snpmatrix'); if(!host)return;
-  var M=R.snp_matrix;
-  if(!(M&&M.rows&&M.rows.length)){ if(sec)sec.style.display='none'; var nv=el('nav-snpmx'); if(nv)nv.style.display='none'; return; }
+  var SM=R.snp_matrix, IM=R.indel_matrix;
+  var hasS=!!(SM&&SM.rows&&SM.rows.length), hasI=!!(IM&&IM.rows&&IM.rows.length);
+  if(!hasS&&!hasI){ if(sec)sec.style.display='none'; var nv=el('nav-snpmx'); if(nv)nv.style.display='none'; return; }
   if(sec)sec.style.display='';
+  if(snpmxMode==='indel'&&!hasI) snpmxMode='snp';
+  if(snpmxMode==='snp'&&!hasS) snpmxMode='indel';
+  var isIndel=snpmxMode==='indel', M=isIndel?IM:SM, unit=isIndel?'indel(s)':'SNP site(s)';
   var samples=M.samples, MAXR=400;
   var meta=(R.sample_meta&&R.sample_meta.fields&&R.sample_meta.fields.length)?R.sample_meta:null;
   var metaMaps={}, metaVals={};
@@ -34,9 +41,13 @@ function renderSnpMatrix(){
     meta.fields.map(function(f){ return '<label class="snpmx-fsel">'+esc(f)+' <select data-f="'+esc(f)+'"><option value="">all</option>'+
       metaVals[f].map(function(v){return '<option value="'+esc(v)+'"'+(snpmxFilter[f]===v?' selected':'')+'>'+esc(v)+'</option>';}).join('')+'</select></label>'; }).join('')+
     '<button class="dyn-btn" id="snpmxfclear">clear</button></div>'):'';
+  // SNPs or indels: each view is its own matrix, in the same sample columns
+  var modeUI=(hasS&&hasI)?('<div class="seg snpmx-mode" role="group" aria-label="variants shown">'+
+      '<button type="button" data-mode="snp" class="'+(isIndel?'':'on')+'" aria-pressed="'+(!isIndel)+'">SNPs <span class="snpmx-modec">'+SM.rows.length+'</span></button>'+
+      '<button type="button" data-mode="indel" class="'+(isIndel?'on':'')+'" aria-pressed="'+isIndel+'">Indels <span class="snpmx-modec">'+IM.rows.length+'</span></button></div>'):'';
   host.innerHTML=
-    '<div class="snpmx-controls">'+
-      '<input id="snpmxq" class="dyn-search" type="search" placeholder="filter by gene / position / amino acid...">'+
+    '<div class="snpmx-controls">'+modeUI+
+      '<input id="snpmxq" class="dyn-search" type="search" placeholder="filter by gene / position / '+(isIndel?'HGVS':'amino acid')+'...">'+
       '<label class="snpmx-toggle"><input type="checkbox" id="snpmxdp" checked> show depth</label>'+
       '<button class="dyn-btn" id="snpmxdl" title="Download the matrix as the report holds it (all samples, called cells only) as a wide TSV. The run&#39;s own snp_matrix.tsv also gives AF 0 and the depth where a sample was read without the allele">'+icon('download')+'download matrix (TSV)</button>'+
       '<button class="dyn-btn snpmx-allbtn" id="snpmxallbtn" style="display:none" title="Toggle between the top most-shared sites and a scrollable view of every site"></button>'+
@@ -54,7 +65,7 @@ function renderSnpMatrix(){
     });
     var nfilt=0; for(var kf in snpmxFilter){ if(snpmxFilter[kf]) nfilt++; }
     var total=rows.length, capped=total>MAXR, ncol=1+vi.length;
-    el('snpmxcount').innerHTML=total+' SNP site(s) &#215; '+vi.length+' sample(s)'+(nfilt?' (filtered)':'')+
+    el('snpmxcount').innerHTML=total+' '+unit+' &#215; '+vi.length+' sample(s)'+(nfilt?' (filtered)':'')+
       (capped?(' &#183; showing <b>'+(snpmxAll?('all '+total):(MAXR+' of '+total))+'</b>'):'')+(M.truncated?' &#183; full matrix in the TSV':'');
     var allbtn=el('snpmxallbtn');   // a real button toggles between the top sites and the full virtualized scroll
     if(allbtn){ if(capped){ allbtn.style.display=''; allbtn.classList.toggle('on',snpmxAll);
@@ -67,16 +78,31 @@ function renderSnpMatrix(){
         vi.map(function(i){var s=samples[i],v=(meta.rows[s]||{})[f]||'',bg=metaColor(f,v); return '<th class="snpmx-metacell" style="top:'+(k*mh)+'px;background:'+bg+';color:'+metaText(bg)+'" title="'+esc(f)+': '+esc(v||'-')+'">'+esc(v)+'</th>';}).join('')+'</tr>';
     }).join(''):'';
     var stop=nf*mh;
-    var refhdr='SNP'+((M.reference||REFNAME)?(' <span class="snpmx-refn" title="reference the samples were mapped against">'+esc(M.reference||REFNAME)+'</span>'):'')+(AA2LBL?(' <span class="snpmx-refn2" title="reference of interest - amino acids (and the position, when the pipeline lifts it) are also shown in this reference&#39;s coordinates in brackets when they differ">&#8596; '+esc(AA2LBL)+'</span>'):'');
+    var refhdr=(isIndel?'Indel':'SNP')+((M.reference||REFNAME)?(' <span class="snpmx-refn" title="reference the samples were mapped against">'+esc(M.reference||REFNAME)+'</span>'):'')+(AA2LBL?(' <span class="snpmx-refn2" title="reference of interest - amino acids (and the position, when the pipeline lifts it) are also shown in this reference&#39;s coordinates in brackets when they differ">&#8596; '+esc(AA2LBL)+'</span>'):'');
     var nameRow='<tr><th class="snpmx-info snpmx-corner" style="top:'+stop+'px">'+refhdr+'</th>'+
       vi.map(function(i){var s=samples[i];return '<th class="snpmx-hcell" style="top:'+stop+'px" title="'+esc(s)+'"><span class="snpmx-h">'+esc(s)+'</span></th>';}).join('')+'</tr>';
     function rowHTML(r){
-      var lbl='<b>'+esc(r.gene||r.contig)+'</b>'+geneRvTag(r.gene)+' '+refPos(r.pos,r.pos_h37rv)+' '+esc(r.ref)+'&#8594;'+esc(r.alt)+(r.aa?(' <span class="snpmx-aa">'+aaDual(r.aa,r.aa_h37rv)+'</span>'):'');
+      var lbl;
+      if(isIndel){
+        // an indel: its length (a frameshift when it is not a multiple of three in a gene) and HGVS
+        var fs=/frameshift/.test(r.eff||''), ln=(r.len>0?'+':'&#8722;')+Math.abs(r.len)+' bp';
+        lbl='<b>'+esc(r.gene||r.contig)+'</b>'+geneRvTag(r.gene)+' '+refPos(r.pos)+' '+mxAllele(r.ref)+'&#8594;'+mxAllele(r.alt)+
+          ' <span class="snpmx-len'+(fs?' fs':'')+'" title="'+esc((r.eff||'').replace(/_/g,' ')+(r.hgvs_c?(' · '+r.hgvs_c):''))+'">'+ln+(fs?' &#183; frameshift':'')+'</span>'+
+          (r.aa?(' <span class="snpmx-aa">'+esc(r.aa)+'</span>'):'');
+      } else {
+        lbl='<b>'+esc(r.gene||r.contig)+'</b>'+geneRvTag(r.gene)+' '+refPos(r.pos,r.pos_h37rv)+' '+esc(r.ref)+'&#8594;'+esc(r.alt)+(r.aa?(' <span class="snpmx-aa">'+aaDual(r.aa,r.aa_h37rv)+'</span>'):'');
+        // this SNP and another one of its codon are on the same reads: the codon read whole is the change
+        if(r.mnv_aa) lbl+=' <span class="snpmx-mnv" title="Codon-level change: in '+Object.keys(r.mnv).length+' sample(s) this SNP and another one of its codon are on the same reads (get_MNV), so the amino-acid change is the codon read whole: '+esc(r.mnv_aa.join(', '))+(r.aa?(', not '+esc(r.aa)+' alone'):'')+'">codon '+esc(r.mnv_aa.join(', '))+'</span>';
+      }
       var cells=vi.map(function(i){var c=r.cells[i], s=samples[i];
-        if(!c) return '<td class="snpmx-cell snpmx-empty" title="'+esc(s)+' - no call here. The run&#39;s snp_matrix.tsv says whether it was read: AF 0 with the depth, or empty with no read"></td>';
+        if(!c) return '<td class="snpmx-cell snpmx-empty" title="'+esc(s)+' - no call here. The run&#39;s '+(isIndel?'indel':'snp')+'_matrix.tsv says whether it was read: AF 0 with the depth, or empty with no read"></td>';
         var afTxt=c[0].toFixed(2).replace(/^0/,'').replace(/^1\.00$/,'1');
         var dpTxt=(showDP&&c[1]!=null)?('<span class="snpmx-dp">'+c[1]+'</span>'):'';
-        return '<td class="snpmx-cell'+(showDP?' wdp':'')+'" style="background:'+snpAfColor(c[0])+'" title="'+esc(s)+'  AF='+c[0].toFixed(3)+(c[1]!=null?('  DP='+c[1]):'')+'"><span class="snpmx-af">'+afTxt+'</span>'+dpTxt+'</td>';
+        var low=isIndel&&c[2]&&c[2]!=='PASS', mv=!isIndel&&r.mnv&&r.mnv[i];
+        var tip=esc(s)+'  AF='+c[0].toFixed(3)+(c[1]!=null?('  DP='+c[1]):'')+
+          (low?('  - called below the rule SNPs are held to ('+esc(c[2])+'): kept, because a minority indel starts that way'):'')+
+          (mv?('  - codon-level change '+esc(mv[0])+(mv[1]!=null?(', on '+Math.round(mv[1]*100)+'% of the reads'):'')+(mv[2]?(' (alone: '+esc(mv[2])+')'):'')+(mv[3]==='MNV-masked'?' - the SNPs alone would name another consequence':'')):'');
+        return '<td class="snpmx-cell'+(showDP?' wdp':'')+(low?' snpmx-low':'')+(mv?' snpmx-inmnv':'')+'" style="background:'+snpAfColor(c[0])+'" title="'+tip+'"><span class="snpmx-af">'+afTxt+'</span>'+dpTxt+'</td>';
       }).join('');
       return '<tr><td class="snpmx-info">'+lbl+'</td>'+cells+'</tr>';
     }
@@ -113,15 +139,19 @@ function renderSnpMatrix(){
     el('snpmxfclear').onclick=function(){ snpmxFilter={}; Array.prototype.forEach.call(host.querySelectorAll('.snpmx-filters select'),function(s){s.value='';}); draw(); };
   }
   el('snpmxdl').onclick=function(){
-    var hdr=['reference','contig','pos','pos_'+AA2LBL,'ref_allele','alt_allele','gene','effect','aa_change'];
-    samples.forEach(function(s){hdr.push(s+'|AF');hdr.push(s+'|DP');});
+    var hdr=isIndel?['reference','contig','pos','ref_allele','alt_allele','length','gene','effect','hgvs_c','hgvs_p']
+                   :['reference','contig','pos','pos_'+AA2LBL,'ref_allele','alt_allele','gene','effect','aa_change','codon_change'];
+    var per=isIndel?3:2;
+    samples.forEach(function(s){hdr.push(s+'|AF');hdr.push(s+'|DP');if(isIndel)hdr.push(s+'|FT');});
     var lines=[hdr.join('\t')];
-    if(meta){ meta.fields.forEach(function(f){ var row=['# '+f,'','','','','','','']; samples.forEach(function(s){row.push((meta.rows[s]||{})[f]||'');row.push('');}); lines.push(row.join('\t')); }); }
-    M.rows.forEach(function(r){var row=[M.reference||r.contig,r.contig,r.pos,(r.pos_h37rv?(''+r.pos_h37rv).split(':').pop():''),r.ref,r.alt,r.gene,r.eff,r.aa];
-      samples.forEach(function(s,i){var c=r.cells[i]; if(c){row.push(c[0].toFixed(4));row.push(c[1]==null?'':c[1]);}else{row.push('');row.push('');}});
+    if(meta){ meta.fields.forEach(function(f){ var row=['# '+f]; while(row.length<hdr.length-samples.length*per)row.push(''); samples.forEach(function(s){row.push((meta.rows[s]||{})[f]||''); for(var k=1;k<per;k++)row.push('');}); lines.push(row.join('\t')); }); }
+    M.rows.forEach(function(r){var row=isIndel?[M.reference||r.contig,r.contig,r.pos,r.ref,r.alt,r.len,r.gene,r.eff,r.hgvs_c||'',r.aa]
+                                              :[M.reference||r.contig,r.contig,r.pos,(r.pos_h37rv?(''+r.pos_h37rv).split(':').pop():''),r.ref,r.alt,r.gene,r.eff,r.aa,(r.mnv_aa||[]).join(';')];
+      samples.forEach(function(s,i){var c=r.cells[i]; if(c){row.push(c[0].toFixed(4));row.push(c[1]==null?'':c[1]);if(isIndel)row.push(c[2]||'');}else{row.push('');row.push('');if(isIndel)row.push('');}});
       lines.push(row.join('\t'));});
-    dl(lines.join('\n')+'\n','snp_matrix.tsv','text/tab-separated-values');
+    dl(lines.join('\n')+'\n',isIndel?'indel_matrix.tsv':'snp_matrix.tsv','text/tab-separated-values');
   };
+  Array.prototype.forEach.call(host.querySelectorAll('.snpmx-mode button'),function(b){ b.onclick=function(){ var m=b.getAttribute('data-mode'); if(m!==snpmxMode){ snpmxMode=m; snpmxAll=false; renderSnpMatrix(); } }; });
   window.__snpmxDraw=draw;   // let the header 'all sites' shortcut repaint the matrix without a full re-render
   draw();
   if(window.__syncSitesBtn)window.__syncSitesBtn();
diff --git a/bin/report_assets/shell.html b/bin/report_assets/shell.html
index 4101d9d..a4c09df 100644
--- a/bin/report_assets/shell.html
+++ b/bin/report_assets/shell.html
@@ -10,7 +10,7 @@
 <div class="toc-group" data-page="qc"><a class="toc-gh" href="#p-qc"><span class="toc-n">2</span>Sample QC<span class="toc-badge" id="bdg-qc"></span></a><div class="toc-items"><a class="toc-link" href="#flagged">Flagged samples</a><a class="toc-link" href="#gstats">All samples</a><a class="toc-link" href="#linsum" id="nav-lin">Lineages</a><a class="toc-link" href="#kraken" id="nav-kraken">Contamination</a><a class="toc-link" href="#dist">Distributions</a><a class="toc-link" href="#adna" id="nav-adna">aDNA damage</a></div></div>
 <div class="toc-group" data-page="genome"><a class="toc-gh" href="#p-genome"><span class="toc-n">3</span>Genome &amp; genes</a><div class="toc-items"><a class="toc-link" href="#cons">Consensus</a><a class="toc-link" href="#genome" id="nav-genome">Genome landscape</a><a class="toc-link" href="#deletions" id="nav-del">Deletions</a><a class="toc-link" href="#function" id="nav-function">Functional impact</a><a class="toc-link" href="#geneburden" id="nav-geneburden">Gene burden</a><a class="toc-link" href="#hotspots" id="nav-hot">Variable genes</a><a class="toc-link" href="#pnps" id="nav-pnps">Selection (dN/dS)</a></div></div>
 <div class="toc-group" data-page="related"><a class="toc-gh" href="#p-related"><span class="toc-n">4</span>Relatedness<span class="toc-badge" id="bdg-related"></span></a><div class="toc-items"><a class="toc-link" href="#reldist" id="nav-reldist">Distances</a><a class="toc-link" href="#relclus" id="nav-relclus">Clusters</a><a class="toc-link" href="#relgroup" id="nav-relgroup">Outside their group</a></div></div>
-<div class="toc-group" data-page="variants"><a class="toc-gh" href="#p-variants"><span class="toc-n">5</span>Variants over time<span class="toc-badge" id="bdg-variants"></span></a><div class="toc-items"><a class="toc-link" href="#gains" id="nav-gains">Gained per series</a><a class="toc-link" href="#minority" id="nav-minority">Minority variants</a><a class="toc-link" href="#dynamics" id="nav-dyn">SNP dynamics</a><a class="toc-link" href="#epistasis" id="nav-epi">Co-varying pairs</a><a class="toc-link" href="#snpmatrix" id="nav-snpmx">SNP matrix</a><a class="toc-link" href="#vardose" id="nav-vardose">Variant &#215; dose</a></div></div>
+<div class="toc-group" data-page="variants"><a class="toc-gh" href="#p-variants"><span class="toc-n">5</span>Variants over time<span class="toc-badge" id="bdg-variants"></span></a><div class="toc-items"><a class="toc-link" href="#gains" id="nav-gains">Gained per series</a><a class="toc-link" href="#minority" id="nav-minority">Minority variants</a><a class="toc-link" href="#dynamics" id="nav-dyn">SNP dynamics</a><a class="toc-link" href="#epistasis" id="nav-epi">Co-varying pairs</a><a class="toc-link" href="#snpmatrix" id="nav-snpmx">Variant matrix</a><a class="toc-link" href="#vardose" id="nav-vardose">Variant &#215; dose</a></div></div>
 <div class="toc-group" data-page="drug"><a class="toc-gh" href="#p-drug"><span class="toc-n">6</span>Resistance<span class="toc-badge" id="bdg-drug"></span></a><div class="toc-items"><a class="toc-link" href="#drug" id="nav-drug">Resistance mutations</a><a class="toc-link" href="#drseries" id="nav-drseries">Along each series</a></div></div>
 <div class="toc-group" data-page="gconv"><a class="toc-gh" href="#p-gconv"><span class="toc-n">7</span>Gene conversion<span class="toc-badge" id="bdg-gconv"></span></a><div class="toc-items"><a class="toc-link" href="#gconv" id="nav-gconv">Candidate tracts</a></div></div>
 <div class="toc-group" data-page="diag"><a class="toc-gh" href="#p-diag"><span class="toc-n">8</span>Diagnostics</a><div class="toc-items"><a class="toc-link" href="#corr">Metric pairs</a><a class="toc-link" href="#corrmatrix">Correlation matrix</a><a class="toc-link" href="#qcpca" id="nav-pca">QC space</a><a class="toc-link" href="#divcomp" id="nav-divcomp">Divergence</a><a class="toc-link" href="#dosetx" id="nav-dosetx">Dose &#215; treatment</a><a class="toc-link" href="#temporal" id="nav-temporal">Sampling dates</a></div></div>
@@ -123,7 +123,7 @@
 <div class="panel pad" id="dyn_body"></div></section>
 <section id="epistasis"><h2>Co-varying pairs <span class="c">variants whose trajectories move together (concordant) or against each other (discordant) within a series</span></h2>
 <div class="panel pad" id="epi_body"></div></section>
-<section id="snpmatrix"><h2>SNP matrix <span class="c">every SNP site &#215; sample, as allele frequency; filter by gene or position and download the full matrix</span></h2>
+<section id="snpmatrix"><h2>Variant matrix <span class="c">every SNP site, or every indel, &#215; sample, as allele frequency; filter by gene or position and download the full matrix</span></h2>
 <div class="panel pad" id="snpmx_body"></div></section>
 <section id="vardose"><h2>Variant &#215; dose <span class="c">per-variant allele frequency against dose (Spearman &#961; with Benjamini&#8211;Hochberg FDR); click a row to plot it</span>
 <input class="gsearch" id="vardoseq" type="search" placeholder="search gene / position / AA" style="margin-left:auto"><button class="exp-h" data-panel="vardosePanel" data-render="vardose"><span data-ic="maximize"></span>full</button></h2>
diff --git a/docs/configuration.md b/docs/configuration.md
index 2014028..3731582 100644
--- a/docs/configuration.md
+++ b/docs/configuration.md
@@ -76,6 +76,8 @@ table below.
 | :--- | :---: | :--- |
 | `--make_qc_report` | :octicons-check-circle-fill-16:{ .green } on | The interactive HTML QC report + `qc_flags.tsv` |
 | `--make_snp_matrix` | :octicons-check-circle-fill-16:{ .green } on | The cohort master SNP-matrix TSV |
+| `--make_indel_matrix` | :octicons-check-circle-fill-16:{ .green } on | The cohort master indel-matrix TSV |
+| `--run_mnv` | :octicons-check-circle-fill-16:{ .green } on | Codon-level changes with get_MNV (per-sample `mnv/`, `<samplesheet>_mnv.tsv`, the SNP matrix's `codon_change`) |
 | `--make_consensus` | :octicons-check-circle-fill-16:{ .green } on | The per-sample consensus FASTA |
 | `--exclude_repeats` | :octicons-check-circle-fill-16:{ .green } on | `nucmer` self-alignment repeat masking |
 | `--dynamic_read_filter` | :octicons-check-circle-fill-16:{ .green } on | The length-aware (mappability) read filter |
@@ -196,6 +198,12 @@ Two documented, ready-to-use **organism** configs live in
 | | `--het_min_frac` | `0.10` | Frequency between 0.10 and 0.90 is **Heterozygous**. |
 | | `--filter_min_dp` | `30` | Minimum depth required to call a variant. |
 | | `--min_alt_fwd/rev` | `2` | Min variant supporting reads in FWD and REV strands. |
+| **Indels & codons** | `--make_indel_matrix` | `true` | Emit `<samplesheet>_indel_matrix.tsv`: every indel a sample calls with a PASS (the rules above), each sample's AF, depth and filter. |
+| | `--run_mnv` | `true` | Read the SNPs that share a codon whole, on the reads that carry them ([get_MNV](https://github.com/PathoGenOmics-Lab/get_MNV)): the amino-acid change a per-base annotation misses. |
+| | `--mnv_min_reads` | `2` | Reads that must carry every change of a codon for it to count as one change in `<samplesheet>_mnv.tsv`. |
+| | `--mnv_gff_features` | `CDS` | GFF feature types get_MNV reads codons from. |
+| | `--mnv_translation_table` | `11` | NCBI genetic code get_MNV translates with (11: bacteria and archaea). |
+| | `--mnv_container` | `quay.io/biocontainers/get_mnv` 1.1.5, by digest | The image get_MNV runs in until the pipeline image carries it; empty runs it in the pipeline image. |
 | **Consensus** | `--make_consensus` | `true` | Generate a consensus FASTA for each sample. |
 | | `--consensus_min_dp` | `7` | Depth threshold below which a base becomes "No Call". |
 | | `--consensus_mask_char` | `X` | Character for masked/low-quality sites. |
diff --git a/docs/outputs.md b/docs/outputs.md
index dfb2387..ca517d2 100644
--- a/docs/outputs.md
+++ b/docs/outputs.md
@@ -14,6 +14,8 @@ results_bampiro/
 ├── samples_legio_qc_flags.tsv            # Per-sample PASS/WARN/FAIL verdicts
 ├── samples_legio_summary.tsv             # Cohort metrics table (feeds the report)
 ├── samples_legio_snp_matrix.tsv          # Master SNP matrix (site × sample AF & depth)
+├── samples_legio_indel_matrix.tsv        # Master indel matrix (indel × sample AF, depth & filter; see below)
+├── samples_legio_mnv.tsv                 # 🧩 Codon-level changes: SNPs of one codon on the same reads (get_MNV; see below)
 ├── samples_legio_deletions.tsv           # Stretches some samples have no reads for (see below)
 ├── samples_legio_snp_distances.tsv       # SNPs between every two consensus sequences (square matrix)
 ├── samples_legio_snp_distances_pairs.tsv # The same, one row per pair, with how much each rests on
@@ -34,8 +36,10 @@ results_bampiro/
     ├── MP00091.LENS.filtered.bam.bai       # BAM Index  (or .cram with --output_cram)
     │                                       # (final.bam is published instead when dynamic_read_filter=false, or in addition with --publish_prefilter_bam)
     │
-    ├── MP00091.LENS.ann.vcf.gz             # 🎯 MAIN OUTPUT: Annotated Variants (SNPs/Indels)
+    ├── MP00091.LENS.ann.vcf.gz             # 🎯 MAIN OUTPUT: Annotated SNPs (one record per base)
     ├── MP00091.LENS.ann.vcf.gz.tbi         # Index for the main VCF
+    ├── MP00091.LENS.indels.ann.vcf.gz      # 🎯 Annotated indels: every indel call, FILTER PASS or LowSupport (see below)
+    ├── MP00091.LENS.indels.ann.vcf.gz.tbi  # Index for the indel VCF
     ├── 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)
@@ -47,6 +51,12 @@ results_bampiro/
     ├── MP00091.LENS.var.homo.SNPs.ann...   # 📂 Split VCFs: Subset of Homozygous SNPs (Annotated)
     ├── MP00091.LENS.var.het.SNPs.ann...    # 📂 Split VCFs: Subset of Heterozygous SNPs (Annotated)
     ├── MP00091.LENS.var.homo.indel...      # 📂 SPLIT VCF: Homozygous Indels only
+    ├── MP00091.LENS.var.het.indel...       # 📂 SPLIT VCF: Heterozygous (minority) Indels only
+    │
+    ├── mnv/                                # 🧩 Codon-level changes (get_MNV; only if --run_mnv, on by default)
+    │   ├── MP00091.LENS.MNV.tsv            # -> One row per codon holding a call: the codon read whole, its reads
+    │   ├── MP00091.LENS.MNV.vcf.gz         # -> The same as VCF: a record per SNP and one per MNV (INFO/MR = its reads; bgzipped, not indexed)
+    │   └── MP00091.LENS.mnv.summary.json / .manifest.json   # -> Counts, command, version and checksums of the run
     │
     ├── lineage/                            # 🧬 Pathotypr typing (only if --run_pathotypr)
     │   └── MP00091__<runId>.pathotypr.lineage_summary.tsv   # -> Pathotypr lineage / sub-lineage call (__<runId> avoids multi-lane collisions)
@@ -97,6 +107,55 @@ genotype's dosage passed off as a fraction.
 Rows are keyed on contig and position, so two references that share a contig name cannot be told
 apart; the matrix step warns when it sees one.
 
+With the codon-level table (`--run_mnv`, on by default) two columns follow `aa_change`:
+`codon_change`, the amino-acid change of every codon the SNP is part of together with another SNP
+on the same reads, and `codon_change_samples`, the samples whose reads carry it and 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 by a run without the table. Where they are filled,
+`aa_change` is what the SNP names alone, which is not the change the codon undergoes.
+
+## Indel matrix
+
+`<samplesheet>_indel_matrix.tsv` holds the run's indels the way the SNP matrix holds its SNPs.
+Indels are called with the definition a SNP has (the depth, strand support and allele fraction of
+`--filter_min_dp`, `--min_alt_fwd`/`--min_alt_rev`, `--het_min_frac` and `--hom_threshold`), and
+each sample's `<sample>.<ref>.indels(.ann).vcf.gz` keeps every indel FreeBayes made, normalised
+(left-aligned, one allele per record): the ones that meet the rule have `FILTER` `PASS` and the
+rest `LowSupport`.
+
+A row is an indel (contig, position, REF, ALT) that at least one sample calls with a `PASS`, with its
+`length` (positive for an insertion), `gene`, `effect`, `hgvs_c`, `hgvs_p` and `n_pass` (how many
+samples call it with a `PASS`), then three columns per sample:
+
+| `AF` | `DP` | `FT` | Meaning |
+| :--- | :--- | :--- | :--- |
+| a fraction | the depth | `PASS` | The sample calls the indel and the call meets the rule |
+| a fraction | the depth | `LowSupport` | The sample calls it below the rule: kept, because a minority frameshift starts that way |
+| `0` | the depth | empty | Read at the indel's anchor base, and no such indel called |
+| empty | `0` | empty | No read covers the anchor base in that sample |
+| empty | empty | empty | Not in that sample's genome, or no all-positions VCF |
+
+The depth of a `0` cell is the reads that show a base at the anchor position, so for a long indel it
+is a floor: a read that ends inside the event is not counted. Indels are not lifted to canonical
+(H37Rv) coordinates.
+
+## Codon-level changes (MNVs)
+
+Two SNPs in one codon change one amino acid, and annotated one base at a time they name two changes
+that are not there: CAC>GAG is His>Glu, where each base alone says His>Asp and His>Gln, and GAG>TTG
+is Glu>Leu, where the first base alone says a stop. [get_MNV](https://github.com/PathoGenOmics-Lab/get_MNV)
+groups each sample's passing SNPs and indels by codon and reads the codon whole, on the reads that
+span it (with the base and mapping-quality floors FreeBayes called with).
+
+Its per-sample table writes a combined change for every codon holding two calls, also when no read
+carries them together (two SNPs on different molecules of a mixed population: `MNV Reads` 0).
+`<samplesheet>_mnv.tsv` keeps a codon only where at least `--mnv_min_reads` reads carry every change
+at once: `sample`, `reference`, `contig`, `gene`, `positions`, `ref_bases`, `alt_bases`, `ref_codon`,
+`mnv_codon`, `aa_change` (the codon read whole), `snp_aa_changes` (what each SNP says alone),
+`change_type`, `consequence_shift` (`MNV-masked` when the SNPs alone name another consequence, such
+as a stop), `mnv_reads`, `total_reads`, `mnv_frequency` and `phasing_support`. The same codons mark
+the SNP matrix (`codon_change`) and the report's variant matrix.
+
 ## SNP distances
 
 `<samplesheet>_snp_distances.tsv` is the square matrix of SNPs between every two consensus
diff --git a/docs/qc-report.md b/docs/qc-report.md
index 69cba3d..f43430a 100644
--- a/docs/qc-report.md
+++ b/docs/qc-report.md
@@ -30,7 +30,7 @@ on *Sample QC*, the samples with resistance mutations their lineage does not sha
 | **2 · Sample QC** | Which samples can be trusted: the verdicts, the live thresholds, the **flagged samples** with the value behind every flag, the **exclusion list** and its exports, the table of all samples, lineages, **contamination** (Kraken2), distributions and aDNA damage. |
 | **3 · Genome & genes** | Consensus completeness, the genome landscape (missing calls, **deletions**, SNPs per bin and **per callable kb**, het, indels), the **deletions** found by comparing each sample's stretches without reads with the rest of the cohort, functional impact (snpEff), gene burden, variable genes and dN/dS. A cohort mapped against several references is shown one reference at a time, each with its own samples, length, genes and masked regions. |
 | **4 · Relatedness** | SNP distances between the consensus sequences, over the positions both samples called: a heatmap in the order that keeps each cluster together, the **clusters** at a threshold you can move, and the samples far from the rest of their own group (patient, line, series). |
-| **5 · Variants over time** | What each series **gained since its first time point** (new, risen, unknown, lost, charted per treatment), the **minority variants** (calls below fixation, the reads they rest on and, with libraries of the same DNA, how often another library reproduces them), SNP dynamics, co-varying pairs (epistasis), the SNP matrix and variant &#215; dose. |
+| **5 · Variants over time** | What each series **gained since its first time point** (new, risen, unknown, lost, charted per treatment), the **minority variants** (calls below fixation, the reads they rest on and, with libraries of the same DNA, how often another library reproduces them), SNP dynamics, co-varying pairs (epistasis), the variant matrix (SNPs, with the ones that are one codon-level change marked, or indels) and variant &#215; dose. |
 | **6 · Resistance** | WHO-catalogue mutations grouped by mutation, the sample &#215; drug matrix, every call, and each series' mutations **per time point**, those acquired since the first one set apart. |
 | **7 · Gene conversion** | Candidate tracts and the evidence behind each verdict. |
 | **8 · Diagnostics** | Metric pairs, the correlation matrix, QC space, divergence vs completeness, dose &#215; treatment and sampling dates: views for digging into a problem, none of which flags a sample. |
diff --git a/main.nf b/main.nf
index 5712ff5..2782b7e 100644
--- a/main.nf
+++ b/main.nf
@@ -16,6 +16,7 @@ include { READ_QC }            from './subworkflows/read_qc'              // 2-3
 include { LINEAGE_TYPING }     from './subworkflows/lineage_typing'       // 4. Pathotypr
 include { MAP_READS }          from './subworkflows/map_reads'            // 5. Mapping, merge, read filter
 include { CALL_VARIANTS }      from './subworkflows/call_variants'        // 6. Variant calling
+include { CODON_CHANGES }      from './subworkflows/codon_changes'        // 6b. Codon-level changes (get_MNV)
 include { MAKE_CONSENSUS }     from './subworkflows/make_consensus'       // 7. Consensus
 include { ANNOTATE }           from './subworkflows/annotate'             // 8. Annotation
 include { LEGACY_STATS }       from './subworkflows/legacy_stats'         // 9. Legacy stats
@@ -267,6 +268,9 @@ Run Pathotypr    : ${params.run_pathotypr}
     // 6. Variant Calling
     def variants = CALL_VARIANTS(bams.variant_base)
 
+    // 6b. Codon-level changes: get_MNV reads the SNPs of one codon whole, on the reads that carry them
+    def codons = CODON_CHANGES(variants.codon_calls, bams.variant_base, refGffMap)
+
     // 12. Gene conversion (opt-in). Reads the PRE-FILTER bam and the reference self-alignment,
     // because the masking the rest of the pipeline applies removes exactly this signal.
     def gconv = GENE_CONVERSION(refs.delta, bams.dedup_bam,
@@ -278,8 +282,8 @@ Run Pathotypr    : ${params.run_pathotypr}
                                    variants.allpos, variants.snps, variants.mask_sites)
 
     // 8. Annotation
-    def annotated = ANNOTATE(variants.homo_snp, variants.het_snp, variants.homo_indel,
-                             variants.raw_fb, variants.main_vcf, variants.allpos, refs.snpeff)
+    def annotated = ANNOTATE(variants.homo_snp, variants.het_snp, variants.homo_indel, variants.het_indel,
+                             variants.indels, variants.raw_fb, variants.main_vcf, variants.allpos, refs.snpeff)
 
     // 9. Legacy Stats Generation
     def legacy = LEGACY_STATS(refs.bundle, reads.pe_json, reads.se_json,
@@ -292,5 +296,5 @@ Run Pathotypr    : ${params.run_pathotypr}
     // 11. Cohort-level outputs (QC report + master SNP matrix).
     COHORT_REPORT(legacy.legacy_log, consensus.fasta, refs.bundle, reads.kraken_reports,
                   typing.dr_mutations, annotated.stats_vcf, annotated.freebayes_vcf,
-                  variants.allpos, gconv.cohort_tsv, refMap, refGffMap, tsv_name)
+                  variants.allpos, annotated.indels_vcf, codons.tsv, gconv.cohort_tsv, refMap, refGffMap, tsv_name)
 }
diff --git a/modules/mnv.nf b/modules/mnv.nf
new file mode 100644
index 0000000..1d58df8
--- /dev/null
+++ b/modules/mnv.nf
@@ -0,0 +1,92 @@
+nextflow.enable.dsl=2
+
+include { getSavePath; getSampleDir } from './utils'
+
+/* ====================================================================
+    CODON-LEVEL CHANGES (get_MNV)
+    Two SNPs in one codon change one amino acid. Annotated a base at a time they name two changes that
+    are not there: CAC>GAG is His>Glu, where each base alone says His>Asp and His>Gln, and GAG>TTG is
+    Glu>Leu where the first base alone says a stop. get_MNV groups a sample's calls by codon and reads
+    the codon whole on the reads that span it, so two SNPs on different molecules of a mixed population
+    are not merged into a change no molecule carries.
+==================================================================== */
+
+process GET_MNV {
+    tag "MNV: ${sampleId}"
+    publishDir path: { "${params.outdir}/${getSampleDir(sampleId, params)}" }, mode: params.publish_mode, saveAs: { filename -> getSavePath(filename, params) }
+    // Reads the BAM only around the called SNPs and indels: light.
+    cpus 2
+    memory { 2.GB * task.attempt }
+
+    input:
+    // The SNPs and indels that pass (CALL_FREEBAYES codon_calls), the alignment they were called on,
+    // and the reference with its GFF
+    tuple val(sampleId), val(refId), path(calls), path(calls_tbi), path(bam), path(bai), path(ref_fa), path(gff)
+
+    output:
+    tuple val(sampleId), val(refId), path("${sampleId}.${refId}.MNV.tsv"), emit: tsv
+    // BGZF-compressed but not indexed: get_MNV indexes with tabix, which its image does not carry
+    // (`tabix -p vcf` indexes the file as it is).
+    tuple val(sampleId), val(refId), path("${sampleId}.${refId}.MNV.vcf.gz"), emit: vcf
+    path("${sampleId}.${refId}.mnv.summary.json"),  emit: summary
+    path("${sampleId}.${refId}.mnv.manifest.json"), emit: manifest
+
+    script:
+    """
+    set -euo pipefail
+    # get_MNV names its outputs after its input: give the calls the sample's name.
+    ln -s ${calls} ${sampleId}.${refId}.vcf.gz
+    ln -s ${calls_tbi} ${sampleId}.${refId}.vcf.gz.tbi
+    # get_MNV reads the GFF as plain text, and the samplesheet may give it gzipped.
+    GFF=${gff}
+    if [[ "${gff}" == *.gz ]]; then
+        gzip -cd "${gff}" > get_mnv.gff
+        GFF=get_mnv.gff
+    fi
+    # The same base and mapping quality floors FreeBayes called with, so a codon's reads are the ones
+    # its calls rest on. No read or frequency floor here: the calls already passed theirs, and the
+    # run's table (MNV_TABLE) keeps a codon change only where reads carry it.
+    get_mnv --vcf ${sampleId}.${refId}.vcf.gz --bam ${bam} --fasta ${ref_fa} --gff \$GFF \\
+        --gff-features ${params.mnv_gff_features} --translation-table ${params.mnv_translation_table} \\
+        --quality ${params.freebayes_min_base_qual} --min-mapq ${params.freebayes_min_map_qual} \\
+        --threads ${task.cpus} --both --vcf-gz \\
+        --summary-json ${sampleId}.${refId}.mnv.summary.json \\
+        --run-manifest ${sampleId}.${refId}.mnv.manifest.json
+    """
+
+    stub:
+    """
+    touch ${sampleId}.${refId}.MNV.tsv
+    touch ${sampleId}.${refId}.MNV.vcf.gz
+    touch ${sampleId}.${refId}.mnv.summary.json
+    touch ${sampleId}.${refId}.mnv.manifest.json
+    """
+}
+
+process MNV_TABLE {
+    tag "MNV table"
+    publishDir "${params.outdir}", mode: params.publish_mode
+    cpus 1
+    memory { 2.GB * task.attempt }
+
+    input:
+    path(manifest)                          // sample<TAB>reference<TAB>get_MNV TSV, one line per sample
+    path(tsvs, stageAs: 'mnv/*')            // every sample's get_MNV TSV
+    val(basename)
+
+    output:
+    path("${basename}_mnv.tsv"), emit: table
+
+    script:
+    """
+    set -euo pipefail
+    awk -F'\\t' 'BEGIN { OFS = "\\t" } NF >= 3 { print \$1, \$2, "mnv/" \$3 }' ${manifest} > manifest.staged.tsv
+    python3 ${projectDir}/bin/collect_mnv.py --manifest manifest.staged.tsv \\
+        --min-reads ${params.mnv_min_reads} -o ${basename}_mnv.tsv
+    """
+
+    stub:
+    """
+    touch ${basename}_mnv.tsv
+    """
+}
diff --git a/modules/qc.nf b/modules/qc.nf
index a979c88..1e03521 100644
--- a/modules/qc.nf
+++ b/modules/qc.nf
@@ -346,6 +346,10 @@ process DUMP_VERSIONS {
       ver python    'Python '            python3   --version
       ver multiqc   'multiqc, version '  multiqc   --version
       ver pathotypr 'pathotypr '         pathotypr --version
+      # get_MNV runs in its own image until the pipeline image carries it: NA here then, and each
+      # sample's mnv/<sample>.<ref>.mnv.manifest.json records the version that ran.
+      ver get_mnv   'get_mnv '           get_mnv   --version
+      echo "get_mnv container: !{params.mnv_container ?: params.container}"
 
       # The marker catalogue is provenance too, and it is the half that cannot be recovered from
       # the results. Catalogue v1.0.0 assigned each variant a single drug inherited from its gene;
diff --git a/modules/report.nf b/modules/report.nf
index f21135a..74b4caa 100644
--- a/modules/report.nf
+++ b/modules/report.nf
@@ -135,6 +135,8 @@ process QC_REPORT {
     path(depth_profiles, stageAs: 'depth/*')   // per-sample DEPTH_PROFILE tables -> deletions + SNPs per callable kb (may be empty)
     path(snp_distances)     // SNP_DISTANCES pairs TSV -> relatedness page + GROUP_MISMATCH flag (may be NO_FILE)
     path(snp_matrix)        // SNP_MATRIX TSV -> what each series gained since its first time point (may be NO_FILE)
+    path(indel_matrix)      // INDEL_MATRIX TSV -> the indel view of the matrix panel (may be NO_FILE)
+    path(mnv)               // MNV_TABLE TSV -> the SNPs that are one codon-level change (may be NO_FILE)
     val(provenance)         // pre-quoted provenance tokens (container=..., reference=...)
     val(basename)
 
@@ -163,6 +165,8 @@ process QC_REPORT {
     DEL_ARG=""; [ -n "${depth_profiles}" ] && DEL_ARG="--depth-profiles ${depth_profiles} --out-deletions ${basename}_deletions.tsv"
     DIST_ARG=""; case "${snp_distances}" in ""|NO_FILE*) ;; *) [ -s "${snp_distances}" ] && DIST_ARG="--snp-distances ${snp_distances}";; esac
     MX_ARG="";  case "${snp_matrix}" in ""|NO_FILE*) ;; *) [ -s "${snp_matrix}" ] && MX_ARG="--snp-matrix ${snp_matrix}";; esac
+    IX_ARG="";  case "${indel_matrix}" in ""|NO_FILE*) ;; *) [ -s "${indel_matrix}" ] && IX_ARG="--indel-matrix ${indel_matrix}";; esac
+    MNV_ARG=""; case "${mnv}" in ""|NO_FILE*) ;; *) [ -s "${mnv}" ] && MNV_ARG="--mnv ${mnv}";; esac
     python3 ${projectDir}/bin/qc_report.py \\
         --summary ${summary} \\
         ${cons_arg} \\
@@ -170,7 +174,7 @@ process QC_REPORT {
         --ref-gff ${refPairs(ref_ids, gffs)} \\
         --ref-mask ${refPairs(ref_ids, masks)} \\
         --ref-fai ${refPairs(ref_ids, fais)} \\
-        \$LC_ARG \$MD_ARG \$VCF_ARG \$VH_ARG \$PL_ARG \$DR_ARG \$GC_ARG \$KRK_ARG \$DEL_ARG \$DIST_ARG \$MX_ARG \\
+        \$LC_ARG \$MD_ARG \$VCF_ARG \$VH_ARG \$PL_ARG \$DR_ARG \$GC_ARG \$KRK_ARG \$DEL_ARG \$DIST_ARG \$MX_ARG \$IX_ARG \$MNV_ARG \\
         --cluster-snps ${params.snp_cluster_threshold} \\
         --min-dp ${params.consensus_min_dp} \\
         --deletion-min-len ${params.deletion_min_len} \\
@@ -284,6 +288,7 @@ process SNP_MATRIX {
     // annotate_main_vcf both off the VCFs above ARE these files, and one task cannot stage two
     // inputs under one name.
     path(depth_vcfs, stageAs: 'depth/*')
+    path(mnv)               // the run's codon-level table (MNV_TABLE) -> codon_change columns (may be NO_FILE)
     val(reference)          // reference id the samples were mapped against
     val(basename)
 
@@ -294,10 +299,11 @@ process SNP_MATRIX {
     """
     set -euo pipefail
     DEPTH_ARG=""; [ -n "${depth_vcfs}" ] && DEPTH_ARG="--depth-vcfs ${depth_vcfs}"
+    MNV_ARG="";  case "${mnv}" in ""|NO_FILE*) ;; *) [ -s "${mnv}" ] && MNV_ARG="--mnv ${mnv}";; esac
     python3 ${projectDir}/bin/build_snp_matrix.py \\
         --vcfs ${vcfs} \\
         --reference "${reference}" \\
-        \$DEPTH_ARG \\
+        \$DEPTH_ARG \$MNV_ARG \\
         --threads ${task.cpus} \\
         -o ${basename}_snp_matrix.tsv
     """
@@ -307,3 +313,37 @@ process SNP_MATRIX {
     touch ${basename}_snp_matrix.tsv
     """
 }
+
+process INDEL_MATRIX {
+    tag "Indel matrix"
+    publishDir "${params.outdir}", mode: params.publish_mode
+    // Like SNP_MATRIX, the cost is reading every sample's all-positions VCF, in parallel.
+    cpus 4
+    memory { 4.GB * task.attempt }
+
+    input:
+    path(vcfs)                              // every sample's indel VCF (annotated when the run annotates)
+    path(depth_vcfs, stageAs: 'depth/*')    // every sample's all-positions VCF: the depth where a sample calls none
+    val(reference)
+    val(basename)
+
+    output:
+    path("${basename}_indel_matrix.tsv"), emit: matrix
+
+    script:
+    """
+    set -euo pipefail
+    DEPTH_ARG=""; [ -n "${depth_vcfs}" ] && DEPTH_ARG="--depth-vcfs ${depth_vcfs}"
+    python3 ${projectDir}/bin/build_indel_matrix.py \\
+        --vcfs ${vcfs} \\
+        --reference "${reference}" \\
+        \$DEPTH_ARG \\
+        --threads ${task.cpus} \\
+        -o ${basename}_indel_matrix.tsv
+    """
+
+    stub:
+    """
+    touch ${basename}_indel_matrix.tsv
+    """
+}
diff --git a/modules/utils.nf b/modules/utils.nf
index 416ebc3..f95050b 100644
--- a/modules/utils.nf
+++ b/modules/utils.nf
@@ -123,6 +123,7 @@ def getSavePath(filename, params) {
         lower.contains('backbone')  ||
         lower.contains('header_template') ||
         lower.contains('valid_snps_formatted') ||
+        lower.contains('codon_calls') ||
         lower == 'reference.fa' ||          
         lower.endsWith('.fai')  ||          
         lower.endsWith('.amb')  ||          
@@ -157,6 +158,12 @@ def getSavePath(filename, params) {
         return name
     }
 
+    // A4. get_MNV's codon-level tables, VCF and run records go to a folder of their own. Before the
+    //     annotation rules below, which would drop an un-annotated .vcf.gz.
+    if (lower.contains('.mnv.')) {
+        return "mnv/${name}"
+    }
+
     // B. ANNOTATION LOGIC (Filter Raw VCFs if annotation is enabled)
     if (asBool(params.annotate_legacy_vcfs)) {
         if ((lower.contains('.var.') || lower.contains('freebayes.raw')) && !lower.contains('.ann.')) {
diff --git a/modules/variants.nf b/modules/variants.nf
index 5944e83..078b455 100644
--- a/modules/variants.nf
+++ b/modules/variants.nf
@@ -30,6 +30,11 @@ process CALL_FREEBAYES {
     tuple val(sampleId), val(refId), path("${sampleId}.${refId}.var.homo.SNPs.vcf"),  emit: homo_snp
     tuple val(sampleId), val(refId), path("${sampleId}.${refId}.var.het.SNPs.vcf"),   emit: het_snp
     tuple val(sampleId), val(refId), path("${sampleId}.${refId}.var.homo.indel.vcf"), emit: homo_indel
+    tuple val(sampleId), val(refId), path("${sampleId}.${refId}.var.het.indel.vcf"),  emit: het_indel
+    // Every indel call, normalised: PASS where it meets the hom or het rule, LowSupport where not
+    tuple val(sampleId), val(refId), path("${sampleId}.${refId}.indels.vcf.gz"), path("${sampleId}.${refId}.indels.vcf.gz.tbi"), emit: indels
+    // The SNPs and indels that pass, in one VCF: what get_MNV groups by codon (never published)
+    tuple val(sampleId), val(refId), path("codon_calls.vcf.gz"), path("codon_calls.vcf.gz.tbi"), emit: codon_calls
 
     shell:
     """
@@ -101,6 +106,20 @@ process CALL_FREEBAYES {
     bcftools view \$EXCL_ARG --types snps,mnps -i "\$hom_rule" -Ov -o !{sampleId}.!{refId}.var.homo.SNPs.vcf  normalized.vcf
     bcftools view \$EXCL_ARG --types snps,mnps -i "\$het_rule" -Ov -o !{sampleId}.!{refId}.var.het.SNPs.vcf   normalized.vcf
     bcftools view \$EXCL_ARG --types indels    -i "\$hom_rule" -Ov -o !{sampleId}.!{refId}.var.homo.indel.vcf normalized.vcf
+    bcftools view \$EXCL_ARG --types indels    -i "\$het_rule" -Ov -o !{sampleId}.!{refId}.var.het.indel.vcf  normalized.vcf
+
+    # 7b. Indels, with the same definition of a call as a SNP. One file holds every indel FreeBayes
+    # made, normalised (left-aligned, one allele per record): the ones that meet the hom or het rule
+    # are PASS and the rest LowSupport, so a minority indel that falls short here keeps its fraction
+    # where the cohort matrix reads it instead of vanishing. Masked regions are left out, as for SNPs.
+    bcftools view \$EXCL_ARG --types indels -Ou normalized.vcf \
+      | bcftools filter -s LowSupport -i "\$hom_rule || \$het_rule" -Oz -o !{sampleId}.!{refId}.indels.vcf.gz
+    safe_tabix !{sampleId}.!{refId}.indels.vcf.gz
+
+    # 7c. The calls get_MNV groups by codon: the SNPs and indels that pass, in one VCF. Kept apart from
+    # the published files; the codon-level table is the deliverable.
+    bcftools view \$EXCL_ARG --types snps,indels -i "\$hom_rule || \$het_rule" -Oz -o codon_calls.vcf.gz normalized.vcf
+    safe_tabix codon_calls.vcf.gz
 
     # 8. Format for Backbone Integration
     # Adds ADP, WT, HET, HOM, NC tags to INFO/FORMAT for the consensus step. INFO is tokenised on ';'
@@ -123,6 +142,11 @@ process CALL_FREEBAYES {
     touch ${sampleId}.${refId}.var.homo.SNPs.vcf
     touch ${sampleId}.${refId}.var.het.SNPs.vcf
     touch ${sampleId}.${refId}.var.homo.indel.vcf
+    touch ${sampleId}.${refId}.var.het.indel.vcf
+    touch ${sampleId}.${refId}.indels.vcf.gz
+    touch ${sampleId}.${refId}.indels.vcf.gz.tbi
+    touch codon_calls.vcf.gz
+    touch codon_calls.vcf.gz.tbi
     """
 }
 
diff --git a/nextflow.config b/nextflow.config
index fa32e53..212f648 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -168,6 +168,22 @@ params {
     min_alt_fwd   = 2
     min_alt_rev   = 2
 
+    // --- Indels & codon-level changes ---
+    // Indels are called with the definition a SNP has (the rules above); every one of them is kept per
+    // sample in <sample>.<ref>.indels(.ann).vcf.gz, marked PASS or LowSupport, and the run's in the matrix.
+    make_indel_matrix     = true       // <samplesheet>_indel_matrix.tsv: every indel a sample calls with a PASS, each sample's AF, depth and filter
+    // Two SNPs in one codon are one amino-acid change, and annotated a base at a time they name two that
+    // are not there (CAC>GAG is His>Glu, not His>Asp and His>Gln). get_MNV reads such a codon whole, on the
+    // reads that carry both changes: per-sample tables under mnv/, the run's in <samplesheet>_mnv.tsv, and
+    // the SNP matrix's codon_change columns.
+    run_mnv               = true       // codon-level changes with get_MNV (needs the BAM, the reference and its GFF)
+    mnv_min_reads         = 2          // reads that must carry every change of a codon for it to count as one change
+    mnv_gff_features      = 'CDS'      // GFF feature types get_MNV reads codons from
+    mnv_translation_table = 11         // NCBI genetic code get_MNV translates with (11: bacteria and archaea)
+    // get_MNV runs in its bioconda image, pinned by digest, until the pipeline image carries it too; set it
+    // empty to run it in the pipeline image.
+    mnv_container         = "docker://quay.io/biocontainers/get_mnv@sha256:d443e2453b1bb831045325a563337c895d444b6135abd7b926569bf265221853"   // get_mnv 1.1.5 (bioconda)
+
     // --- Flags & Options ---
     exclude_repeats      = true
     // H37Rv Illumina "blind spots" (Zenodo 3701840): 3.62% of positions where short-read calling is
@@ -302,6 +318,8 @@ process {
 
     // Container definition (pinned digest via params.container -> reproducible builds)
     container      = params.container
+    // get_MNV has an image of its own until the pipeline image carries it (params.mnv_container)
+    withName: 'GET_MNV' { container = params.mnv_container ?: params.container }
 
     // Error handling & Shell settings
     // Capture exit codes from upstream processes when piping
@@ -415,6 +433,9 @@ profiles {
         // anything. The default keeps the scheme, because that is the path every cluster run
         // takes; this profile is the one that has to drop it.
         process.container    = params.container.replaceFirst(/^docker:\/\//, '')
+        process {
+            withName: 'GET_MNV' { container = (params.mnv_container ?: params.container).replaceFirst(/^docker:\/\//, '') }
+        }
     }
     // 'generic' = a NON-tuberculosis organism. BAMpiro's defaults are TB-tuned (diploid calling for mixed
     // infections + the MTBC-only features); this profile flips them to sane defaults for a clonal haploid
diff --git a/subworkflows/annotate.nf b/subworkflows/annotate.nf
index ff948ae..72e533e 100644
--- a/subworkflows/annotate.nf
+++ b/subworkflows/annotate.nf
@@ -6,6 +6,7 @@
 */
 
 include { ANNOTATE_LEGACY_VCF; ANNOTATE_MAIN_VCF } from '../modules/annotation'
+include { ANNOTATE_LEGACY_VCF as ANNOTATE_INDELS } from '../modules/annotation'
 include { asBool } from '../modules/utils'
 
 workflow ANNOTATE {
@@ -14,6 +15,8 @@ workflow ANNOTATE {
     homo_snp     // (sId, rId, var.homo.SNPs.vcf)
     het_snp      // (sId, rId, var.het.SNPs.vcf)
     homo_indel   // (sId, rId, var.homo.indel.vcf)
+    het_indel    // (sId, rId, var.het.indel.vcf)
+    indels       // (sId, rId, indels.vcf.gz, tbi): every indel call, PASS or LowSupport
     raw_fb       // (sId, rId, freebayes.raw.vcf.gz, tbi)
     main_vcf     // (sId, rId, vcf.gz, tbi)
     allpos_vcf   // (sId, rId, all.pos.vcf.gz, tbi)
@@ -29,6 +32,7 @@ workflow ANNOTATE {
             .mix(homo_snp.map   { sId, rId, vcf -> tuple(rId, sId, "var.homo.SNPs", vcf) })
             .mix(het_snp.map    { sId, rId, vcf -> tuple(rId, sId, "var.het.SNPs", vcf) })
             .mix(homo_indel.map { sId, rId, vcf -> tuple(rId, sId, "var.homo.indel", vcf) })
+            .mix(het_indel.map  { sId, rId, vcf -> tuple(rId, sId, "var.het.indel", vcf) })
             .mix(raw_fb.map     { sId, rId, vcf, tbi -> tuple(rId, sId, "freebayes.raw", vcf) })
 
         def ann_leg_in = leg_inputs.combine(snpeff_db, by: 0)
@@ -45,8 +49,14 @@ workflow ANNOTATE {
     // Annotate Main VCF
     def vcf_for_stats = Channel.empty()
     def ch_snpeff_stats = Channel.empty()
+    // The indels: annotated beside the main VCF (their gene, effect and HGVS), as they are otherwise
+    def indels_out = indels.map { sId, rId, vcf, tbi -> tuple(sId, rId, vcf) }
 
     if (asBool(params.annotate_main_vcf)) {
+        def ann_ind = ANNOTATE_INDELS(indels.map { sId, rId, vcf, tbi -> tuple(rId, sId, "indels", vcf) }
+                                            .combine(snpeff_db, by: 0)
+                                            .map { rId, sId, lbl, vcf, cfg, dat -> tuple(sId, rId, lbl, vcf, cfg, dat) })
+        indels_out = ann_ind.out.map { sId, rId, vcf, tbi -> tuple(sId, rId, vcf) }
         def main_pre = main_vcf.map { sId, rId, vcf, tbi -> tuple(rId, sId, vcf, tbi) }
         def ann_main_in = main_pre.combine(snpeff_db, by: 0)
             .map { rId, sId, vcf, tbi, cfg, dat -> tuple(sId, rId, vcf, tbi, cfg, dat) }
@@ -63,4 +73,5 @@ workflow ANNOTATE {
     freebayes_vcf = freebayes_ann     // annotated freebayes.raw; empty when annotate_legacy_vcfs is off
     stats_vcf     = vcf_for_stats     // (sId, rId, vcf) feeding the legacy stats and the cohort report
     snpeff_stats  = ch_snpeff_stats   // SnpEff CSV for MultiQC; empty when annotate_main_vcf is off
+    indels_vcf    = indels_out        // (sId, rId, indels(.ann).vcf.gz) feeding the indel matrix
 }
diff --git a/subworkflows/call_variants.nf b/subworkflows/call_variants.nf
index 4b6c317..91b5d80 100644
--- a/subworkflows/call_variants.nf
+++ b/subworkflows/call_variants.nf
@@ -33,6 +33,9 @@ workflow CALL_VARIANTS {
     homo_snp   = fb_out.homo_snp     // legacy split VCFs
     het_snp    = fb_out.het_snp
     homo_indel = fb_out.homo_indel
+    het_indel  = fb_out.het_indel
+    indels     = fb_out.indels       // (sId, rId, indels.vcf.gz, tbi): every indel call, PASS or LowSupport
+    codon_calls = fb_out.codon_calls // (sId, rId, codon_calls.vcf.gz, tbi): the passing SNPs + indels, for get_MNV
     allpos     = vcf_ch.allpos       // (sId, rId, all.pos.vcf.gz, tbi)
     main_vcf   = vcf_ch.main_vcf     // (sId, rId, vcf.gz, tbi)
 }
diff --git a/subworkflows/codon_changes.nf b/subworkflows/codon_changes.nf
new file mode 100644
index 0000000..c55fef6
--- /dev/null
+++ b/subworkflows/codon_changes.nf
@@ -0,0 +1,29 @@
+/*
+===============================================================================
+ Stage 6b: Codon-level changes
+===============================================================================
+ get_MNV reads the SNPs that share a codon whole, on the reads that carry them:
+ the amino-acid change a per-base annotation misses.
+*/
+
+include { GET_MNV } from '../modules/mnv'
+include { asBool } from '../modules/utils'
+
+workflow CODON_CHANGES {
+
+    take:
+    codon_calls   // (sId, rId, codon_calls.vcf.gz, tbi): the passing SNPs + indels
+    vbase         // (sId, rId, bam, bai, ref_fa, exclude_txt): the alignment they were called on
+    refGffMap     // refId -> refGff
+
+    main:
+    def mnv_tsv = Channel.empty()
+    if (asBool(params.run_mnv)) {
+        def mnv_in = codon_calls.join(vbase, by: [0, 1])
+            .map { sId, rId, vcf, tbi, bam, bai, ref_fa, excl -> tuple(sId, rId, vcf, tbi, bam, bai, ref_fa, file(refGffMap[rId])) }
+        mnv_tsv = GET_MNV(mnv_in).tsv
+    }
+
+    emit:
+    tsv = mnv_tsv // (sId, rId, <sample>.<ref>.MNV.tsv); empty when run_mnv is off
+}
diff --git a/subworkflows/cohort_report.nf b/subworkflows/cohort_report.nf
index 31c5130..0bf8457 100644
--- a/subworkflows/cohort_report.nf
+++ b/subworkflows/cohort_report.nf
@@ -6,7 +6,8 @@
 */
 
 include { ANNOTATE_CANONICAL } from '../modules/annotation'
-include { COLLECT_SUMMARY; COLLECT_DR; DEPTH_PROFILE; LIFT_VARIANTS; QC_REPORT; SNP_DISTANCES; SNP_MATRIX } from '../modules/report'
+include { COLLECT_SUMMARY; COLLECT_DR; DEPTH_PROFILE; INDEL_MATRIX; LIFT_VARIANTS; QC_REPORT; SNP_DISTANCES; SNP_MATRIX } from '../modules/report'
+include { MNV_TABLE } from '../modules/mnv'
 include { asBool } from '../modules/utils'
 
 workflow COHORT_REPORT {
@@ -20,6 +21,8 @@ workflow COHORT_REPORT {
     vcf_for_stats      // (sId, rId, vcf)
     freebayes_ann      // annotated freebayes.raw VCFs; empty when annotate_legacy_vcfs is off
     allpos_vcf         // (sId, rId, all.pos.vcf.gz, tbi): the depth of every site, for the SNP matrix
+    indels_vcf         // (sId, rId, indels(.ann).vcf.gz): every indel call, PASS or LowSupport
+    mnv_tsv            // (sId, rId, <sample>.<ref>.MNV.tsv): get_MNV's codons; empty when run_mnv is off
     gconv_cohort       // cohort gene-conversion TSV; empty when --find_gene_conversion is off
     refMap             // refId -> refFasta
     refGffMap          // refId -> refGff
@@ -36,14 +39,31 @@ workflow COHORT_REPORT {
                       .collect().ifEmpty([])
     def ref_name = refMap.keySet().join(',')
 
+    // 11a0. The run's codon-level changes: every sample's get_MNV codons that reads carry whole. The
+    // manifest names each file's sample and reference, so neither is read back out of a file name.
+    def mnv_table = asBool(params.run_mnv)
+        ? MNV_TABLE(mnv_tsv.map { sId, rId, f -> "${sId}\t${rId}\t${f.name}" }
+                           .collectFile(name: 'mnv_manifest.tsv', newLine: true),
+                    mnv_tsv.map { sId, rId, f -> f }.collect(), tsv_name).table
+              .ifEmpty(file("${projectDir}/assets/NO_FILE_MNV"))
+        : file("${projectDir}/assets/NO_FILE_MNV")
+
     // 11a. Master SNP matrix: rows = SNP sites, columns = reference/annotation + per-sample AF & depth.
     // The all-positions VCFs say what a sample's reads show at a site it has no call at. Built before
     // the report, which reads it to tell a site absent at a series' first time point from one not read.
+    // The codon-level table says which SNPs are one amino-acid change together (codon_change columns).
     def snp_matrix = asBool(params.make_snp_matrix)
         ? SNP_MATRIX(report_vcfs, allpos_vcf.map { sId, rId, vcf, tbi -> vcf }.collect().ifEmpty([]),
-                     ref_name, tsv_name).matrix
+                     mnv_table, ref_name, tsv_name).matrix
         : file("${projectDir}/assets/NO_FILE_MATRIX")
 
+    // 11a'. Master indel matrix: every indel a sample calls with a PASS, each sample's AF, depth and filter.
+    def indel_matrix = asBool(params.make_indel_matrix)
+        ? INDEL_MATRIX(indels_vcf.map { sId, rId, vcf -> vcf }.collect().ifEmpty([]),
+                       allpos_vcf.map { sId, rId, vcf, tbi -> vcf }.collect().ifEmpty([]),
+                       ref_name, tsv_name).matrix
+        : file("${projectDir}/assets/NO_FILE_INDELS")
+
     // 11b. What each sample's reads cover, from its all-positions VCF and its reference's GFF: the
     // report reads them against each other for deletions and SNPs per callable kb, and the tables
     // are published per sample whether or not a report is made.
@@ -134,7 +154,8 @@ workflow COHORT_REPORT {
         def report_gconv = gconv_cohort.ifEmpty([])
         QC_REPORT(summ.summary, summ.gene_burden, cons_files, report_gffs, report_masks, report_fais,
                   report_ref_ids, report_meta, report_vcfs, report_vcfs_h37rv, pos_liftover, dr_report,
-                  report_gconv, report_kraken, depth_profiles, snp_dist, snp_matrix, provenance, tsv_name)
+                  report_gconv, report_kraken, depth_profiles, snp_dist, snp_matrix, indel_matrix, mnv_table,
+                  provenance, tsv_name)
     }
 
 }
diff --git a/tests/pipeline/test_container_scheme_per_profile.py b/tests/pipeline/test_container_scheme_per_profile.py
index e937fae..a48915f 100644
--- a/tests/pipeline/test_container_scheme_per_profile.py
+++ b/tests/pipeline/test_container_scheme_per_profile.py
@@ -6,7 +6,8 @@
 `docker` profile is the one that strips it.
 
 Both halves are asserted here. Dropping the scheme from the default would fix a laptop and break
-every cluster, which is the more expensive direction to get wrong.
+every cluster, which is the more expensive direction to get wrong. get_MNV runs in an image of its
+own (`--mnv_container`), and the same two runtimes read that reference too.
 """
 
 from __future__ import annotations
@@ -23,13 +24,13 @@
 CONTAINER = re.compile(r"^\s*container\s*=\s*'([^']+)'", re.M)
 
 
-def effective_container(nextflow, env, profile):
-    """What `process.container` resolves to under a profile, as Nextflow itself reports it."""
+def effective_container(nextflow, env, profile, image="bampiro"):
+    """What a process's container resolves to under a profile, as Nextflow itself reports it."""
     out = subprocess.run(
         [nextflow, "config", "-profile", profile, str(REPO_ROOT)],
         capture_output=True, text=True, cwd=REPO_ROOT, env=env, timeout=300)
     assert out.returncode == 0, out.stderr[-2000:]
-    found = [m.group(1) for m in CONTAINER.finditer(out.stdout) if "bampiro" in m.group(1)]
+    found = [m.group(1) for m in CONTAINER.finditer(out.stdout) if image in m.group(1)]
     assert found, f"no container resolved under -profile {profile}"
     return found[-1]
 
@@ -57,3 +58,13 @@ def test_the_two_forms_are_the_same_image(nextflow, nextflow_env):
     laptop = effective_container(nextflow, nextflow_env, "local,docker")
 
     assert cluster == "docker://" + laptop
+
+
+def test_get_mnv_image_takes_the_scheme_each_runtime_wants(nextflow, nextflow_env):
+    """GET_MNV's container is set on its own, so the docker profile has to strip it on its own."""
+    cluster = effective_container(nextflow, nextflow_env, "garnatxa", image="get_mnv")
+    laptop = effective_container(nextflow, nextflow_env, "local,docker", image="get_mnv")
+
+    assert cluster.startswith("docker://"), cluster
+    assert "get_mnv@sha256:" in laptop, "the digest pin has to survive the rewrite"
+    assert cluster == "docker://" + laptop
diff --git a/tests/unit/test_build_indel_matrix.py b/tests/unit/test_build_indel_matrix.py
new file mode 100644
index 0000000..314bfb3
--- /dev/null
+++ b/tests/unit/test_build_indel_matrix.py
@@ -0,0 +1,132 @@
+"""Unit tests for bin/build_indel_matrix.py.
+
+The run's indels in one wide matrix: one row per indel that at least one sample calls with a PASS,
+three columns per sample (AF, DP, FT). A call below the rule keeps its fraction, marked LowSupport;
+a sample without a call reads AF 0 at the depth of the indel's anchor base where it was read.
+"""
+
+from __future__ import annotations
+
+import subprocess
+import sys
+
+from conftest import BIN, load_script
+
+load_script("build_snp_matrix")          # build_indel_matrix reads VCFs with its helpers
+bim = load_script("build_indel_matrix")
+
+HEADER = (
+    "##fileformat=VCFv4.2\n"
+    "##contig=<ID=chr1,length=1000>\n"
+    "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"g\">\n"
+    "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"d\">\n"
+    "##FORMAT=<ID=RO,Number=1,Type=Integer,Description=\"r\">\n"
+    "##FORMAT=<ID=AO,Number=A,Type=Integer,Description=\"a\">\n"
+    "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample}\n")
+ANN = ("ANN=TC|frameshift_variant|HIGH|Rv0678|gene-Rv0678|transcript|t|protein_coding|1/1|"
+       "c.144dupC|p.Glu49fs|145/498|145/498|49/165||")
+
+
+def _vcf(path, sample, rows):
+    """rows: (pos, ref, alt, filter, ro, ao, info)."""
+    body = "".join(f"chr1\t{pos}\t.\t{ref}\t{alt}\t50\t{flt}\t{info or '.'}\tGT:DP:RO:AO\t0/1:{ro + ao}:{ro}:{ao}\n"
+                   for pos, ref, alt, flt, ro, ao, info in rows)
+    path.write_text(HEADER.format(sample=sample) + body)
+    return path
+
+
+def _allpos(path, sample, depth):
+    """An all-positions VCF with the given depth at every position 1..300."""
+    body = "".join(f"chr1\t{p}\t.\tA\t.\t.\t.\t.\tGT:DP:AD\t0:{depth}:{depth},0\n" for p in range(1, 301))
+    path.write_text(HEADER.format(sample=sample).replace(
+        "##FORMAT=<ID=RO", "##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"a\">\n##FORMAT=<ID=RO") + body)
+    return path
+
+
+def _run(tmp_path, vcfs, depth_vcfs=()):
+    out = tmp_path / "indel_matrix.tsv"
+    args = [sys.executable, str(BIN / "build_indel_matrix.py"), "--vcfs", *map(str, vcfs), "--reference", "ref",
+            "-o", str(out)]
+    if depth_vcfs:
+        args += ["--depth-vcfs", *map(str, depth_vcfs)]
+    proc = subprocess.run(args, capture_output=True, text=True)
+    assert proc.returncode == 0, proc.stderr
+    lines = out.read_text().splitlines()
+    header = lines[0].split("\t")
+    return header, [dict(zip(header, ln.split("\t"))) for ln in lines[1:]], proc.stderr
+
+
+def test_a_row_needs_a_pass_in_at_least_one_sample(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN),
+                                       (200, "CG", "C", "LowSupport", 38, 2, None)])
+    header, rows, err = _run(tmp_path, [a])
+    assert [r["pos"] for r in rows] == ["100"], "an indel no sample calls with a PASS is left out"
+    assert "1 more only below the rule" in err
+
+
+def test_a_call_below_the_rule_keeps_its_fraction_marked_low_support(tmp_path):
+    """A minority frameshift starts below the rule: blanking it would say the sample does not have it."""
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN)])
+    b = _vcf(tmp_path / "b.vcf", "B", [(100, "T", "TC", "LowSupport", 36, 4, ANN)])
+    _, rows, _ = _run(tmp_path, [a, b])
+    assert (rows[0]["B|AF"], rows[0]["B|DP"], rows[0]["B|FT"]) == ("0.1000", "40", "LowSupport")
+    assert (rows[0]["A|AF"], rows[0]["A|FT"]) == ("1.0000", "PASS")
+    assert rows[0]["n_pass"] == "1"
+
+
+def test_a_sample_without_the_indel_reads_zero_at_the_anchor_depth(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN)])
+    b = _vcf(tmp_path / "b.vcf", "B", [])
+    _, rows, _ = _run(tmp_path, [a, b], [_allpos(tmp_path / "a.all.vcf", "A", 40), _allpos(tmp_path / "b.all.vcf", "B", 55)])
+    assert (rows[0]["B|AF"], rows[0]["B|DP"], rows[0]["B|FT"]) == ("0", "55", "")
+
+
+def test_a_sample_nobody_read_there_is_not_zero(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN)])
+    b = _vcf(tmp_path / "b.vcf", "B", [])
+    _, rows, _ = _run(tmp_path, [a, b], [_allpos(tmp_path / "b.all.vcf", "B", 0)])
+    assert (rows[0]["B|AF"], rows[0]["B|DP"]) == ("", "0"), "absence where nothing was read says nothing"
+
+
+def test_without_depths_an_absent_cell_stays_blank(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN)])
+    b = _vcf(tmp_path / "b.vcf", "B", [])
+    _, rows, _ = _run(tmp_path, [a, b])
+    assert (rows[0]["B|AF"], rows[0]["B|DP"], rows[0]["B|FT"]) == ("", "", "")
+
+
+def test_snps_are_left_to_the_snp_matrix(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(50, "A", "G", "PASS", 0, 40, None), (60, "AC", "GT", "PASS", 0, 40, None),
+                                       (100, "T", "TC", "PASS", 0, 40, None)])
+    _, rows, _ = _run(tmp_path, [a])
+    assert [r["pos"] for r in rows] == ["100"]
+
+
+def test_the_annotation_and_the_length_are_carried(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, ANN), (200, "CGAT", "C", "PASS", 0, 40, None)])
+    header, rows, _ = _run(tmp_path, [a])
+    assert header[:11] == ["reference", "contig", "pos", "ref_allele", "alt_allele", "length", "gene", "effect",
+                           "hgvs_c", "hgvs_p", "n_pass"]
+    ins, dele = rows
+    assert (ins["length"], ins["gene"], ins["effect"], ins["hgvs_c"], ins["hgvs_p"]) == \
+           ("1", "Rv0678", "frameshift_variant", "c.144dupC", "p.Glu49fs")
+    assert dele["length"] == "-3"
+
+
+def test_an_unfiltered_vcf_reads_as_passing(tmp_path):
+    """An indel VCF written without the soft filter (FILTER '.') is every call the caller kept."""
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", ".", 0, 40, None)])
+    _, rows, _ = _run(tmp_path, [a])
+    assert rows[0]["A|FT"] == "PASS" and rows[0]["n_pass"] == "1"
+
+
+def test_two_indels_at_one_position_are_two_rows(tmp_path):
+    a = _vcf(tmp_path / "a.vcf", "A", [(100, "T", "TC", "PASS", 0, 40, None)])
+    b = _vcf(tmp_path / "b.vcf", "B", [(100, "T", "TCC", "PASS", 0, 40, None)])
+    _, rows, _ = _run(tmp_path, [a, b])
+    assert [(r["alt_allele"], r["A|AF"], r["B|AF"]) for r in rows] == [("TC", "1.0000", ""), ("TCC", "", "1.0000")]
+
+
+def test_ann_reads_gene_effect_and_both_hgvs():
+    assert bim._ann(ANN) == ("Rv0678", "frameshift_variant", "c.144dupC", "p.Glu49fs")
+    assert bim._ann("DP=10") == ("", "", "", "")
diff --git a/tests/unit/test_build_snp_matrix.py b/tests/unit/test_build_snp_matrix.py
index 9e4f577..cab7130 100644
--- a/tests/unit/test_build_snp_matrix.py
+++ b/tests/unit/test_build_snp_matrix.py
@@ -579,3 +579,58 @@ def test_end_to_end_an_unreadable_depth_file_leaves_that_sample_unknown(tmp_path
     assert "could not read depths" in err
     assert rows["chr1:100"][10:12] == ["", ""]
     assert rows["chr1:200"][8:10] == ["0", "38"], "the other samples are unaffected"
+
+
+# --------------------------------------------------------------------------
+# --mnv: the SNPs that are one codon-level change with another one on the same reads
+# --------------------------------------------------------------------------
+
+def _mnv_table(path, rows):
+    header = ["sample", "reference", "contig", "gene", "positions", "aa_change", "mnv_frequency"]
+    path.write_text("\t".join(header) + "\n" + "".join("\t".join(r) + "\n" for r in rows))
+    return path
+
+
+def _matrix_with(tmp_path, extra):
+    vcf = tmp_path / "S1.vcf"
+    vcf.write_text("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n"
+                   "chr1\t139\t.\tC\tG\t50\t.\t.\tGT:AD:DP\t1/1:0,40:40\n"
+                   "chr1\t141\t.\tC\tG\t50\t.\t.\tGT:AD:DP\t1/1:0,40:40\n"
+                   "chr1\t155\t.\tC\tT\t50\t.\t.\tGT:AD:DP\t1/1:0,40:40\n")
+    out = tmp_path / "m.tsv"
+    proc = subprocess.run([sys.executable, str(bsm.__file__), "--vcfs", str(vcf), "-o", str(out), *extra],
+                          capture_output=True, text=True)
+    assert proc.returncode == 0, proc.stderr
+    lines = out.read_text().splitlines()
+    return lines[0].split("\t"), [ln.split("\t") for ln in lines[1:]], proc.stderr
+
+
+def test_without_the_codon_table_the_layout_is_the_old_one(tmp_path):
+    header, _, _ = _matrix_with(tmp_path, [])
+    assert header == ["reference", "contig", "pos", "ref_allele", "alt_allele", "gene", "effect", "aa_change",
+                      "S1|AF", "S1|DP"]
+
+
+def test_a_snp_of_a_phased_codon_carries_the_codon_change(tmp_path):
+    table = _mnv_table(tmp_path / "mnv.tsv", [["S1", "ref", "chr1", "rpoB", "139,141", "His445Glu", "0.9904"]])
+    header, rows, err = _matrix_with(tmp_path, ["--mnv", str(table)])
+    assert header[8:10] == ["codon_change", "codon_change_samples"]
+    by = {r[2]: r for r in rows}
+    assert by["139"][8:10] == ["His445Glu", "His445Glu=S1:0.9904"]
+    assert by["141"][8:10] == ["His445Glu", "His445Glu=S1:0.9904"]
+    assert by["155"][8:10] == ["", ""], "a SNP that stands alone"
+    assert by["155"][10:] == ["1.0000", "40"], "the sample columns follow"
+    assert "2 of them part of a codon-level change" in err
+
+
+def test_a_placeholder_table_leaves_the_layout_alone(tmp_path):
+    ph = tmp_path / "NO_FILE_MNV"
+    ph.write_text("")
+    header, _, _ = _matrix_with(tmp_path, ["--mnv", str(ph)])
+    assert "codon_change" not in header
+
+
+def test_codon_cells_list_every_change_and_its_samples():
+    assert bsm._codon_cells(None) == ("", "")
+    assert bsm._codon_cells({"His445Glu": [("A", "0.99"), ("B", "")], "His445Lys": [("C", "0.5")]}) == \
+           ("His445Glu;His445Lys", "His445Glu=A:0.99,B;His445Lys=C:0.5")
diff --git a/tests/unit/test_collect_mnv.py b/tests/unit/test_collect_mnv.py
new file mode 100644
index 0000000..e5bb0f4
--- /dev/null
+++ b/tests/unit/test_collect_mnv.py
@@ -0,0 +1,106 @@
+"""Unit tests for bin/collect_mnv.py: the run's codon-level changes from get_MNV.
+
+get_MNV writes a combined amino acid for every codon holding more than one SNV, also when no read
+carries them together (two SNVs on different molecules). A codon change is kept only where reads
+carry every change of the codon at once.
+"""
+
+from __future__ import annotations
+
+import csv
+import subprocess
+import sys
+
+from conftest import BIN, load_script
+
+cm = load_script("collect_mnv")
+
+COLS = ["Chromosome", "Gene", "Positions", "Reference Bases", "Base Changes", "AA Changes", "SNP AA Changes",
+        "Variant Type", "Change Type", "Reference Codon", "SNP Codon", "MNV Codon", "SNP Reads", "MNV Reads",
+        "Total Reads", "SNP Frequencies", "MNV Frequencies", "MNV Consequence Shift", "MNV Phasing Support"]
+# Rows as get_MNV 1.1.5 writes them for a simulated sample (tests/unit data from a real run of the tool).
+ROWS = [
+    ["win", "rpoB_Rv0667", "6155", "C", "T", "Ser450Leu", "Ser450Leu", "SNP", "Non-synonymous", "TCG", "TTG", "TTG",
+     "102", "0", "103", "0.9903", "0.0000", "-", "-"],
+    # on the same reads: His445Glu, where each base alone says His445Asp and His445Gln
+    ["win", "rpoB_Rv0667", "6139, 6141", "C, C", "G, G", "His445Glu", "His445Asp, His445Gln", "SNP/MNV",
+     "Non-synonymous", "CAC", "GAC, CAG", "GAG", "2, 1", "103", "104", "0.0189, 0.0096", "0.9904", "Concordant",
+     "1.0000"],
+    # two SNVs of one codon on different molecules: get_MNV still writes a combined change, carried by no read
+    ["win", "rpoB_Rv0667", "6094, 6096", "C, G", "A, C", "Leu430Ile", "Leu430Met, Leu430Leu", "SNP/MNV",
+     "Non-synonymous", "CTG", "ATG, CTC", "ATC", "89, 26", "0", "115", "0.7739, 0.2261", "0.0000", "Concordant",
+     "0.0000"],
+    # a minority MNV whose bases alone would name a stop
+    ["win", "rpoB_Rv0667", "6184, 6185", "G, A", "T, T", "Glu460Leu", "Glu460*, Glu460Val", "SNP/MNV",
+     "Non-synonymous", "GAG", "TAG, GTG", "TTG", "0, 0", "23", "99", "0.0000, 0.0000", "0.2323", "MNV-masked",
+     "1.0000"],
+    ["win", "NP_215192.1_Rv0678", "24130", "T", "TC", "Glu49Argfs", "-", "INDEL", "Frameshift Indel", "GAT", "-",
+     "GAT", "-", "-", "-", "-", "-", "-", "-"],
+]
+
+
+def _tsv(path, rows=ROWS):
+    with open(path, "w", newline="") as fh:
+        w = csv.writer(fh, delimiter="\t", lineterminator="\n")
+        w.writerow(COLS)
+        w.writerows(rows)
+    return path
+
+
+def _run(tmp_path, manifest_lines, min_reads=2):
+    manifest = tmp_path / "manifest.tsv"
+    manifest.write_text("".join("\t".join(x) + "\n" for x in manifest_lines))
+    out = tmp_path / "mnv.tsv"
+    proc = subprocess.run([sys.executable, str(BIN / "collect_mnv.py"), "--manifest", str(manifest),
+                           "--min-reads", str(min_reads), "-o", str(out)], capture_output=True, text=True)
+    assert proc.returncode == 0, proc.stderr
+    with open(out, newline="") as fh:
+        return list(csv.DictReader(fh, delimiter="\t")), proc.stderr
+
+
+def test_only_codons_that_reads_carry_whole_are_kept(tmp_path):
+    rows, err = _run(tmp_path, [("S1", "win", str(_tsv(tmp_path / "S1.win.MNV.tsv")))])
+    assert [r["aa_change"] for r in rows] == ["His445Glu", "Glu460Leu"]
+    assert "Leu430Ile" not in {r["aa_change"] for r in rows}, "no read carries it: a change nobody has"
+    assert "1 of them masked" in err
+
+
+def test_a_row_says_what_the_bases_say_alone_and_how_many_reads_carry_it(tmp_path):
+    rows, _ = _run(tmp_path, [("S1", "win", str(_tsv(tmp_path / "S1.win.MNV.tsv")))])
+    his = rows[0]
+    assert (his["sample"], his["reference"], his["gene"], his["positions"]) == ("S1", "win", "rpoB_Rv0667", "6139,6141")
+    assert (his["ref_codon"], his["mnv_codon"], his["snp_aa_changes"]) == ("CAC", "GAG", "His445Asp, His445Gln")
+    assert (his["mnv_reads"], his["total_reads"], his["mnv_frequency"], his["phasing_support"]) == \
+           ("103", "104", "0.9904", "1.0000")
+    assert rows[1]["consequence_shift"] == "MNV-masked"
+
+
+def test_the_read_floor_is_the_one_given(tmp_path):
+    rows, _ = _run(tmp_path, [("S1", "win", str(_tsv(tmp_path / "S1.win.MNV.tsv")))], min_reads=50)
+    assert [r["aa_change"] for r in rows] == ["His445Glu"]
+
+
+def test_samples_and_references_come_from_the_manifest_not_the_file_name(tmp_path):
+    a, b = _tsv(tmp_path / "x.tsv"), _tsv(tmp_path / "y.tsv", ROWS[1:2])
+    rows, _ = _run(tmp_path, [("A.1", "refA", str(a)), ("B", "refB", str(b))])
+    assert sorted({(r["sample"], r["reference"]) for r in rows}) == [("A.1", "refA"), ("B", "refB")]
+
+
+def test_an_empty_or_missing_file_is_skipped_and_counted(tmp_path):
+    empty = tmp_path / "empty.tsv"
+    empty.write_text("")
+    rows, err = _run(tmp_path, [("S1", "win", str(empty)), ("S2", "win", str(tmp_path / "nope.tsv"))])
+    assert rows == []
+    assert "2 empty or missing file(s)" in err
+
+
+def test_rows_are_ordered_by_position_then_sample(tmp_path):
+    rows, _ = _run(tmp_path, [("S2", "win", str(_tsv(tmp_path / "a.tsv"))), ("S1", "win", str(_tsv(tmp_path / "b.tsv")))])
+    assert [(r["positions"], r["sample"]) for r in rows] == [("6139,6141", "S1"), ("6139,6141", "S2"),
+                                                             ("6184,6185", "S1"), ("6184,6185", "S2")]
+
+
+def test_helpers():
+    assert cm._list("6139, 6141") == "6139,6141"
+    assert cm._clean("-") == "" and cm._clean(".") == "" and cm._clean(" x ") == "x"
+    assert cm._int("12") == 12 and cm._int("-") is None
diff --git a/tests/unit/test_get_mnv_script.py b/tests/unit/test_get_mnv_script.py
new file mode 100644
index 0000000..183f52d
--- /dev/null
+++ b/tests/unit/test_get_mnv_script.py
@@ -0,0 +1,75 @@
+"""The GET_MNV script hands get_MNV a GFF it can read.
+
+get_MNV reads the GFF as plain text and fails on a gzipped one ("stream did not contain valid
+UTF-8"), while the samplesheet takes either: the SnpEff build and every other GFF reader of the
+pipeline open both. A stub run skips the script block, so the script is rendered here and run in
+bash, with a get_mnv on PATH that records what it was given.
+"""
+
+from __future__ import annotations
+
+import gzip
+import os
+import re
+import stat
+import subprocess
+
+from conftest import REPO_ROOT
+
+MODULE = REPO_ROOT / "modules" / "mnv.nf"
+GFF = "##gff-version 3\nc1\tRefSeq\tCDS\t1\t300\t.\t+\t0\tID=cds-1;gene=rpoB\n"
+
+
+def get_mnv_script():
+    text = MODULE.read_text()
+    body = text[text.index("process GET_MNV"):]
+    m = re.search(r'script:\s*"""(.*?)"""', body, re.S)
+    assert m, "GET_MNV has no script block"
+    return m.group(1)
+
+
+def render(block, values):
+    """What Nextflow runs: ${...} interpolated, \\$ and \\\\ unescaped."""
+    def sub(m):
+        t = m.group(0)
+        if t == "\\\\":
+            return "\\"
+        if t == "\\$":
+            return "$"
+        return str(values[t[2:-1]])
+    return re.sub(r"\\\\|\\\$|\$\{[^}]+\}", sub, block)
+
+
+def run_with(tmp_path, gff_name):
+    for f in ("codon_calls.vcf.gz", "codon_calls.vcf.gz.tbi", "S1.bam", "S1.bam.bai", "reference.fa"):
+        (tmp_path / f).write_text("")
+    bindir = tmp_path / "bin"
+    bindir.mkdir()
+    fake = bindir / "get_mnv"
+    fake.write_text('#!/usr/bin/env bash\nprintf \'%s\\n\' "$@" > get_mnv.args\n')
+    fake.chmod(fake.stat().st_mode | stat.S_IEXEC)
+    values = {"calls": "codon_calls.vcf.gz", "calls_tbi": "codon_calls.vcf.gz.tbi", "sampleId": "S1", "refId": "r",
+              "bam": "S1.bam", "ref_fa": "reference.fa", "gff": gff_name, "task.cpus": 2,
+              "params.mnv_gff_features": "CDS", "params.mnv_translation_table": 11,
+              "params.freebayes_min_base_qual": 20, "params.freebayes_min_map_qual": 30}
+    env = dict(os.environ, PATH=f"{bindir}{os.pathsep}{os.environ['PATH']}")
+    proc = subprocess.run(["bash", "-c", render(get_mnv_script(), values)], cwd=tmp_path, env=env,
+                          capture_output=True, text=True)
+    assert proc.returncode == 0, proc.stderr
+    args = (tmp_path / "get_mnv.args").read_text().splitlines()
+    return args, tmp_path / args[args.index("--gff") + 1]
+
+
+def test_a_gzipped_gff_reaches_get_mnv_as_plain_text(tmp_path):
+    with gzip.open(tmp_path / "genes.gff.gz", "wt") as fh:
+        fh.write(GFF)
+    args, gff = run_with(tmp_path, "genes.gff.gz")
+    assert gff.name != "genes.gff.gz", "get_MNV was handed the gzipped GFF, which it cannot read"
+    assert gff.read_text() == GFF
+
+
+def test_a_plain_gff_is_handed_over_as_given(tmp_path):
+    (tmp_path / "genes.gff").write_text(GFF)
+    args, gff = run_with(tmp_path, "genes.gff")
+    assert gff.name == "genes.gff"
+    assert args[args.index("--vcf") + 1] == "S1.r.vcf.gz", "get_MNV names its outputs after the calls it reads"
diff --git a/tests/unit/test_qc_report_indels_mnv.py b/tests/unit/test_qc_report_indels_mnv.py
new file mode 100644
index 0000000..c46ec23
--- /dev/null
+++ b/tests/unit/test_qc_report_indels_mnv.py
@@ -0,0 +1,91 @@
+"""The report's indel view and its codon-level marks.
+
+parse_indel_matrix turns the run's indel matrix into the matrix panel's second view, in the SNP view's
+column order; parse_mnv_table and build_snp_matrix mark the SNPs that are one amino-acid change with
+another one of their codon on the same reads.
+"""
+
+from __future__ import annotations
+
+from qcreport import panels as qc_panels
+from qcreport import parsers as qc_parsers
+
+INDEL_HEADER = ["reference", "contig", "pos", "ref_allele", "alt_allele", "length", "gene", "effect", "hgvs_c",
+                "hgvs_p", "n_pass", "A|AF", "A|DP", "A|FT", "B|AF", "B|DP", "B|FT"]
+
+
+def write(path, lines):
+    path.write_text("".join("\t".join(map(str, ln)) + "\n" for ln in lines), encoding="utf-8")
+    return str(path)
+
+
+def test_the_indel_view_keeps_only_called_cells_and_their_filter(tmp_path):
+    p = write(tmp_path / "ind.tsv", [INDEL_HEADER,
+              ["r", "c1", 100, "T", "TC", 1, "Rv0678", "frameshift_variant", "c.144dupC", "p.Glu49fs", 1,
+               "1.0000", 40, "PASS", "0.1000", 40, "LowSupport"],
+              ["r", "c1", 200, "CGA", "C", -2, "g", "frameshift_variant", "", "", 1, "0.9000", 30, "PASS", "0", 50, ""]])
+    m = qc_parsers.parse_indel_matrix(p)
+    assert m["samples"] == ["A", "B"]
+    first, second = m["rows"]
+    assert (first["gene"], first["len"], first["aa"], first["hgvs_c"]) == ("Rv0678", 1, "p.Glu49fs", "c.144dupC")
+    assert first["cells"] == {0: [1.0, 40, "PASS"], 1: [0.1, 40, "LowSupport"]}
+    assert second["cells"] == {0: [0.9, 30, "PASS"]}, "AF 0 is read without the indel, not a call"
+    assert first["n"] == 2 and second["n"] == 1
+
+
+def test_the_indel_view_follows_the_snp_view_column_order(tmp_path):
+    p = write(tmp_path / "ind.tsv", [INDEL_HEADER,
+              ["r", "c1", 100, "T", "TC", 1, "", "", "", "", 1, "1.0000", 40, "PASS", "1.0000", 40, "PASS"]])
+    m = qc_parsers.parse_indel_matrix(p, order=["B", "Z", "A"])
+    assert m["samples"] == ["B", "A"]
+    assert m["rows"][0]["cells"] == {0: [1.0, 40, "PASS"], 1: [1.0, 40, "PASS"]}
+
+
+def test_an_indel_matrix_without_calls_or_a_placeholder_hides_the_view(tmp_path):
+    empty = write(tmp_path / "ind.tsv", [INDEL_HEADER])
+    assert qc_parsers.parse_indel_matrix(empty) is None
+    ph = tmp_path / "NO_FILE_INDELS"
+    ph.write_text("")
+    assert qc_parsers.parse_indel_matrix(str(ph)) is None
+    assert qc_parsers.parse_indel_matrix(None) is None
+
+
+def test_the_codon_table_maps_every_position_of_a_codon(tmp_path):
+    p = write(tmp_path / "mnv.tsv", [["sample", "reference", "contig", "gene", "positions", "aa_change",
+                                      "snp_aa_changes", "consequence_shift", "mnv_frequency"],
+                                     ["S1", "r", "c1", "rpoB", "139,141", "His445Glu", "His445Asp, His445Gln",
+                                      "Concordant", "0.9904"]])
+    m = qc_parsers.parse_mnv_table(p)
+    entry = ["His445Glu", 0.9904, "His445Asp, His445Gln", "Concordant"]
+    assert m == {"S1": {"c1:139": entry, "c1:141": entry}}
+    assert qc_parsers.parse_mnv_table(None) == {}
+
+
+def test_the_snp_view_marks_only_the_sample_whose_reads_carry_the_codon_change():
+    variants = {"S1": {"c1:139": {"ref": "C", "alt": "G", "af": 1.0, "dp": 40, "aa": "p.His445Asp"}},
+                "S2": {"c1:139": {"ref": "C", "alt": "G", "af": 1.0, "dp": 38, "aa": "p.His445Asp"}}}
+    mnv = {"S1": {"c1:139": ["His445Glu", 0.99, "His445Asp, His445Gln", "Concordant"]}}
+    row = qc_panels.build_snp_matrix(variants, mnv=mnv)["rows"][0]
+    assert row["mnv"] == {0: ["His445Glu", 0.99, "His445Asp, His445Gln", "Concordant"]}
+    assert row["mnv_aa"] == ["His445Glu"]
+
+
+def test_a_snp_view_without_codon_changes_carries_no_marks():
+    variants = {"S1": {"c1:139": {"ref": "C", "alt": "G", "af": 1.0, "dp": 40}}}
+    row = qc_panels.build_snp_matrix(variants)["rows"][0]
+    assert "mnv" not in row and "mnv_aa" not in row
+
+
+def test_an_mnp_record_is_every_snp_it_changes(tmp_path):
+    """FreeBayes writes two changes of one codon on the same reads as one record (CAC>GAG): each base
+    it changes is a SNP at its own position, with the record's fraction and its codon-level annotation."""
+    vcf = tmp_path / "s.vcf"
+    vcf.write_text("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n"
+                   "c1\t139\t.\tCAC\tGAG\t50\t.\tANN=GAG|missense_variant|MODERATE|rpoB|g|transcript|t|protein_coding|1/1|c.1333_1335delCACinsGAG|"
+                   "p.His445Glu|||||\tGT:AD:DP\t1/1:1,103:104\n"
+                   "c1\t200\t.\tAT\tA\t50\t.\t.\tGT:AD:DP\t1/1:0,40:40\n")
+    v = qc_parsers.parse_vcfs([str(vcf)])["S1"]
+    assert sorted(v) == ["c1:139", "c1:141"], "the unchanged middle base is no SNP, and the indel is not one"
+    assert (v["c1:139"]["ref"], v["c1:139"]["alt"], v["c1:141"]["ref"], v["c1:141"]["alt"]) == ("C", "G", "C", "G")
+    assert v["c1:139"]["aa"] == v["c1:141"]["aa"] == "p.His445Glu"
+    assert v["c1:141"]["af"] == round(103 / 104, 4)
diff --git a/tests/unit/test_vcf_filter_rules.py b/tests/unit/test_vcf_filter_rules.py
index 827306c..8c61ebb 100644
--- a/tests/unit/test_vcf_filter_rules.py
+++ b/tests/unit/test_vcf_filter_rules.py
@@ -125,6 +125,31 @@ def test_the_boundary_fraction_is_homozygous_and_not_heterozygous(self, tmp_path
         assert self._select(vcf, het()) == []
 
 
+@pytest.mark.skipif(not shutil.which("bcftools"), reason="bcftools is not installed")
+class TestIndelSoftFilter:
+    """CALL_FREEBAYES keeps every indel and marks the ones short of the rule instead of dropping them:
+    `bcftools view --types indels | bcftools filter -s LowSupport -i "<hom> || <het>"`, as in the process."""
+
+    def test_every_indel_is_kept_and_only_the_short_ones_are_marked(self, tmp_path):
+        header = TestAgainstRealBcftools._vcf(tmp_path, []).read_text()
+        rows = [
+            "c\t100\t.\tA\tAT\t50\t.\tDP=50;RO=1;AO=49;SAF=25;SAR=24\tGT:DP\t1/1:50\n",    # hom indel
+            "c\t200\t.\tAT\tA\t50\t.\tDP=50;RO=40;AO=10;SAF=5;SAR=5\tGT:DP\t0/1:50\n",    # 20% indel
+            "c\t300\t.\tAT\tA\t50\t.\tDP=50;RO=46;AO=4;SAF=4;SAR=0\tGT:DP\t0/1:50\n",     # one strand
+            "c\t400\t.\tA\tT\t50\t.\tDP=50;RO=1;AO=49;SAF=25;SAR=24\tGT:DP\t1/1:50\n",     # a SNP
+        ]
+        vcf = tmp_path / "calls.vcf"
+        vcf.write_text(header + "".join(rows))
+        rule = f"{hom()} || {het()}"
+        view = subprocess.run(["bcftools", "view", "--types", "indels", "-Ou", str(vcf)],
+                              capture_output=True, check=True)
+        out = subprocess.run(["bcftools", "filter", "-s", "LowSupport", "-i", rule, "-Ov", "-"],
+                             input=view.stdout, capture_output=True, check=True)
+        got = [(ln.split(b"\t")[1].decode(), ln.split(b"\t")[6].decode())
+               for ln in out.stdout.splitlines() if not ln.startswith(b"#")]
+        assert got == [("100", "PASS"), ("200", "PASS"), ("300", "LowSupport")]
+
+
 class TestCli:
     """The processes call this as a command and substitute the result into bcftools."""