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:

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:

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.