bio-consensus-sequences

Version Compatibility

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

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.

• CLI: bcftools consensus -f reference.fa input.vcf.gz
• Python: cyvcf2 + Bio.SeqIO for simple SNP-only cases

Basic Usage

Generate Consensus

bcftools consensus -f reference.fa input.vcf.gz > consensus.fa

Specify Sample

bcftools consensus -f reference.fa -s sample1 input.vcf.gz > sample1.fa

Output to File

bcftools consensus -f reference.fa -o consensus.fa input.vcf.gz

Haplotype Selection

First Haplotype Only

bcftools consensus -f reference.fa -H 1 input.vcf.gz > haplotype1.fa

Second Haplotype Only

bcftools consensus -f reference.fa -H 2 input.vcf.gz > haplotype2.fa

Haplotype Options

Option	Description
-H 1	First haplotype
-H 2	Second haplotype
-H A	Apply all ALT alleles
-H R	Apply REF alleles where heterozygous
-I	Apply IUPAC ambiguity codes (separate flag)

IUPAC Codes for Heterozygous Sites

bcftools consensus -f reference.fa -I input.vcf.gz > consensus_iupac.fa

Heterozygous sites encoded with IUPAC ambiguity codes:

• A/G → R
• C/T → Y
• A/C → M
• G/T → K
• A/T → W
• C/G → S

Missing Data Handling

Mark Missing as N

bcftools consensus -f reference.fa -M N input.vcf.gz > consensus.fa

Mark Low Coverage as N

Using a mask BED file:

# Create mask from depth
samtools depth input.bam | awk '$3<10 {print $1"\t"$2-1"\t"$2}' > low_coverage.bed

# Apply mask
bcftools consensus -f reference.fa -m low_coverage.bed input.vcf.gz > consensus.fa

Mask Options

Option	Description
-m FILE	Mask regions in BED file with N
-M CHAR	Character for masked regions (default N)

Region Selection

Specific Region

bcftools consensus -f reference.fa -r chr1:1000-2000 input.vcf.gz > region.fa

Multiple Regions

Use with BED file to extract multiple regions.

Chain Files

Generate Chain File

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

Chain files map coordinates between reference and consensus:

• Useful for liftover of annotations
• Required when indels change sequence length

Chain File Format

chain score ref_name ref_size ref_strand ref_start ref_end query_name query_size query_strand query_start query_end id

Sample-Specific Consensus

For Each Sample

for sample in $(bcftools query -l input.vcf.gz); do
    bcftools consensus -f reference.fa -s "$sample" input.vcf.gz > "${sample}.fa"
done

Both Haplotypes

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"

Filtering Before Consensus

PASS Variants Only

bcftools view -f PASS input.vcf.gz | \
    bcftools consensus -f reference.fa > consensus.fa

High-Quality Variants Only

bcftools filter -i 'QUAL>=30 && INFO/DP>=10' input.vcf.gz | \
    bcftools consensus -f reference.fa > consensus.fa

SNPs Only

bcftools view -v snps input.vcf.gz | \
    bcftools consensus -f reference.fa > consensus_snps.fa

Sequence Naming

Default Naming

Output uses reference sequence names.

Custom Prefix

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

Sequences named: sample1_chr1, sample1_chr2, etc.

Common Workflows

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

Viral Genome Assembly

# Apply high-quality variants only
bcftools filter -i 'QUAL>=30 && INFO/DP>=20' variants.vcf.gz | \
    bcftools view -f PASS | \
    bcftools consensus -f reference.fa -M N > consensus.fa

Gene-Specific Consensus

# Extract gene region
bcftools consensus -f reference.fa -r chr1:1000000-1010000 \
    -s sample1 variants.vcf.gz > gene.fa

Masked Low-Coverage Regions

# Create mask from coverage
samtools depth -a input.bam | \
    awk '$3<5 {print $1"\t"$2-1"\t"$2}' | \
    bedtools merge > low_coverage.bed

# Generate consensus with mask
bcftools consensus -f reference.fa -m low_coverage.bed \
    variants.vcf.gz > consensus.fa

Verify Consensus

Check Differences

# Align consensus to reference
minimap2 -a reference.fa consensus.fa | samtools view -bS > alignment.bam

# Or simple comparison
diff <(grep -v "^>" reference.fa) <(grep -v "^>" consensus.fa) | head

Count Changes

# Number of differences
bcftools view -H input.vcf.gz | wc -l

Handling Overlapping Variants

bcftools consensus handles overlapping variants automatically:

• Applies variants in order
• Warns about conflicts

Check for warnings:

bcftools consensus -f reference.fa input.vcf.gz 2>&1 | grep -i warn

cyvcf2 Consensus (Simple Cases)

Manual Consensus Generation

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.

Quick Reference

Task	Command
Basic consensus	bcftools consensus -f ref.fa in.vcf.gz
Specific sample	bcftools consensus -f ref.fa -s sample in.vcf.gz
Haplotype 1	bcftools consensus -f ref.fa -H 1 in.vcf.gz
IUPAC codes	bcftools consensus -f ref.fa -I in.vcf.gz
With mask	bcftools consensus -f ref.fa -m mask.bed in.vcf.gz
Generate chain	bcftools consensus -f ref.fa -c chain.txt in.vcf.gz
Specific region	bcftools consensus -f ref.fa -r chr1:1-1000 in.vcf.gz

Common Errors

Error	Cause	Solution
not indexed	VCF not indexed	Run bcftools index
sequence not found	Chromosome mismatch	Check chromosome names
overlapping records	Variants overlap	Usually OK, check warnings
REF does not match	Wrong reference	Use same reference as caller

Related Skills

• variant-calling - Generate VCF for consensus
• filtering-best-practices - Filter variants before consensus
• variant-normalization - Normalize indels first
• alignment-files/reference-operations - Reference manipulation