🧬 Variant Calling Workflow using BWA, SAMtools, and BCFtools

This guide outlines the steps to identify genetic variants in E. coli using publicly available tools and data. It is meant for educational use, particularly for school students learning about bioinformatics pipelines.

---

📌 Workflow Overview

---

🛠️ Step 1: Install Required Tools

sudo apt install bwa samtools bcftools

Tools installed:

• bwa: For aligning sequencing reads.
• samtools: For handling SAM/BAM files.
• bcftools: For variant calling.

---

📂 Step 2: Download Reference Genome

mkdir -p /home/data/ref_genome

curl -L -o ecoli_rel606.fasta.gz \
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/017/985/GCA_000017985.1_ASM1798v1/GCA_000017985.1_ASM1798v1_genomic.fna.gz

gunzip ecoli_rel606.fasta.gz

---

📦 Step 3: Download Sample Read Data

curl -L -o sub.tar.gz https://ndownloader.figshare.com/files/14418248
tar xvf sub.tar.gz
mv sub/ ~/

This dataset contains paired-end FASTQ files.

---

🧬 Step 4: Index the Reference Genome

bwa index ecoli_rel606.fasta

Indexing prepares the reference for alignment.

---

🧾 Step 5: Align Reads with BWA

bwa mem ecoli_rel606.fasta \
SRR2584866_1.trim.sub.fastq \
SRR2584866_2.trim.sub.fastq \
> SRR2584866.aligned.sam

Output: SRR2584866.aligned.sam

---

🔄 Step 6: Convert SAM to BAM

samtools view -S -b SRR2584866.aligned.sam > SRR2584866.aligned.bam

BAM files are binary and compressed.

---

📚 Step 7: Sort and index the BAM File

samtools sort -o SRR2584866.aligned.sorted.bam SRR2584866.aligned.bam

 samtools index SRR2584866.aligned.sorted.bam

Sorting and indexing improve processing speed for the next steps.

---

📊 Step 8: Check BAM File Stats

samtools flagstat SRR2584866.aligned.sorted.bam

This gives summary metrics such as total reads and mapping rates.

---

🔍 Step 9: Generate Raw BCF File

bcftools mpileup -O b -o SRR2584866_raw.bcf \
-f ecoli_rel606.fasta \
SRR2584866.aligned.sorted.bam

---

🧬 Step 10: Call Variants

bcftools call --ploidy 1 -m -v -o SRR2584866_variants.vcf SRR2584866_raw.bcf

• --ploidy 1: Because bacteria are haploid.
• Output: VCF file listing SNPs/indels.

---

👁️ Step 11: View Variants

less -S SRR2584866_variants.vcf

Use arrow keys to scroll horizontally if needed.

---

✅ Conclusion

This workflow covers a full variant calling pipeline starting from raw reads through to variant identification. It introduces essential tools and formats (FASTQ, FASTA, SAM, BAM, BCF, VCF) in genome analysis.

---