domains / science / vcf-consensus-fasta / consensus-sequences

consensus-sequences

Skill card

combined

Both copies describe the same core task (applying VCF variants to a reference with bcftools consensus), but the GPTomics copy adds material constraints and workflows the FreedomIntelligence copy lacks: a governing principle about silent reference emission at zero-coverage sites, the chimeric-haplotype footgun on unphased VCFs, version-compatibility checks, and an iVar-based viral/amplicon consensus path. The FreedomIntelligence copy adds specific chain-file generation, custom sequence prefixing, and a broader set of filtering examples. Crowning one would hide the other's real alternative steps and safety constraints.

Repositories

Other copies

Combined skill

# Consensus Sequences

Generate consensus FASTA sequences by applying VCF variants to a reference using bcftools consensus, or build viral/amplicon consensus with iVar.

## 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, introspect the installed package or CLI help and adapt the example to match the actual API rather than retrying. Note: the -H argument vocabulary 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.

## 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. This drives three traps:

1. A consensus silently emits REFERENCE wherever the VCF is silent, including positions with zero coverage. Mask no-coverage sites 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. Verify | phasing before selecting a haplotype.
3. A single FASTA cannot faithfully represent a diploid genome. Every projection is lossy; for phase-sensitive work keep the VCF, not the consensus.

The input VCF must be bgzipped and indexed (bgzip + bcftools index/tabix). The REF bases in the VCF must match the FASTA exactly or bcftools warns and skips those records. Normalize first if needed.

## Basic Usage

```bash
bcftools index input.vcf.gz
bcftools consensus -f reference.fa input.vcf.gz > consensus.fa
bcftools consensus -f reference.fa -o consensus.fa input.vcf.gz
```

For a multi-sample VCF, always pass -s:

```bash
bcftools query -l input.vcf.gz
bcftools consensus -f reference.fa -s sample1 input.vcf.gz > sample1.fa
```

Restrict to a region with -r:

```bash
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 -I / --iupac-codes) | IUPAC ambiguity code | Retain heterozygosity in one sequence |
| -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. Verify phasing before any -H 1/-H 2:

```bash
bcftools query -f '%CHROM\t%POS[\t%GT]\n' input.vcf.gz | head
```

If genotypes are unphased, phase first (read-backed WhatsHap/HapCUT2, trio, statistical SHAPEIT/Eagle, or native long-read phasing).

## IUPAC Codes for Heterozygous Sites

```bash
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 and Masking

Mark missing as N:

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

Mark low coverage as N using a mask BED file:

```bash
samtools depth input.bam | awk '$3<10 {print $1"\t"$2-1"\t"$2}' > low_coverage.bed
bcftools consensus -f reference.fa -m low_coverage.bed input.vcf.gz > consensus.fa
```

Mask options: -m FILE masks regions in a BED file with N; -M CHAR sets the character for masked regions (default N).

## Region Selection

Specific region:

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

Multiple regions: use with a BED file to extract multiple regions.

## Chain Files

Generate a chain file to map coordinates between reference and consensus:

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

Chain files are useful for liftover of annotations and required when indels change sequence length.

## Sample-Specific Consensus

For each sample:

```bash
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 (only if phased):

```bash
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:

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

High-quality variants only:

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

SNPs only:

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

## Sequence Naming

Default naming uses reference sequence names. Custom prefix:

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

Sequences named: sample1_chr1, sample1_chr2, etc.

## Common Workflows

Phylogenetic analysis preparation:

```bash
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
cat consensus/*.fa > all_samples.fa
```

Viral genome assembly (bcftools path):

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

Viral/amplicon consensus (iVar path):

```bash
samtools mpileup -f reference.fa input.bam | ivar consensus -p out
```

Gene-specific consensus:

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

Masked low-coverage regions:

```bash
samtools depth -a input.bam | \
    awk '$3<5 {print $1"\t"$2-1"\t"$2}' | \
    bedtools merge > low_coverage.bed
bcftools consensus -f reference.fa -m low_coverage.bed variants.vcf.gz > consensus.fa
```

## Verify Consensus

Check differences:

```bash
minimap2 -a reference.fa consensus.fa | samtools view -bS > alignment.bam
diff <(grep -v "^>" reference.fa) <(grep -v "^>" consensus.fa) | head
```

Count changes:

```bash
bcftools view -H input.vcf.gz | wc -l
```

## Handling Overlapping Variants

bcftools consensus handles overlapping variants automatically, applying variants in order and warning about conflicts. Check for warnings:

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