Module 2.4 bcftools and vcf filtering - mestato/EPP622_2026 GitHub Wiki
Call variants on one sample
We're going to start with a slow, step by step example of one file so that we can see what is happening. Lets set up.
Get bcftools installed.
conda install bioconda::bcftools
From your analyses directory.
mkdir 04_bcftools
cd 04_bcftools
Link to a single bam file and its index file. (You may need to change name below to get one of your files). Link to the reference.
ln -s ../03_bwamem2/mapped/21_FP_OH_FSB_F_10.sorted.bam .
ln -s ../03_bwamem2/mapped/21_FP_OH_FSB_F_10.sorted.bam.bai .
We start with mpileup. It walks across every position in the reference, collects the bases stacked there from read alignments, and computes the likelihood of each possible genotype.
bcftools mpileup \
-f /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.fasta \
-a AD,ADF,ADR,DP,SP \
-q 5 \
-d 250 \
-Ou \
-o 21_FP_OH_FSB_F_10.pileup.bcf \
21_FP_OH_FSB_F_10.sorted.bam
Flag decoder:
| Flag | Meaning |
|---|---|
| -f | reference FASTA (indexed) |
| -a | extra annotations to carry forward: allelic depth, strand-split AD, total depth, strand bias |
| -q 5 | skip reads with mapping quality < 5 (drops multi-mappers) |
| -d 250 | max depth per sample; caps runtime in pileup-y regions |
| -Ou | uncompressed BCF — fastest when you're about to pipe (we aren't doing that here but will be soon) |
We can take a look at what mpileup is doing
bcftools index 21_FP_OH_FSB_F_10.pileup.bcf
bcftools view -r Chr01:10000-200000 21_FP_OH_FSB_F_10.pileup.bcf | tail
Now we turn the likelihoods into actual genotypes.
bcftools call \
-m \
-v \
-Ob \
-o 21_FP_OH_FSB_F_10.raw.bcf \
21_FP_OH_FSB_F_10.pileup.bcf
Lets take a look at some genotypes
bcftools index 21_FP_OH_FSB_F_10.raw.bcf
bcftools view -r Chr01:10000-2000000 21_FP_OH_FSB_F_10.raw.bcf | tail
Normally we do not save the mpileup results. Best practice is to pipe the mpileup straight to the bcf file.
How many variants did we get?
bcftools stats 21_FP_OH_FSB_F_10.raw.bcf > 21_FP_OH_FSB_F_10.raw.stats.txt
OPTIONAL - at this point you could use bcftools norm. Indels get left-aligned and multiallelic records get split so one line = one ALT allele. If you aren't throwing away indels and multiallelic records, this is probably a good idea.
OPTIONAL - at this point you could use bcftools +fill-tags. You may wish to add the INFO fields you'll filter on later on. For example the caller does not write AF, MAF, or Hardy-Weinberg stats.
Its typical to do filtering at this point. The filtering you do will be based on your project and your experimental design.
bcftools filter -i 'TYPE="snp" & INFO/MQ>=40 & FORMAT/DP>=10' -o 21_FP_OH_FSB_F_10.filteredMQ40DP10.bcf 21_FP_OH_FSB_F_10.raw.bcf
bcftools stats 21_FP_OH_FSB_F_10.filteredMQ40DP10.bcf > 21_FP_OH_FSB_F_10.filteredMQ40DP10.stats.txt
Remember INFO is SNP level, FORMAT is genotype level. so this is SNPs only, with a SNP quality of 40 or more and sample read depth of 10 or more.
| Filter | Meaning |
|---|---|
| QUAL | Variant Quality Score = The probability that a variant exists at that site. Low value means the evidence for a mutation is weak (often due to low sequencing depth or sequencing errors). |
| MQ | Mapping Quality = The confidence that the reads belong to that specific genomic location. Low value means the reads are mapping to multiple places in the genome (e.g., repeats, duplications). |
| INFO/DP | Site level depth. Total raw or combined read depth across all samples combined at this specific genomic coordinate. |
| FORMAT/DP | Individual dpeth. Read depth for an individual sample at this specific coordinate. |
Bcftools filter is a bit of a nightmare so pay close attention to what you are doing.
- -i includes, -e excludes
- Inside FORMAT expressions, | and & are per-sample operators, while || and && collapse across all samples.
- you actually use bcftools view (not filter) to get biallelic snps
bcftools view --types snps -m2 -M2 -Oz -o b21_FP_OH_FSB_F_10.filteredMQ40DP10Biallelic.vcf.gz 21_FP_OH_FSB_F_10.filteredMQ40DP10.bcf
bcftools stats b21_FP_OH_FSB_F_10.filteredMQ40DP10Biallelic.vcf.gz > b21_FP_OH_FSB_F_10.filteredMQ40DP10Biallelic.stats.txt
| Flag | Meaning |
|---|---|
| --types snps | Keeps only Single Nucleotide Polymorphisms, removing indels and other variant types. |
| -m2 (or --min-alleles 2) | Sets the minimum number of alleles to 2 (removes monomorphic/reference sites). |
| -M2 (or --max-alleles 2) | Sets the maximum number of alleles to 2 (removes multiallelic sites with >3 alleles) |
| -Oz | Outputs compressed VCF format (instead of bcf) |
Next, we are going to call variants on multiple samples. We aren't going to use a loop (surprise!). We are going to give all our bam files to bcftools as a list (why would bcftools want to look at all our bamfiles together???)
We will also pipe straight from mpileup to call, because writing mpileup files is unecessary and they take up a lot of room.
Make file of bams
ls ../03_bwamem2/*bam > bam_list.txt
Run snp calling
bcftools mpileup \
-Ou \
-q 5 \
-b bam_list.txt \
-f /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.fasta \
-a AD,ADF,ADR,DP,SP | \
bcftools call \
-m \
-v \
-Oz \
-o module2_data.raw.vcf.gz
bcftools stats module2_data.raw.vcf.gz > module2_data.raw.stats.txt
bcftools filter -i 'TYPE="snp" & INFO/MQ>=40 & FORMAT/DP>=10' -o module2_data.filteredMQ40DP10.vcf.gz module2_data.raw.vcf.gz
bcftools stats module2_data.filteredMQ40DP10.vcf.gz > module2_data.filteredMQ40DP10.stats.txt
done