From 83e402fbed3da1ae11d8d8b71bb0007eead63665 Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Thu, 24 Sep 2026 19:51:39 +0200 Subject: [PATCH 1/2] Deliver indels, and read two SNPs of one codon as one change The main VCF holds SNPs only, and of the indels only the homozygous ones reached a filtered file, the legacy var.homo.indel.vcf. A minority indel was in none, and a minority frameshift is how resistance through Rv0678 often starts. CALL_FREEBAYES now writes each sample's indels, normalised, with FILTER PASS where they meet the hom or het rule a SNP has and LowSupport where not; annotated, they are ..indels.ann.vcf.gz. INDEL_MATRIX puts the run's indels in _indel_matrix.tsv: every indel a sample calls with a PASS, its length, gene, effect and HGVS, and per sample the fraction, the depth and the filter, with AF 0 at the anchor's depth where a sample was read without it. The report's matrix panel gains an Indels view, and the legacy split set gains var.het.indel.vcf. Two SNPs in one codon change one amino acid, and annotated a base at a time they name two changes that are not there. GET_MNV runs get_MNV 1.1.5 on each sample's passing SNPs and indels with the BAM they were called on, and reads every codon whole on the reads that span it. MNV_TABLE keeps in _mnv.tsv the codons at least --mnv_min_reads reads carry whole: get_MNV also writes a combined change for two SNPs of one codon on different molecules, which no read carries. With it, the SNP matrix gains codon_change and codon_change_samples. On reads simulated from H37Rv, rpoB codon 445 reads His445Glu on 103 of 104 reads (base by base, His445Asp and His445Gln), a pair in codon 460 on 23% of the reads reads Glu460Leu where its first base alone says a stop, and two SNPs of codon 430 on different molecules stay two changes. get_MNV runs in its bioconda image, pinned by digest, until the pipeline image is rebuilt with it; the Dockerfile now installs it. The report read single-base records only, so every SNP FreeBayes wrote as an MNP (CAC>GAG: two changes of one codon on the same reads) was in the SNP matrix and in none of the report's tables. It now reads such a record as the SNPs it is made of. --- .github/dockerfile/Dockerfile | 4 + CHANGELOG.md | 37 ++++ README.md | 2 + assets/NO_FILE.README.md | 2 + assets/NO_FILE_INDELS | 0 assets/NO_FILE_MNV | 0 bin/build_indel_matrix.py | 167 ++++++++++++++++++ bin/build_snp_matrix.py | 48 ++++- bin/collect_mnv.py | 109 ++++++++++++ bin/qc_report.py | 13 +- bin/qcreport/panels.py | 15 +- bin/qcreport/parsers.py | 121 ++++++++++++- bin/qcreport/render.py | 7 +- bin/report_assets/css/04_panels.css | 10 ++ bin/report_assets/js/12_snpmatrix.js | 60 +++++-- bin/report_assets/shell.html | 4 +- docs/configuration.md | 8 + docs/outputs.md | 61 ++++++- docs/qc-report.md | 2 +- main.nf | 10 +- modules/mnv.nf | 91 ++++++++++ modules/qc.nf | 4 + modules/report.nf | 44 ++++- modules/utils.nf | 7 + modules/variants.nf | 24 +++ nextflow.config | 21 +++ subworkflows/annotate.nf | 11 ++ subworkflows/call_variants.nf | 3 + subworkflows/codon_changes.nf | 29 +++ subworkflows/cohort_report.nf | 27 ++- .../test_container_scheme_per_profile.py | 19 +- tests/unit/test_build_indel_matrix.py | 132 ++++++++++++++ tests/unit/test_build_snp_matrix.py | 55 ++++++ tests/unit/test_collect_mnv.py | 106 +++++++++++ tests/unit/test_get_mnv_script.py | 75 ++++++++ tests/unit/test_qc_report_indels_mnv.py | 91 ++++++++++ tests/unit/test_vcf_filter_rules.py | 25 +++ 37 files changed, 1403 insertions(+), 41 deletions(-) create mode 100644 assets/NO_FILE_INDELS create mode 100644 assets/NO_FILE_MNV create mode 100644 bin/build_indel_matrix.py create mode 100644 bin/collect_mnv.py create mode 100644 modules/mnv.nf create mode 100644 subworkflows/codon_changes.nf create mode 100644 tests/unit/test_build_indel_matrix.py create mode 100644 tests/unit/test_collect_mnv.py create mode 100644 tests/unit/test_get_mnv_script.py create mode 100644 tests/unit/test_qc_report_indels_mnv.py 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 + `..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. `_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 `_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.** `..ann.vcf.gz` holds SNPs, one + record per base; the indels are in `..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 _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))+'…'+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 ''; }).join('')+ ''):''; + // SNPs or indels: each view is its own matrix, in the same sample columns + var modeUI=(hasS&&hasI)?('
'+ + ''+ + '
'):''; host.innerHTML= - '
'+ - ''+ + '
'+modeUI+ + ''+ ''+ ''+ ''+ @@ -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) × '+vi.length+' sample(s)'+(nfilt?' (filtered)':'')+ + el('snpmxcount').innerHTML=total+' '+unit+' × '+vi.length+' sample(s)'+(nfilt?' (filtered)':'')+ (capped?(' · showing '+(snpmxAll?('all '+total):(MAXR+' of '+total))+''):'')+(M.truncated?' · 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 ''+esc(v)+'';}).join('')+''; }).join(''):''; var stop=nf*mh; - var refhdr='SNP'+((M.reference||REFNAME)?(' '+esc(M.reference||REFNAME)+''):'')+(AA2LBL?(' ↔ '+esc(AA2LBL)+''):''); + var refhdr=(isIndel?'Indel':'SNP')+((M.reference||REFNAME)?(' '+esc(M.reference||REFNAME)+''):'')+(AA2LBL?(' ↔ '+esc(AA2LBL)+''):''); var nameRow=''+refhdr+''+ vi.map(function(i){var s=samples[i];return ''+esc(s)+'';}).join('')+''; function rowHTML(r){ - var lbl=''+esc(r.gene||r.contig)+''+geneRvTag(r.gene)+' '+refPos(r.pos,r.pos_h37rv)+' '+esc(r.ref)+'→'+esc(r.alt)+(r.aa?(' '+aaDual(r.aa,r.aa_h37rv)+''):''); + 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?'+':'−')+Math.abs(r.len)+' bp'; + lbl=''+esc(r.gene||r.contig)+''+geneRvTag(r.gene)+' '+refPos(r.pos)+' '+mxAllele(r.ref)+'→'+mxAllele(r.alt)+ + ' '+ln+(fs?' · frameshift':'')+''+ + (r.aa?(' '+esc(r.aa)+''):''); + } else { + lbl=''+esc(r.gene||r.contig)+''+geneRvTag(r.gene)+' '+refPos(r.pos,r.pos_h37rv)+' '+esc(r.ref)+'→'+esc(r.alt)+(r.aa?(' '+aaDual(r.aa,r.aa_h37rv)+''):''); + // 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+=' codon '+esc(r.mnv_aa.join(', '))+''; + } var cells=vi.map(function(i){var c=r.cells[i], s=samples[i]; - if(!c) return ''; + if(!c) return ''; var afTxt=c[0].toFixed(2).replace(/^0/,'').replace(/^1\.00$/,'1'); var dpTxt=(showDP&&c[1]!=null)?(''+c[1]+''):''; - return ''+afTxt+''+dpTxt+''; + 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 ''+afTxt+''+dpTxt+''; }).join(''); return ''+lbl+''+cells+''; } @@ -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.length2Sample QC
Flagged samplesAll samplesLineagesContaminationDistributionsaDNA damage
3Genome & genes
ConsensusGenome landscapeDeletionsFunctional impactGene burdenVariable genesSelection (dN/dS)
4Relatedness
DistancesClustersOutside their group
-
5Variants over time
Gained per seriesMinority variantsSNP dynamicsCo-varying pairsSNP matrixVariant × dose
+
5Variants over time
Gained per seriesMinority variantsSNP dynamicsCo-varying pairsVariant matrixVariant × dose
6Resistance
Resistance mutationsAlong each series
7Gene conversion
Candidate tracts
8Diagnostics
Metric pairsCorrelation matrixQC spaceDivergenceDose × treatmentSampling dates
@@ -123,7 +123,7 @@

Co-varying pairs variants whose trajectories move together (concordant) or against each other (discordant) within a series

-

SNP matrix every SNP site × sample, as allele frequency; filter by gene or position and download the full matrix

+

Variant matrix every SNP site, or every indel, × sample, as allele frequency; filter by gene or position and download the full matrix

Variant × dose per-variant allele frequency against dose (Spearman ρ with Benjamini–Hochberg FDR); click a row to plot it

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/`, `_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 `_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 `_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..9e82032 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(.tbi) # -> The same as VCF: a record per SNP and one per MNV (INFO/MR = its reads) + │ └── MP00091.LENS.mnv.summary.json / .manifest.json # -> Counts, command, version and checksums of the run │ ├── lineage/ # 🧬 Pathotypr typing (only if --run_pathotypr) │ └── MP00091__.pathotypr.lineage_summary.tsv # -> Pathotypr lineage / sub-lineage call (__ 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 + +`_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 `..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). +`_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 `_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 × 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 × dose. | | **6 · Resistance** | WHO-catalogue mutations grouped by mutation, the sample × 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 × 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..dd5f52f --- /dev/null +++ b/modules/mnv.nf @@ -0,0 +1,91 @@ +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 + tuple val(sampleId), val(refId), path("${sampleId}.${refId}.MNV.vcf.gz"), path("${sampleId}.${refId}.MNV.vcf.gz.tbi"), 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 --index-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.vcf.gz.tbi + 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) // samplereferenceget_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/..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 ..indels(.ann).vcf.gz, marked PASS or LowSupport, and the run's in the matrix. + make_indel_matrix = true // _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 _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, ..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, ..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=\n" + "##FORMAT=\n" + "##FORMAT=\n" + "##FORMAT=\n" + "##FORMAT=\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=\n##FORMAT= 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 " || "`, 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.""" From d9a667c37804ea58f39c7ba85071caaf4181e30b Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Thu, 24 Sep 2026 20:34:17 +0200 Subject: [PATCH 2/2] Publish the MNV VCF unindexed: get_MNV's image has no tabix The first real run of GET_MNV failed on every sample. get_MNV writes its VCF and indexes it with tabix, which the bioconda image it runs in does not carry, so it warned, skipped the index and exited 0, and Nextflow failed the task over the .tbi it had been told to expect. A stub run touches every declared output, so only a real run could see it; an end-to-end run did. The VCF is now published bgzipped and unindexed (tabix -p vcf indexes it as it is), and the codon table the run reads never used the index. --- docs/outputs.md | 2 +- modules/mnv.nf | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/outputs.md b/docs/outputs.md index 9e82032..ca517d2 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -55,7 +55,7 @@ results_bampiro/ │ ├── 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(.tbi) # -> The same as VCF: a record per SNP and one per MNV (INFO/MR = 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) diff --git a/modules/mnv.nf b/modules/mnv.nf index dd5f52f..1d58df8 100644 --- a/modules/mnv.nf +++ b/modules/mnv.nf @@ -25,7 +25,9 @@ process GET_MNV { output: tuple val(sampleId), val(refId), path("${sampleId}.${refId}.MNV.tsv"), emit: tsv - tuple val(sampleId), val(refId), path("${sampleId}.${refId}.MNV.vcf.gz"), path("${sampleId}.${refId}.MNV.vcf.gz.tbi"), emit: vcf + // 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 @@ -47,7 +49,7 @@ process GET_MNV { 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 --index-vcf-gz \\ + --threads ${task.cpus} --both --vcf-gz \\ --summary-json ${sampleId}.${refId}.mnv.summary.json \\ --run-manifest ${sampleId}.${refId}.mnv.manifest.json """ @@ -56,7 +58,6 @@ process GET_MNV { """ touch ${sampleId}.${refId}.MNV.tsv touch ${sampleId}.${refId}.MNV.vcf.gz - touch ${sampleId}.${refId}.MNV.vcf.gz.tbi touch ${sampleId}.${refId}.mnv.summary.json touch ${sampleId}.${refId}.mnv.manifest.json """