Skip to content
BranchesTags

Repository files navigation

MitoQuest: Human Mitochondrial sequencing data Analysis Toolkit

MitoQuest is a cross-platform, efficient and practical bioinformatics toolkit written in C++ (C++17) for analyzing human mitochondrial DNA (mtDNA) from whole-genome / mtDNA-targeted sequencing (WGS) data. It calls mitochondrial SNVs and Indels, quantifies heteroplasmy and homoplasmy, estimates mtDNA copy number, and provides downstream utilities for VCF sub-sampling, annotation, and quality recalibration.

MitoQuest is built on top of htslib (vendored as a submodule) and is designed to scale from a single-sample clinical workflow to population-level cohorts with thousands of samples.

mitoquest: Human Mitochondrial sequencing data Analysis Toolkit
Version: 1.11.0

Usage: mitoquest <command> [options]
Commands:
  caller       Mitochondrial variants and heteroplasmy/homoplasmy caller.
  subsam       Extract mitochondrial variants for specified samples from VCF files and output a new VCF file.
  copynum      Estimate per-chromosome (incl. mtDNA) relative copy number from a BAM/CRAM file.
  trans-prep   Extract mother-child mtDNA allele transmission pairs from a multi-sample VCF + FAM file.
  ne-estimate  Estimate the mtDNA bottleneck size (Ne) from transmission pairs via Beta-Binomial MMLE (Maximum Marginal Likelihood).

In addition to the main mitoquest binary, the project ships:

  • tools/ — a suite of Python helper scripts for VCF QC, annotation, pipeline assembly, and downstream analysis.
  • data/ — curated reference resources (population databases, in-silico predictors, RNA/protein domain annotations, blacklist, Phylotree variants, …) used by the annotation and QC tools.

Installation and quick start

Pre-built static binaries are published on the GitHub Releases page — most user shoild simply download the binary and run.

Platform Download Notes
Linux (x86_64) mitoquest-linux-static Requires glibc ≥ 2.35 (see below)
macOS (arm64 / Intel) mitoquest-macos-static Requires macOS 12+

System requirements for mitoquest-linux-static

The Linux binary is a partial-static build, produced on Ubuntu 22.04 (glibc 2.35) in CI. It bundles libstdc++, libgcc, htslib, zlib, bzip2, xz, and openssl statically — only the system C library (glibc) is linked dynamically. Because glibc symbol versions are forward-compatible only, the binary requires the host glibc version to be ≥ 2.35.

Quick check on your machine:

# If the printed glibc version is >= 2.35, mitoquest-linux-static will run.
ldd --version | head -1

If you see this — or you are on CentOS / RHEL / Rocky / AlmaLinux / older Ubuntu / older Debian — you can compile from source instead (see below).

# Linux
wget https://github.com/ShujiaHuang/mitoquest/releases/latest/download/mitoquest-linux-static
chmod +x mitoquest-linux-static
mv mitoquest-linux-static motiquest
./mitoquest --help
# macOS
curl -LO https://github.com/ShujiaHuang/mitoquest/releases/latest/download/motiquest-macos-static
chmod +x mitoquest-linux-static motiquest
mv mitoquest-linux-static motiquest
./mitoquest --help

Important

Rename the downloaded binary to mitoquest for convenience. You may also move it to a directory in your $PATH (e.g. /usr/local/bin) for system-wide access.

Tip

If the pre-built binary does not work on your system (e.g. glibc < 2.35 on older Linux), you can compile from source instead. Detailed build instructions (CMake build, static build, and manual g++ fallback) are available in docs/INSTALL_FROM_SOURCE.md.

Commands overview

Usage: mitoquest <command> [options]

Commands:
  caller       Call mitochondrial variants (SNVs/Indels) and quantify
               heteroplasmy / homoplasmy from BAM/CRAM files.
  subsam       Extract a subset of samples from an mtDNA VCF and recompute
               INFO fields.
  copynum      Estimate per-chromosome (incl. mtDNA) relative copy number
               from a BAM/CRAM file.
  trans-prep   Extract mother-child mtDNA allele transmission pairs from a
               multi-sample VCF + a PLINK FAM file.
  ne-estimate  Estimate the mtDNA bottleneck size (Ne) from transmission
               pairs via the Beta-Binomial Maximum Marginal Likelihood Estimator (MMLE).

mitoquest caller — Variant calling and heteroplasmy/homoplasmy detection

mitoquest caller reads aligned reads from one or more BAM/CRAM files and emits a multi-sample VCF containing per-sample heteroplasmy fractions (HF) and per-site INFO fields suitable for downstream filtering and annotation.

Full parameter reference of caller

Usage: mitoquest caller [options] -f ref.fa -o output.vcf.gz in1.bam [in2.bam ...]

Required options:
  -f, --reference FILE       Reference FASTA file
  -o, --output    FILE       Output VCF file (use .vcf.gz for bgzipped output)

Optional options:
  -b, --bam-list FILE        List of input BAM/CRAM filenames, one per line.
  -Q, --min-BQ INT           Skip bases with base quality smaller than INT
                             (default: 20).
  -q, --min-MQ INT           Skip alignments with mapQ smaller than INT
                             (default: 20).
  -r, --regions REG[,...]    Comma-separated regions to process
                             (default: entire mtDNA contig).
                             Format: chr | chr:start | chr:start-end
                             Example: chrM or chrM:1-1000,chrM:8000-8200
  -p, --pairs-map-only       Only use paired reads where mate maps to the
                             same chromosome.
  -P, --proper-pairs-only    Only use properly paired reads (SAM flag 0x2).
  --filename-has-samplename  When BAM/CRAM filenames are 'SampleID.xxxx.bam',
                             skip reading sample names from the BAM header.
                             Saves significant time on large cohorts.
  -j, --het-threshold FLOAT  Heteroplasmy fraction threshold below which an
                             allele is treated as reference (default: 0.01).
  -c, --chunk INT            Chunk size (bp) for parallel region processing
                             (default: 1000).
  -t, --threads INT          Number of threads (default: all available CPUs).
  -h, --help                 Show this help message and exit.

Usage examples

Single-sample variant calling on the entire mtDNA contig:

mitoquest caller \
    -f reference.fasta \
    -o sample.mt.vcf.gz \
    sample.bam

Multi-sample call from a list of BAM/CRAM files (with quality filters):

mitoquest caller \
    -f reference.fasta \
    -o cohort.mt.vcf.gz \
    -b bamfile.list \
    -Q 30 -q 30 -t 16

Recommended call with sample-name optimisation (large cohorts):

# When BAMs are named `SampleID.bam` / `SampleID.cram`, --filename-has-samplename
# avoids opening every BAM just to read @RG SM tags.
mitoquest caller \
    -f reference.fasta \
    -o cohort.mt.vcf.gz \
    -Q 30 -q 30 -t 24 \
    --filename-has-samplename \
    -b bamfile.list

Restrict calling to a specific region (e.g., the D-loop / control region):

mitoquest caller \
    -f reference.fasta \
    -o cohort.dloop.vcf.gz \
    -r chrM:1-576,chrM:16024-16569 \
    -Q 30 -q 30 -t 16 \
    --filename-has-samplename \
    -b bamfile.list

Stricter heteroplasmy detection (lower the HF threshold):

# Default heteroplasmy threshold is 0.01 (1%).  Lower it to 0.005 to capture
# very low-frequency heteroplasmy (use with caution; requires high coverage).
mitoquest caller \
    -f reference.fasta \
    -o cohort.lowhet.vcf.gz \
    -j 0.005 \
    -Q 30 -q 30 -t 24 \
    --filename-has-samplename \
    -b bamfile.list

Use only properly paired reads (more stringent for NUMT filtering):

mitoquest caller \
    -f reference.fasta \
    -o cohort.proper.vcf.gz \
    -P --filename-has-samplename \
    -Q 30 -q 30 -t 16 \
    -b bamfile.list

Mix a file list with extra BAMs on the command line:

mitoquest caller \
    -f reference.fasta \
    -o cohort.vcf.gz \
    -b bamfile.list \
    -Q 30 -q 30 -t 16 \
    extra1.cram extra2.bam

mitoquest subsam — Extract samples from an mtDNA VCF

Extract a subset of samples from a multi-sample VCF/BCF file, optionally recomputing INFO fields (AC/AN/AF/NS, …) so they reflect only the kept samples.

Full parameter reference of subsam

Usage: mitoquest subsam [options] -i <input.vcf> -o <output.vcf> [-s <samplelist>] [<sample1> ...]

Options:
  -i, --input FILE     Input VCF/BCF file (required).
  -o, --output FILE    Output VCF/BCF file (required).
  -s, --sample FILE    File listing sample names to keep, one per line.
  -O, --output-type    Output file type: v|z|b|u
                       (v: VCF, z: bgzipped VCF, b: BCF, u: uncompressed BCF).
                       Default: inferred from output filename extension.
  --no-update-info     Do not recalculate AC/AN/AF/NS in the INFO column.
  --keep-all-site      Retain sites that become reference-only after
                       subsetting (default: drop them).
  -h, --help           Show this help message and exit.

Usage examples for extracting samples

Extract samples listed in a file (output bgzipped VCF):

mitoquest subsam \
    -i cohort.mt.vcf.gz \
    -o subset.mt.vcf.gz \
    -s sample_names.txt

Extract two specific samples directly from the command line:

mitoquest subsam \
    -i cohort.mt.vcf.gz \
    -o pair.mt.vcf \
    -O v \
    SampleA SampleB

Keep all sites (including those that become ref-only) and skip INFO updates:

# Useful when you want to preserve a fixed multi-cohort site list, e.g. for
# downstream join-call merging.
mitoquest subsam \
    -i cohort.mt.vcf.gz \
    -o subset.mt.vcf.gz \
    -s sample_names.txt \
    --keep-all-site --no-update-info

mitoquest copynum — mtDNA copy-number estimation

mitoquest copynum estimates per-chromosome relative copy number from a sorted/indexed BAM or CRAM file. The autosomal chromosomes serve as the diploid baseline (CN = 2); the mitochondrial chromosome is reported on the same scale, i.e. the expected number of mtDNA molecules per diploid cell. The output is a TSV with fragment counts, GC content, length-normalized fragment ratio, and the copy-number mean + 95% confidence interval for every contig in the BAM header.

Full parameter reference

Usage: mitoquest copynum [options] <input.bam/cram>

Options:
  -r, --reference FILE   Reference genome FASTA file (required; needed for
                         CRAM decoding and GC content calculation).
  -o, --output    FILE   Output TSV file (default: stdout).
  -q, --mapq      INT    Minimum mapping quality score [0].
  -t, --threads   INT    Number of worker threads [hardware_concurrency].
  -s, --seqtype   STR    Sequencing type: auto | pe | se [auto].
                         pe = paired-end, se = single-end.
  -L, --regions   STR    Restrict counting / GC / length normalization to
                         the listed regions, given either as a comma-
                         separated list (e.g. 'chrM:1-300,chrM:16000-16569')
                         or as a path to a file (one region per line, in
                         either 'chr:start-end' samtools form or BED-style
                         'chr<TAB>start<TAB>end' triples; '#' starts a
                         comment). The intervals REPLACE the whole-
                         chromosome window for the chromosomes they cover;
                         chromosomes not mentioned keep their full-length
                         behaviour. Typical use: exclude NUMT-affected zones
                         from the mtDNA copy-number estimate.
  -h, --help             Show this help message and exit.

Usage examples for estimating copynum

Estimate copy numbers from a BAM file (paired-end auto-detected):

mitoquest copynum -r reference.fasta sample.bam > sample.cn.tsv

Estimate from a CRAM file with a stricter MAPQ filter and 8 threads:

mitoquest copynum \
    -r reference.fasta \
    -q 30 -t 8 \
    -o sample.cn.tsv \
    sample.cram

Force single-end counting (useful for legacy unpaired data):

mitoquest copynum \
    -r reference.fasta \
    -s se -q 30 -t 4 \
    sample.bam > sample.cn.tsv

Exclude NUMT-affected zones via inline regions (recommended for mtCN):

# Restrict chrM measurement to the two stretches least affected by NUMTs;
# autosomes are still measured over their full length, so the diploid
# baseline is unchanged.
mitoquest copynum \
    -r reference.fasta \
    -L 'chrM:1-300,chrM:16000-16569' \
    -q 30 -t 8 \
    sample.bam > sample.cn.tsv

Exclude NUMT zones via a BED file (one region per line):

# my_chrM_regions.bed (0-based half-open, standard BED):
#   chrM    0       300
#   chrM    15999   16569
mitoquest copynum \
    -r reference.fasta \
    -L my_chrM_regions.bed \
    -q 30 -t 8 \
    sample.bam > sample.cn.tsv

The TSV output has one row per contig in the BAM header, with the following columns:

#Chromosome  Fragments  Chrom_Length  GC_Content  Fragment_Normalized_Ratio  CopyNum  CopyNum-CI95-Lower  CopyNum-CI95-Upper  Effective_Length  Regions_Used
  • Chrom_Length always reports the full contig length from the BAM header.
  • Effective_Length is the number of bases that were actually counted (equals Chrom_Length when -L did not target that contig; otherwise it is the sum of the merged user-supplied intervals on that chromosome).
  • Regions_Used is . for unrestricted chromosomes, or a comma-separated list of merged start-end intervals (1-based inclusive) when -L was applied to that chromosome.
  • When -L is supplied, the output also includes a #Regions argument: ... header comment recording the original CLI value for reproducibility.

Auxiliary Python tools (tools/)

The tools/ directory ships a set of Python helper scripts that complement the C++ binaries:

Script Purpose
tools/mito_annotate.py Annotate a mitoquest caller VCF with population, in-silico, and disease databases.
tools/parse_annotatedVCF.py Convert an annotated VCF into a tidy table for downstream statistics.
tools/mtDNA_variant_QC.py Python QC prototype (superseded by mitoquest variant-qc).
tools/mtDNA_vcf_to_table.py Flatten a multi-sample VCF into long-format TSV (one row per sample × site).
tools/parse_vcf.py Generic VCF parsing helper (used by other tools).
tools/filter_mergedVCF.py Apply hard filters on a merged mtDNA VCF.
tools/merge_cr_ncr_vcf.py Merge coding-region and non-coding-region VCFs into a single file.
tools/rewrite_vcf.py Rewrite a VCF (header normalisation, CHROM renaming, etc.).
tools/shift_fasta.py Produce a circularly shifted FASTA (used by the join-region pipeline).
tools/create_join_seq.py Build the joined coding-region / non-coding-region reference for re-alignment.
tools/detect_NUMT_by_mtCN.py Flag potential NUMT contamination using mtCN ratios per sample.
tools/plot_bottleneck_simulation.py Plot per-bin observed drift vs. analytical Kimura prediction p_m(1−p_m)/Ne from ne-estimate --bin-simulation.
tools/plot_ne_profile.py Plot the MMLE and Kimura objective curves over Ne from ne-estimate --ne-profile.
tools/vcf_format_validator.py Sanity-check a VCF for downstream compatibility.

Each script supports -h / --help.

A typical end-to-end workflow

# 1. Variant calling
mitoquest caller \
    -f rCRS.fasta \
    -o cohort.raw.vcf.gz \
    -Q 30 -q 30 -t 24 \
    --filename-has-samplename \
    -b bamfile.list

# 2. Annotation against MITOMAP / HelixMTdb / gnomAD / in-silico predictors
python tools/mito_annotate.py \
    -i cohort.raw.vcf.gz \
    -o cohort.annotated.vcf.gz \
    --resource-dir data

# 3. Bayesian quality control (filter true mutations from artefacts)
mitoquest variant-qc \
    -i cohort.raw.vcf.gz \
    -o cohort.qc.vcf.gz \
    -t cohort.qc.report.tsv \
    --dp-threshold 200 --hq-threshold 30

# 4. mtDNA copy-number estimation per sample
while IFS= read -r bam; do
    mitoquest copynum -r reference.fasta -q 30 -t 4 "$bam"
done < bamfile.list > cohort.mtCN.tsv

# 5. Convert to a long-format table for analysis in R/Python
python tools/mtDNA_vcf_to_table.py \
    -i cohort.qc.vcf.gz \
    -o cohort.long.tsv

Reference resources (data/)

The data/ directory contains curated reference resources used by the annotation and QC tools. See data/README.md and data/UPDATE_LOG.md for sources and update history.

Highlights:

  • Population databases: HelixMTdb, gnomAD v3.1 (chrM), dbSNP (chrM), MITOMAP polymorphisms / disease, ClinVar (chrM SNVs).
  • In-silico predictors: MitoTIP, MitImpact, HmtVar, t-APOGEE, phyloP100way.
  • Functional annotations: tRNA secondary-structure positions, rRNA bridge bases, RNA modifications, UniProt protein domains, Complex-I proton-pump residues.
  • Phylogenetic context: Phylotree variants, human-shifted chimp alignment, mtDNA gene loci.
  • Blacklist: known systematic-error sites (data/blacklist.txt).

Tips and best practices

  • Reference FASTA: Use the rCRS sequence (NC_012920.1) as the reference whenever possible — most downstream databases and tools (MITOMAP, HelixMTdb, MitoTIP) are coordinated to it.
  • --filename-has-samplename: If your BAM files are named {SampleID}.bam or {SampleID}.cram, always set this flag — it avoids reading every BAM header and can save hours on large cohorts.
  • -j / --het-threshold: The default 0.01 (1 %) is a sensible balance between sensitivity and false-positive rate at typical WGS depths (≥30× on chrM). Lower it to 0.005 only when coverage is very high (≥1000×) and you are explicitly targeting low-frequency heteroplasmy.
  • -r / --regions: For large multi-sample analyses, you can split chrM into chunks and run several mitoquest caller instances in parallel; each instance is multithreaded internally via -t.
  • Output compression: Always use .vcf.gz as the output filename — mitoquest automatically writes bgzipped, tabix-ready output when the extension matches.
  • NUMT filtering: For samples with elevated mtCN ratios (mitoquest copynum output much higher than expected), consider either:
    • restricting the mtDNA measurement to NUMT-poor stretches with mitoquest copynum -L 'chrM:1-300,chrM:16000-16569' ... (chrM is measured only on the listed intervals; autosomes stay full-length),
    • and/or running the caller with -P/--proper-pairs-only and inspecting outputs with tools/detect_NUMT_by_mtCN.py.

Development

MitoQuest is under active development. To update to the latest version:

git pull
git submodule update --recursive
cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --parallel

To rebuild from scratch (recommended after a major htslib update):

rm -rf build
(cd htslib && make distclean) || true
cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --parallel

Pull requests, bug reports and feature requests are very welcome at https://github.com/ShujiaHuang/mitoquest.

"What I cannot create, I do not understand." — Richard Feynman

About

Human Mitochondrial sequencing data Analysis Toolkit. Fetching information like variants, heteroplasmy/homoplasmy, etc. on mitochondrial genome from sequencing data

Topics

bioinformaticsbioinformatics-toolheteroplasmyhomoplasmymitochondriamitochondrial-variantsnext-generation-sequencingvariant-callersvariant-callingwhole-genome-sequencing

Resources

Readme
AGPL-3.0 license
Activity

Stars

5 stars

Watchers

1 watching

Forks

0 forks
Report repository

Releases

Packages

Used by

Contributors

Languages