Step 4: Variant Calling (bcftools mpileup call) - aduckwo4-boop/EPP622-2026-SNPCalling GitHub Wiki

Step 4: Variant Calling (bcftools mpileup + call)

Commands

BCFtools was used to identify genetic variants across all four samples by examining the aligned sequencing reads. The mpileup command summarizes the read evidence at each genomic position, while call identifies variants relative to the reference genome. Variants were then filtered based on SNP type, mapping quality, and sequencing depth to retain higher-confidence calls.

Created the variant calling directory and BAM file list:

mkdir 04_bcftools
cd 04_bcftools
ln -s ../03_bwamem2/*.sorted.bam > bam_list.txt
nano run_bcftools.sh

# ln -s: creates symbolic links to the BAM files.
# ls: lists the sorted BAM files and saves their paths in bam_list.txt.
# >: redirects the file list into a text file.

Created run_bcftools.sh:

bcftools mpileup \
  -Ou \
  -q 5 \
  -b bam_list.txt \
  -f ../ref/sacCer3.fa \
  -a AD,ADF,ADR,DP,SP | \
bcftools call \
  -m \
  -v \
  -Oz \
  -o all_samples.raw.vcf.gz

bcftools stats all_samples.raw.vcf.gz > all_samples.raw.stats.txt
bcftools filter -i 'TYPE="snp" & INFO/MQ>=40 & FORMAT/DP>=10' -o all_samples.filtered.vcf.gz all_samples.raw.vcf.gz
bcftools stats all_samples.filteredMQ40DP10.vcf.gz > all_samples.filtered.stats.txt

# bcftools mpileup: calculates genotype likelihoods from aligned sequencing reads.
# -Ou: outputs uncompressed BCF data for piping into bcftools call.
# -q 5: excludes reads with mapping quality below 5.
# -b bam_list.txt: specifies the BAM files used for joint variant calling.
# -f: specifies the reference genome FASTA file.
# -a: requests additional read-level annotations.
# AD: allele depths for reference and alternate alleles.
# ADF: allele depths from forward-strand reads.
# ADR: allele depths from reverse-strand reads.
# DP: sequencing depth.
# SP: strand bias P-value.
# |: pipes the mpileup output directly into bcftools call.
# bcftools call: identifies genetic variants from genotype likelihoods.
# -m: uses the multiallelic variant calling model.
# -v: outputs variant sites only.
# -Oz: outputs a compressed VCF file.
# -o: specifies the output filename.
# bcftools stats: generates summary statistics for the VCF file.
# >: saves the statistics to a text file.
# bcftools filter: filters variants based on specified criteria.
# -i: includes only variants meeting all filtering conditions.
# TYPE="snp": retains only single nucleotide polymorphisms.
# &: logical AND, requiring all conditions to be satisfied.
# INFO/MQ>=40: retains variants with site-level mapping quality of at least 40.
# FORMAT/DP>=10: requires a sample-level sequencing depth of at least 10
# for at least one sample at the site.

Linked the reference directory and ran the script:

ln -s ../ref
bash run_bcftools.sh

Results

bcftools stats report path and name: '/lustre/isaac24/proj/UTK0505/aduckwo4/622_test2/04_bcftools/all_samples.raw.stats.txt`