Step 4 Variant Calling & Filtering (bcftools mpileup call) - alyrcoun/EPP622-2026-SNPCalling GitHub Wiki

Variant Calling with bcftools

Analysis Overview

Now that all four sample read sets are aligned to the S. cerevisiae reference genome, we will perform multi-sample joint variant calling. Rather than calling variants on each sample individually, giving all BAM files to bcftools together allows the caller to leverage population-level information across samples, increasing sensitivity at low-coverage loci. Following raw variant calling, we will apply quality and depth filters using bcftools filter and bcftools view to eliminate false positives caused by sequencing errors or misalignments.

Make sure you have bcftools installed. If you don't, use this code:

conda install bioconda::bcftools

Navigate to your bcftools directory.

cd bcftools_results

# Link to reference genome FASTA
ln -s ../bwamem2_results/sacCer3.fa .

# Link all four sorted BAM files and their .bai index files
ln -s ../bwamem2_results/*.sorted.bam* .

# Create bam_list.txt file so `bcftools mpileup` knows which BAM files to process together when running joint multi-sample variant calling 
ls Sample_*.sorted.bam > bam_list.txt

# Use `cat` if you want to double check it lists samples correctly
cat bam_list.txt

Now let's run the joint variant calling and filtering pipeline by piping bcftools mpileup directly into bcftools call to calculate genotype likelihoods and call raw variants without writing intermediate pileup files to disk.

# Joint multi-sample variant calling pipeline
# bcftools mpileup scans reference positions and collects alignment stacks across all BAM files
# -Ou outputs uncompressed BCF format for maximum computational efficiency when piping
# -q 5 skips reads with mapping quality < 5 (filters out multi-mapping reads)
# -d 250 sets maximum read depth per sample cap to prevent excessive memory usage
# -b bam_list.txt specifies list of input sample BAM files for joint multi-sample calling
# -f sacCer3.fa specifies reference genome FASTA needed to compare read bases against reference
# -a AD,ADF,ADR,DP,SP annotates allelic depth, forward/reverse strand depth, total depth, and strand bias
# bcftools call -m -v uses multiallelic calling model (-m) and outputs variant sites only (-v)
# -Oz -o raw_calls.vcf.gz writes bgzip-compressed VCF file
bcftools mpileup \
  -Ou \
  -q 5 \
  -d 250 \
  -b bam_list.txt \
  -f sacCer3.fa \
  -a AD,ADF,ADR,DP,SP | \
bcftools call \
  -m \
  -v \
  -Oz \
  -o raw_calls.vcf.gz

# Generate summary stats for the raw VCF calls
bcftools stats raw_calls.vcf.gz > raw_calls.stats.txt

# Filter for SNPs with Mapping Quality >= 40 and sample depth (DP) >= 10
bcftools filter -i 'TYPE="snp" & INFO/MQ>=40 & FORMAT/DP>=10' \
  -o filtered_MQ40_DP10.vcf.gz \
  raw_calls.vcf.gz

# Generate summary stats for the filtered VCF calls
bcftools stats filtered_MQ40_DP10.vcf.gz > filtered_MQ40_DP10.stats.txt

# Once those finish, run `ls -la` to ensure all 4 output files exist in your bcftools_results directory: raw_calls.vcf.gz, raw_calls.stats.txt, filtered_MQ40_DP10.vcf.gz, and filtered_MQ40_DP10.stats.txt

Results

bcftools stats report path and name: Path: /lustre/isaac24/proj/UTK0505/alyrcoun/epp_622/test2/test2_analyses/bcftools_results/ Name: raw_calls.stats.txt