Generate consensus FASTA sequences by applying VCF variants to a reference using bcftools consensus. Use when creating sample-specific reference sequences or reconstructing haplotypes.
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.
Consensus Sequences
"Generate a consensus sequence from my VCF" → Apply called variants to a reference FASTA, producing a sample-specific genome with optional haplotype selection and low-coverage masking.
Goal: Generate consensus sequences for downstream analyses like phylogenetics, viral surveillance, or gene-level comparison.
Approach: Filter variants to high-quality calls, apply per-sample consensus generation, mask low-coverage regions with N, then combine for multi-sample workflows.
Phylogenetic Analysis Preparation
# For each sample, generate consensus
mkdir -p consensus
for sample in $(bcftools query -l cohort.vcf.gz); do
bcftools view -s "$sample" cohort.vcf.gz | \
bcftools view -c 1 | \
bcftools consensus -f reference.fa > "consensus/${sample}.fa"
done
# Combine for alignment
cat consensus/*.fa > all_samples.fa
from cyvcf2 import VCF
from Bio import SeqIO
# Load reference
ref_dict = {rec.id: str(rec.seq) for rec in SeqIO.parse('reference.fa', 'fasta')}
# Apply variants (SNPs only, simplified)
vcf = VCF('input.vcf.gz')
changes = {}
for variant in vcf:
if variant.is_snp and len(variant.ALT) == 1:
chrom = variant.CHROM
pos = variant.POS - 1 # 0-based
if chrom not in changes:
changes[chrom] = {}
changes[chrom][pos] = variant.ALT[0]
# Apply changes
for chrom, positions in changes.items():
seq = list(ref_dict[chrom])
for pos, alt in positions.items():
seq[pos] = alt
ref_dict[chrom] = ''.join(seq)
# Write output
with open('consensus.fa', 'w') as f:
for chrom, seq in ref_dict.items():
f.write(f'>{chrom}\n{seq}\n')
Note: Use bcftools consensus for production - handles indels and edge cases properly.
