Codon-level variant annotation
A VCF lists bases.
get_MNV reads codons.
Two substitutions in one codon make one amino acid, not two. get_MNV translates every codon whole, places each indel against its gene and, given a BAM or CRAM, counts the reads that carry the changes together.
- In
- a VCF, BCF or iVar table, a FASTA, a GFF or gene table, and optionally a BAM or CRAM
- Out
- a TSV, a VCF with
ANN, and a report that opens offline
4 of 4 · The whole codon
example/G35894.var.snp.vcf
104941CT
104942CG
104944AG
Rv0095c, as the gene reads minus strand: positions count down
codon 91 · Val1 change
GCC Ala · one SNP row
codon 92 · Gly2 changes
104942 aloneCGGArg
104941 aloneGAGGlu
both, one moleculeCAGGln
codon 93 · Hisno change
- The VCF lists three substitutions, one record each, written on the reference's forward strand.
- Grouped by codon. The gene runs on the minus strand, so C>T at 104941 is G>A in the gene. 104941 and 104942 share codon 92; 104944 sits in codon 91 on its own.
- One base at a time, codon 92 reads Arg and Glu. That is what a
base-by-base annotator reports: the input's own snpEff
ANNsays p.Gly92Arg and p.Gly92Glu. - The whole codon reads Gln: the amino acid the protein carries when both changes sit on one molecule, which a BAM can confirm read by read. get_MNV writes it as one SNP/MNV row, p.Gly92Gln, with the two single readings beside it.
example/G35894.var.snp.vcf.Try it: one codon, two readings
Change any base. Each change is read on its own, the way a base-by-base annotator reads it, and then all together, the way get_MNV reads the codon.
The script that makes this interactive did not run. Here is the first example, worked out.
| Read as | Codon | Amino acid |
|---|---|---|
| Reference | CGT | Arg |
| Position 1 alone, C>T | TGT | Cys missense |
| Position 2 alone, G>C | CCT | Pro missense |
| Both, as get_MNV reads it | TCT | Ser missense |
One base at a time this codon gives Cys and Pro. Read whole, it gives Ser, which neither of them names.
Translated with the bacterial code (table 11), get_MNV's default; nine NCBI tables are supported. This is a codon inside a gene: the first codon of a gene is read as an initiator instead.
What it does with a VCF
Steps 1 to 3 need no reads. A BAM or CRAM only adds step 4, which turns "these two changes are in one codon" into "these two changes are on one molecule".
- 1
Decompose
Each record's REF and ALT into the changes it holds: SNVs, MNVs, insertions, deletions, delins.
- 2
Place
Each change in its gene and codon, from a GFF or a gene table, on either strand and across exons.
- 3
Translate
The whole codon, in one of nine NCBI genetic codes, and each indel against the protein.
- 4
Count
With reads: how many carry each change, on which strand, and whether the changes travel together.
- 5
Write
A TSV, a VCF with
ANN, BCF on request, and a report that opens offline.
How get_MNV tells one molecule from two · Indels and local haplotypes
Install
One binary, no C dependencies, on macOS, Linux and Windows. A desktop app does the same with drag and drop.
Bioconda
conda install -c bioconda get_mnvFrom source
git clone https://github.com/PathoGenOmics-Lab/get_MNV.git
cd get_MNV
cargo install --path .Rust 1.91 or newer.
Desktop app
macOS, Linux and Windows builds on the Releases page, with a guided tutorial.
Not signed with an Apple certificate: on the first launch, right-click the app and choose Open.
Check what you downloaded
Every release carries a SHA256SUMS file covering all of its assets. Download it
next to the file you want and check them together; each file it finds prints OK.
sha256sum --ignore-missing -c SHA256SUMSOn macOS, shasum -a 256 --ignore-missing -c SHA256SUMS does the same.
Your first run
The repository ships a M. tuberculosis sample: a VCF, its reference and a gene table. From a clone:
get_mnv \
--vcf example/G35894.var.snp.vcf \
--fasta example/MTB_ancestor.fas \
--genes example/anot_genes.txt \
--both --output-dir out \
--report report.html- 950 records in, 941 rows out: 797 SNP, 10 SNP/MNV codons and 134 intergenic
out/G35894.var.snp.MNV.tsvandout/G35894.var.snp.MNV.vcfreport.html, the example report
What you get
Codon 92 of Rv0095c from Figure 1, in each of the three outputs.
TSV one row per codon or indel
- Positions
104941, 104942- AA Changes
Gly92Gln- SNP AA Changes
Gly92Glu, Gly92Arg- Variant Type
SNP/MNV- HGVS c.
c.[274G>C;275G>A]- DBS Class
CC>TG
VCF each base, and the codon
104941 C T …AA=Gly92Glu…
104942 C G …AA=Gly92Arg…
104941 CC TG …AA=Gly92Gln;…
ANN=TG|missense_variant|MODERATE|
Rv0095c_Rv0095c|…|p.Gly92Gln|
…|INFO_COMPOUND_ANNOTATIONANN is in snpEff's format, so bcftools +split-vep reads it.
Report one HTML file, no server
Open the example report
Read on
Where to go once the first run works.
Command line tutorial
A first run end to end on the bundled data, the output explained line by line.
Common recipes
A cohort, an iVar table, a compressed VCF, a report from earlier runs.
CLI reference
Every option with its default, and what each one changes about the answer.
Output formats
Each TSV column, VCF INFO key and JSON field, and when a value is absent rather than zero.
Linkage
Telling a real haplotype from two variants that only share a codon, using the reads.
Desktop GUI
The same analysis with drag and drop, a read viewer and the report built in.
Troubleshooting
The input mismatches that stop a run, and what each message asks for.
Glossary
SNV, MNV, codon phase, cis and trans, in the sense these pages use them.