GPTomics/bio-consensus-sequences

Version Compatibility

Reference examples tested with: bcftools 1.19+, samtools 1.19+, bedtools 2.31+, iVar 1.4+, minimap2 2.26+, BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: the -H argument vocabulary (N/R/A/I/LR/LA/SR/SA/NpIu, where I = IUPAC code for all genotypes) has grown across bcftools 1.x releases; IUPAC output is available both as the -H I code and the standalone -I/--iupac-codes flag. Always confirm the accepted letters with bcftools consensus on the installed version before scripting a selection.

Consensus Sequences

"Generate a consensus sequence from my VCF" -> Apply called variants onto a reference FASTA, producing a sample-specific sequence, with a deliberate choice of haplotype projection and no-coverage masking.

• CLI (from a VCF): bcftools consensus -f reference.fa input.vcf.gz
• CLI (viral/amplicon from a BAM): samtools mpileup ... | ivar consensus -p out
• Python: cyvcf2 + Bio.SeqIO for SNP-only prototypes

The governing principle

bcftools consensus walks the reference and substitutes ALT alleles at the positions present in the VCF. Everything else is copied from the reference verbatim -- which drives three traps that ruin more consensus analyses than any tool bug:

1. A consensus silently emits REFERENCE wherever the VCF is silent -- including positions with zero coverage. No data and confidently-reference look identical in the output. An unmasked consensus therefore manufactures false confidence at exactly the sites where the sample was never observed. Mask no-coverage sites (below) or the FASTA lies.
2. -H 1 on an UNPHASED VCF yields a chimeric pseudo-haplotype. Haplotype selection is only meaningful when genotypes are phased; on unphased data it mixes alleles from different real chromosomes into a sequence that exists in no cell. Verify | phasing before selecting a haplotype.
3. A single FASTA cannot faithfully represent a diploid genome. Every projection (-H 1, -I, -H A) is lossy in a different way; for phase-sensitive work keep the VCF, not the consensus.

The input VCF must be bgzipped and indexed (bgzip + bcftools index/tabix); plain-gzip or unindexed input errors out. The REF bases in the VCF must match the FASTA exactly or bcftools warns and skips those records. Normalize first (see Normalization).

Basic Usage

bcftools consensus reads variants from a bgzipped, indexed VCF and writes FASTA:

bcftools index input.vcf.gz                                    # .csi index (or tabix -p vcf)
bcftools consensus -f reference.fa input.vcf.gz > consensus.fa
bcftools consensus -f reference.fa -o consensus.fa input.vcf.gz  # -o instead of redirect

For a multi-sample VCF, always pass -s -- without it, the applied genotypes are undefined:

bcftools query -l input.vcf.gz                                 # list samples
bcftools consensus -f reference.fa -s sample1 input.vcf.gz > sample1.fa

Restrict to a region with -r (the FASTA header is then >chr:from-to):

bcftools consensus -f reference.fa -r chr1:1000000-1010000 -s sample1 input.vcf.gz > gene.fa

Haplotype Selection and the Phasing Trap

-H chooses which allele to apply from FORMAT/GT. The codes are case-insensitive:

| Option | Applies | Use when | |--------|---------|----------| | -H 1 / -H 2 | Allele at GT index 1 or 2 | Emitting one true chromosome -- only valid on PHASED genotypes | | -H A | ALT allele in every genotype | Maximum divergence from reference; a chimera of both chromosomes | | -H R | REF allele at heterozygous sites | Conservative consensus; discards het ALT alleles | | -H I (or the standalone -I / --iupac-codes flag) | IUPAC ambiguity code | Retain heterozygosity in one sequence (see caveat below) | | -H LA/LR/SA/SR | Longer/shorter allele, tie broken by ALT/REF | Length-driven selection; confirm the letter set on the installed version |

The chimeric-haplotype footgun. -H 1/-H 2 are only meaningful when genotypes are phased (0|1, pipe separator). With unphased genotypes (0/1, slash), the assignment of "which allele is haplotype 1" is arbitrary per site, so -H 1 across many heterozygous sites produces a switch-error mosaic that corresponds to no real chromosome -- while looking like a clean haplotype FASTA. This is the single most dangerous consensus mistake. Verify phasing before any -H 1/-H 2:

bcftools query -f '%CHROM\t%POS[\t%GT]\n' input.vcf.gz | head   # phased: 0|1 ; unphased: 0/1

If genotypes are unphased, phase first (read-backed WhatsHap/HapCUT2, trio, statistical SHAPEIT/Eagle -- accurate for common variants, poor for rare/singletons -- or native long-read phasing). See phasing-imputation/haplotype-phasing and variant-calling/vcf-basics for GT interpretation.

What a Consensus Cannot Represent

A single consensus FASTA is a lossy projection of a diploid genome; the right projection depends on the downstream use, and some tasks need the VCF instead:

| Strategy | Flag | Best for | Loses | |----------|------|----------|-------| | Two haplotype sequences | -H 1 + -H 2 (phased) | Allele-specific expression, compound-het, HLA, cis-regulatory haplotypes | Nothing (if correctly phased) | | IUPAC ambiguity codes | -I | Retaining het signal in one sequence | Phase/linkage; many tree/alignment tools read IUPAC as N | | All ALT alleles | -H A | Max divergence, quick draft | Reality -- exists in no cell | | REF at het sites | -H R | Conservative single sequence | Every heterozygous ALT allele |

Two hard boundaries:

• For phase-sensitive work, keep the VCF (or two phased haplotype FASTAs), not a single consensus. Collapsing hets to IUPAC or picking one allele discards linkage that the analysis needs -- treating a consensus FASTA as "the sample's genome" for compound-het or allele-specific analysis is a category error.
• bcftools consensus cannot apply symbolic SV alleles (<DEL>, <INS>, <DUP>, <INV>): those carry no ALT sequence, only INFO fields, so consensus has nothing to substitute. Short-read SV VCFs (Manta/DELLY) are mostly symbolic and are NOT directly consensus-able. Folding SVs into a consensus needs sequence-resolved records (long-read/assembly callers emit these) or an assembly-based approach -- see variant-calling/structural-variant-calling.

For phylogenetics specifically, prefer one clean phased haplotype or a homozygous-ALT-only sequence over IUPAC, because ambiguity codes are silently dropped by many tree builders:

bcftools view -i 'GT="AA"' input.vcf.gz | bcftools consensus -f reference.fa > hom_alt.fa

Masking No-Coverage Sites (the load-bearing footgun)

Because unobserved positions are emitted as reference (trap 1), a consensus must mask sites with insufficient data. -m mask.bed replaces the listed regions (default char N via --mask-with N). The mask must be built from callable depth, and the depth step hides a silent bug:

samtools depth WITHOUT -a OMITS zero-coverage positions from its output -- so those positions never enter the low-depth BED, never get masked, and stay as reference: the exact false-confidence failure the mask was meant to prevent. Always use -a (report all positions) so no-coverage sites are captured:

# Build a mask of every position below the callable-depth threshold. -a is mandatory:
# without it, zero-coverage positions are absent from the output and escape masking.
samtools depth -a aligned.bam | awk '$3 < 10 {print $1"\t"$2-1"\t"$2}' | bedtools merge > lowcov.bed

bcftools consensus -f reference.fa -m lowcov.bed input.vcf.gz > consensus.fa

The < 10 threshold is a minimum-callable-depth policy (10x is a common floor for confident base calls); set it to the depth below which the calls are not trusted. bedtools genomecov -bga -ibam aligned.bam is an equivalent zero-coverage-aware alternative that also emits 0-depth intervals.

Do NOT rely on -M/-a for this: -M N outputs N only for missing ./. genotypes already present in the VCF, and -a N replaces every position absent from the VCF (which N-outs the entire non-variant genome). Neither distinguishes no-coverage from confident-reference -- only a depth-derived mask does.

Normalization Before Consensus

Goal: Apply indels at the correct reference position and sequence.

Approach: Left-align and split multiallelics with bcftools norm so each record matches the reference context; consensus applies records positionally and mis-represented indels corrupt the output.

bcftools norm -f reference.fa input.vcf.gz -Oz -o norm.vcf.gz
bcftools index norm.vcf.gz
bcftools consensus -f reference.fa norm.vcf.gz > consensus.fa

Un-normalized or overlapping indels produce wrong sequence, and bcftools consensus only warns to stderr while still emitting output -- so the corruption is silent unless the stderr is inspected. Even after norm, two records whose REF spans collide remain a hazard; grep the run for warnings and inspect the region. See variant-calling/variant-normalization.

bcftools consensus -f reference.fa norm.vcf.gz 2>&1 >consensus.fa | grep -i 'overlap\|warn'

Viral / Amplicon Consensus with iVar

For amplicon surveillance (SARS-CoV-2 and similar), ivar consensus builds a per-sample consensus directly from a pileup. Its two key thresholds are epidemiological policy decisions, not defaults to accept blindly -- they propagate into lineage assignment and transmission-cluster inference:

# Trim PCR primers FIRST -- primer-derived bases are not sample sequence and, at
# primer-binding-site mutations, cause reference-biased miscalls if left in.
ivar trim -b primers.bed -p trimmed -i aligned.bam
samtools sort -o trimmed.sorted.bam trimmed.bam

# -aa keeps all positions (so no-coverage becomes N), -A keeps orphan mates, -d 0 lifts the depth cap.
samtools mpileup -aa -A -d 0 -B -Q 0 trimmed.sorted.bam | ivar consensus -p sample -q 20 -t 0.5 -m 10 -n N

| Flag | Default | Decision | |------|---------|----------| | -m min depth | 10 | Below this, iVar emits N. Too low -> single-read sequencing errors become "mutations" that corrupt outbreak phylogenies. Too high -> excessive Ns, an unusably fragmented genome. | | -t min frequency to call a base | 0 (majority) | 0 calls the most common base. For a strict majority consensus use 0.5. Too low bakes minority/within-host variants and contamination into the "genome", inflating diversity and creating phantom transmission links. Raise (e.g. 0.03) only deliberately for intrahost variant work, not for a reference consensus. | | -q min base quality | 20 | Bases below this are not counted toward depth/frequency. | | -n no-coverage char | N | Character emitted where depth < -m. |

Always report -m and -t alongside a surveillance consensus -- the genome is only as trustworthy as those two numbers. Alternatives: bcftools consensus from a called VCF, or ViralConsensus (Moshiri 2023) which calls consensus directly from the alignment without an intermediate VCF, faster and lower-memory for large batches.

Filtering Before Consensus

Apply only trusted calls; pipe filtered VCF straight into consensus:

bcftools view -f PASS input.vcf.gz -Oz -o pass.vcf.gz && bcftools index pass.vcf.gz
bcftools consensus -f reference.fa pass.vcf.gz > consensus.fa

bcftools view -v snps input.vcf.gz -Oz -o snps.vcf.gz && bcftools index snps.vcf.gz  # SNPs only

Filtered VCFs must be re-bgzipped and re-indexed before bcftools consensus reads them.

Chain Files and Naming

-c chain.txt writes a liftover chain mapping reference coordinates to consensus coordinates -- needed when indels shift positions and annotations must be lifted. -p PREFIX prepends a string to output sequence names (>sample1_chr1).

bcftools consensus -f reference.fa -c chain.txt -p "sample1_" input.vcf.gz > consensus.fa

cyvcf2 Consensus (SNP-only prototypes)

For a quick SNP-only substitution in Python (production work should use bcftools consensus, which handles indels, phasing, and masking):

from cyvcf2 import VCF
from Bio import SeqIO

ref = {rec.id: list(str(rec.seq)) for rec in SeqIO.parse('reference.fa', 'fasta')}
for v in VCF('input.vcf.gz'):
    if v.is_snp and len(v.ALT) == 1:
        ref[v.CHROM][v.POS - 1] = v.ALT[0]   # POS is 1-based; list index is 0-based
with open('consensus.fa', 'w') as fh:
    for chrom, seq in ref.items():
        fh.write(f'>{chrom}\n{"".join(seq)}\n')

Verify the Consensus

minimap2 -a reference.fa consensus.fa | samtools view -b -o aln.bam   # inspect where it diverges
bcftools view -H input.vcf.gz | wc -l                                 # variants available to apply

Common Errors

| Error / Symptom | Cause | Fix | |-----------------|-------|-----| | the VCF file is not indexed | Plain-gzip or missing index | bgzip then bcftools index (or tabix -p vcf) | | sequence "chr1" not found | Chromosome names differ between FASTA and VCF | bcftools annotate --rename-chrs map.txt | | REF does not match | Different reference than the caller used | Use the exact FASTA used for calling; normalize | | Clean haplotype looks wrong | -H 1 on an unphased VCF -> chimera | Verify | phasing; phase before -H | | Consensus reference-identical over gaps | No-coverage sites emitted as reference | Mask with samtools depth -a derived BED and -m | | Garbled indels, stderr overlap warnings | Un-normalized/overlapping records | bcftools norm -f ref.fa first; inspect warnings | | <DEL>/<INS> not applied | Symbolic SV alleles carry no ALT sequence | Use sequence-resolved SV records; see structural-variant-calling |

Related Skills

• variant-calling/variant-calling - Generate the VCF consensus is built from
• variant-calling/vcf-basics - Interpret GT and phasing (| vs /) before -H
• variant-calling/variant-normalization - Left-align indels before consensus
• variant-calling/filtering-best-practices - Restrict to trusted calls first
• variant-calling/structural-variant-calling - Sequence-resolved SVs for SV-aware consensus
• phasing-imputation/haplotype-phasing - Produce phased genotypes for true haplotypes
• phylogenetics/modern-tree-inference - Build trees from a consensus alignment

References

• Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10(2):giab008. doi:10.1093/gigascience/giab008. (bcftools consensus / norm / mpileup.)
• Grubaugh ND, Gangavarapu K, Quick J, Matteson NL, De Jesus JG, Main BJ, et al. An amplicon-based sequencing framework for accurately measuring intrahost virus diversity using PrimalSeq and iVar. Genome Biology. 2019;20(1):8. doi:10.1186/s13059-018-1618-7. (iVar consensus/trim; depth -m and frequency -t thresholds.)
• Moshiri N. ViralConsensus: a fast and memory-efficient tool for calling viral consensus genome sequences directly from read alignment data. Bioinformatics. 2023;39(5):btad317. doi:10.1093/bioinformatics/btad317.