This guide covers generating sample-specific sequences by applying variants to a reference. bcftools consensus substitutes ALT alleles at VCF positions and copies the reference everywhere else -- which drives the three traps to keep in mind: a consensus emits reference (not N) at no-coverage sites unless they are masked, -H 1 on an unphased VCF builds a chimeric non-haplotype, and a single FASTA cannot faithfully hold a diploid genome. For viral/amplicon surveillance, iVar builds the consensus directly from a pileup and its depth/frequency thresholds are epidemiological policy choices.
conda install -c bioconda bcftools)bcftools index input.vcf.gz)Tell your AI agent what you want to do: - "Generate a consensus FASTA by applying my VCF variants to the reference" - "Extract both haplotypes for a specific gene region from my phased VCF" - "Check whether my VCF is phased before pulling a haplotype consensus" - "Create consensus sequences for all samples in my multi-sample VCF" - "Mask no-coverage regions with N so my consensus does not report false reference" - "Build a SARS-CoV-2 consensus from amplicon reads with iVar"
bcftools consensus modifies a reference sequence by applying variants from a VCF file:
Reference: ATCGATCGATCG
VCF: pos 5: T->C
pos 8: C->G
Consensus: ATCGACCGAGCG
This creates a sequence that represents what the sample's genome looks like.
bcftools consensus -f reference.fa input.vcf.gz > consensus.fa
bcftools consensus -f reference.fa -o consensus.fa input.vcf.gz
Multi-sample VCFs require specifying which sample:
bcftools consensus -f reference.fa -s sample_name input.vcf.gz > sample.fa
bcftools query -l input.vcf.gz
At heterozygous sites, you need to choose which allele to apply:
| Option | Description |
|---|---|
-H 1 |
First haplotype (allele 1) |
-H 2 |
Second haplotype (allele 2) |
-H A |
Apply all ALT alleles |
-H R |
Keep REF at heterozygous sites |
-I |
Use IUPAC ambiguity codes |
-H 1/-H 2 are only meaningful on phased genotypes (0|1, pipe separator). On unphased genotypes (0/1, slash), which allele is "haplotype 1" is arbitrary at each site, so -H 1 across many heterozygous sites produces a chimeric sequence that matches no real chromosome. Always check first:
bcftools query -f '%CHROM\t%POS[\t%GT]\n' input.vcf.gz | head # phased 0|1, unphased 0/1
If unphased, phase first with read-backed (WhatsHap, HapCUT2), trio, statistical (SHAPEIT, Eagle), or long-read methods before extracting a haplotype.
bcftools consensus -f reference.fa -H 1 input.vcf.gz > haplotype1.fa
bcftools consensus -f reference.fa -H 2 input.vcf.gz > haplotype2.fa
For heterozygous sites, use ambiguity codes:
bcftools consensus -f reference.fa -I input.vcf.gz > consensus_iupac.fa
IUPAC codes: - R = A or G - Y = C or T - M = A or C - K = G or T - W = A or T - S = C or G
Because unobserved positions are emitted as reference, masking no-coverage sites is what stops a consensus from claiming false reference calls. Build the mask BED from callable depth. The -a in samtools depth is mandatory: without it, zero-coverage positions are omitted from the output, never enter the BED, and stay as reference:
# -a reports every position (incl. zero coverage); < 10x is a common callable-depth floor
samtools depth -a input.bam | \
awk '$3<10 {print $1"\t"$2-1"\t"$2}' | \
bedtools merge > low_coverage.bed
# Apply mask
bcftools consensus -f reference.fa -m low_coverage.bed input.vcf.gz > consensus.fa
Default mask character is N. Change it with --mask-with (NOT -M, which is --missing and only sets the character for missing ./. genotypes):
bcftools consensus -f reference.fa -m mask.bed --mask-with X input.vcf.gz > consensus.fa
bcftools consensus -f reference.fa -M N input.vcf.gz > consensus.fa
bcftools consensus -f reference.fa -r chr1 input.vcf.gz > chr1_consensus.fa
bcftools consensus -f reference.fa -r chr1:1000000-2000000 input.vcf.gz > region.fa
# Extract gene coordinates from annotation
bcftools consensus -f reference.fa -r chr1:1234567-1245678 input.vcf.gz > gene.fa
bcftools view -f PASS input.vcf.gz | \
bcftools consensus -f reference.fa > consensus.fa
bcftools filter -i 'QUAL>=30 && INFO/DP>=10' input.vcf.gz | \
bcftools consensus -f reference.fa > consensus.fa
bcftools view -v snps input.vcf.gz | \
bcftools consensus -f reference.fa > consensus_snps.fa
bcftools view -V indels input.vcf.gz | \
bcftools consensus -f reference.fa > consensus.fa
When indels change sequence length, coordinate mapping is needed:
bcftools consensus -f reference.fa -c chain.txt input.vcf.gz > consensus.fa
Chain files can be used with liftOver tools to convert coordinates: - BED coordinates - GTF/GFF annotations - Variant positions
for sample in $(bcftools query -l cohort.vcf.gz); do
echo "Processing $sample..."
bcftools consensus -f reference.fa -s "$sample" cohort.vcf.gz > "${sample}.fa"
done
sample="sample1"
bcftools consensus -f reference.fa -s "$sample" -H 1 input.vcf.gz > "${sample}_hap1.fa"
bcftools consensus -f reference.fa -s "$sample" -H 2 input.vcf.gz > "${sample}_hap2.fa"
Output uses reference chromosome names.
bcftools consensus -f reference.fa -p "sample1_" input.vcf.gz > sample1.fa
Results in: >sample1_chr1, >sample1_chr2, etc.
For amplicon surveillance (SARS-CoV-2 and similar), ivar consensus builds a consensus directly from a pileup rather than from a VCF. Trim PCR primers first, then pipe samtools mpileup into ivar consensus:
ivar trim -b primers.bed -p trimmed -i aligned.bam
samtools sort -o trimmed.sorted.bam trimmed.bam
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
The two thresholds are epidemiological policy decisions, not defaults to accept blindly:
- -m (min depth, default 10): below this, iVar emits N. Too low turns single-read errors into "mutations"; too high produces excessive Ns.
- -t (min frequency, default 0 = majority base): too low bakes minority/contaminant variants into the genome and corrupts transmission inference. Use 0.5 for a strict majority consensus; lower values only for deliberate intrahost work.
Report both numbers alongside any surveillance consensus.
bcftools consensus cannot apply symbolic SV alleles (<DEL>, <INS>, <DUP>), because they carry no ALT sequence 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) or an assembly-based approach.
Generate consensus for all samples for multiple sequence alignment:
mkdir -p consensus
# Reference as first sequence
cp reference.fa consensus/reference.fa
# Each sample
for sample in $(bcftools query -l variants.vcf.gz); do
bcftools view -s "$sample" variants.vcf.gz | \
bcftools filter -i 'QUAL>=30' | \
bcftools consensus -f reference.fa > "consensus/${sample}.fa"
done
# Combine all sequences
cat consensus/*.fa > all_samples.fa
# Create depth mask
samtools depth -a aligned.bam | \
awk '$3<20 {print $1"\t"$2-1"\t"$2}' | \
bedtools merge > low_depth.bed
# Generate consensus
bcftools filter -i 'QUAL>=30 && INFO/DP>=20' variants.vcf.gz | \
bcftools view -f PASS | \
bcftools consensus -f reference.fa -m low_depth.bed > consensus.fa
# Report coverage
masked=$(awk '{sum += $3-$2} END {print sum}' low_depth.bed)
total=$(samtools faidx reference.fa -n 1 | awk 'NR==2 {print length}')
echo "Masked $masked of $total bases ($(echo "scale=1; $masked*100/$total" | bc)%)"
gene_region="chr1:1000000-1010000"
sample="patient_001"
# Haplotype 1
bcftools consensus -f reference.fa -r "$gene_region" -s "$sample" -H 1 \
variants.vcf.gz > gene_hap1.fa
# Haplotype 2
bcftools consensus -f reference.fa -r "$gene_region" -s "$sample" -H 2 \
variants.vcf.gz > gene_hap2.fa
# Variants in VCF
bcftools view -H input.vcf.gz | wc -l
# Check for warnings during consensus
bcftools consensus -f reference.fa input.vcf.gz 2>&1 | grep -c "warning"
# Simple character comparison
diff <(grep -v "^>" reference.fa | tr -d '\n') \
<(grep -v "^>" consensus.fa | tr -d '\n') | head
VCF file needs an index:
bcftools index input.vcf.gz
Chromosome names don't match between reference and VCF:
# Check reference chromosomes
grep "^>" reference.fa | head
# Check VCF chromosomes
bcftools view -H input.vcf.gz | cut -f1 | sort -u | head
Fix by renaming:
bcftools annotate --rename-chrs rename.txt input.vcf.gz -Oz -o renamed.vcf.gz
VCF was called with a different reference. The reference FASTA must exactly match what was used for variant calling.
Variants at the same position or overlapping indels. This is usually handled automatically but check the output:
bcftools consensus -f reference.fa input.vcf.gz 2>&1 | grep -i "overlap"
"Generate a consensus FASTA sequence by applying my VCF variants to the reference"
"Extract both haplotypes for a specific gene region from my phased VCF"
"Create consensus sequences for all samples in my multi-sample VCF"
"Mask low-coverage regions in my consensus sequence with N characters"
"Check whether my genotypes are phased before I extract a haplotype consensus"
"Build a SARS-CoV-2 consensus genome from my amplicon BAM with iVar at 10x minimum depth"
bcftools norm before consensus| phasing before -H 1/-H 2; on unphased genotypes it yields a chimeric non-haplotypesamtools depth -a derived BED -- -a is mandatory or zero-depth positions escape masking and stay as reference-I) only when downstream tools support them; many tree builders read IUPAC as Nbcftools consensus cannot apply symbolic SV alleles (<DEL>/<INS>/<DUP>) -- use sequence-resolved records-m (min depth) and -t (min frequency) thresholds -- they are policy choices that affect lineage and transmission inference-c chain.txt to produce a chain file when indels shift coordinates