Skip to content

4. Raw Variant Calling (bcftools)

Jump to bottom
AreliDGS edited this page Oct 9, 2026 · 5 revisions

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.

Clone this wiki locally