4. Raw Variant Calling (bcftools)

This section documents the execution loops for filtering strategies utilized to identify SNPs and Indels across all samples against the Saccharomyces cerevisiae reference genome panel.

---

Workspace and coding

# 1. Make sure you are in the project root directory.
cd /lustre/isaac24/proj/UTK0505/agonza84/EPP622-2026-SNPCalling

# 2. Create and open the variant calling pipeline script
vi scripts/04_variant_calling.sh

# 3. Give execution permissions 
chmod +x scripts/04_variant_calling.sh

---

Execution script code

The script04_variant_calling.sh, performs alignment stacking, joint genotyping, and progressive quality filtrations.

#!/bin/bash

# Exit immediately if any command or pipe fails
set -euo pipefail

# Initialize output subdirectory structure
mkdir -p results/04_variant_calling

# === Phase 1: Generate Explicit BAM Input List ===
echo "results/03_alignment/Sample_A.sorted.bam" > results/04_variant_calling/bam_list.txt
echo "results/03_alignment/Sample_B.sorted.bam" >> results/04_variant_calling/bam_list.txt
echo "results/03_alignment/Sample_C.sorted.bam" >> results/04_variant_calling/bam_list.txt
echo "results/03_alignment/Sample_D.sorted.bam" >> results/04_variant_calling/bam_list.txt

# === Phase 2: Joint Variant Calling (Piping mpileup straight to call) ===
bcftools mpileup \
  -Ou \
  -q 5 \
  -b results/04_variant_calling/bam_list.txt \
  -f data/reference/s_cerevisiae_ref.fasta \
  -a AD,ADF,ADR,DP,SP | \
bcftools call \
  -m \
  -v \
  -Oz \
  -o results/04_variant_calling/yeast_data.raw.vcf.gz

# Extract initial diagnostics baseline
bcftools stats results/04_variant_calling/yeast_data.raw.vcf.gz > results/04_variant_calling/yeast_data.raw.stats.txt

# === Phase 3: Site-level and Genotype Quality Filtering ===
bcftools filter \
  -i 'TYPE="snp" & INFO/MQ>=40 & FORMAT/DP>=10' \
  -Oz \
  -o results/04_variant_calling/yeast_data.filteredMQ40DP10.vcf.gz \
  results/04_variant_calling/yeast_data.raw.vcf.gz

# Extract post-filter diagnostics
bcftools stats results/04_variant_calling/yeast_data.filteredMQ40DP10.vcf.gz > results/04_variant_calling/yeast_data.filteredMQ40DP10.stats.txt

# === Phase 4: Isolation of Strict Biallelic SNPs ===
bcftools view \
  --types snps \
  -m2 \
  -M2 \
  -Oz \
  -o results/04_variant_calling/yeast_data.filteredMQ40DP10Biallelic.vcf.gz \
  results/04_variant_calling/yeast_data.filteredMQ40DP10.vcf.gz

# Extract final biological dataset diagnostics
bcftools stats results/04_variant_calling/yeast_data.filteredMQ40DP10Biallelic.vcf.gz > results/04_variant_calling/yeast_data.filteredMQ40DP10Biallelic.stats.txt

---

Code explanation

bcftools mpileup was used adding the options:

• -f: Indicate the indexed host FASTA genome structure.
• -a AD,ADF,ADR,DP,SP: Extraction of filtering tags (Allelic Depth, Strand-split depths, total raw site depth, and strand bias calculations).
• -q 5: Minimum mapping quality filter threshold.
• -Ou: Directs uncompressed BCF file.

bcftools view was used adding the options:

• --types snps: Strips out insertion/deletion mutation and preserving single nucleotide variations.
• -m2 / -M2: Locks allele states strictly to 2.

---

Usage: 04_variant_calling.sh

# Run the script in the background using nohup
nohup ./scripts/04_variant_calling.sh > results/04_variant_calling.log 2>&1 &

# See the execution progress in real time
tail -f results/04_variant_calling.log

---

Results

• bcftools stats raw report path and name: results/04_variant_calling/yeast_data.raw.stats.txt
• bcftools stats filtered report path and name: results/04_variant_calling/yeast_data.filteredMQ40DP10.stats.txt
• bcftools stats biallelic report path and name: results/04_variant_calling/yeast_data.filteredMQ40DP10Biallelic.stats.txt

The progress tracking table below isolates the number of valid variant sites surviving consecutive bioinformatic stringency filters:

Analysis Stage	Source Target VCF File	Identified SNPs Count	Identified Indels Count
1. Raw Discovery	yeast_data.raw.vcf.gz	4,920	0
2. Quality Filtered	yeast_data.filteredMQ40DP10.vcf.gz	4,802	0
3. Biallelic Selection	yeast_data.filteredMQ40DP10Biallelic.vcf.gz	4,802	0

Findings:

• Indes: Zero indels
• Filter: Applying the threshold of mapping quality (MQ>=40) and sample depth (DP>=10) effectively eliminated 118 low-confidence or false-positive variant positions, just leave the good ones.
• Biallelic integrity: The secondary filter targeting biallelic state constraints (-m2 -M2) showed no change in total SNP counts.