DNA Variant Calling Pipeline with Snakemake

This project presents a reproducible bioinformatics pipeline built using Snakemake to perform DNA variant calling on Escherichia coli samples. It processes paired-end Illumina sequencing data, aligns reads to a reference genome, performs variant calling, applies quality filtering, and concludes with a comparative analysis supported by visualizations.

---

Project Structure

bio_pipeline/
├── data/                     # Input FASTQ files
├── reference/                # Reference genome and indexes
├── results/                  # Output files including BAMs, VCFs, and plots
│   └── visualizations/       # Bar chart and Venn diagram images
├── scripts/                  # Custom Python scripts
├── envs/                     # Environment configurations (optional)
├── Snakefile                 # Main Snakemake workflow
├── config.yaml               # Sample and reference configuration
└── environment.yaml          # Conda environment specification

---

Pipeline Overview

Input Data

• Two paired-end FASTQ datasets from NCBI SRA:
  • ecoli_sample1: Lab-evolved strain (SRR2584863)
  • ecoli_sample2: Field-isolated strain (SRR1212323)
• Reference genome: E. coli K-12 substr. MG1655 (GCF_000005845.2)

Workflow Steps

1. Align reads to reference using BWA
2. Convert SAM to BAM using SAMtools
3. Sort and index BAM files
4. Add read group metadata using GATK
5. Call variants using GATK HaplotypeCaller
6. Compress and index the raw VCFs with bgzip and tabix
7. Filter variants using GATK VariantFiltration with common quality thresholds
8. Compare filtered VCFs using bcftools isec
9. Generate summary plots with matplotlib and matplotlib-venn

---

Visualizations

Two plots are generated to interpret the results:

Bar Chart

Displays the number of high-confidence (PASS) variants in each sample:

Sample	Description	PASS Variants
ecoli_sample1	Lab-evolved (SRR2584863)	34,116
ecoli_sample2	Field-isolated (SRR1212323)	53,404

Venn Diagram

Shows the distribution of PASS variants:

Category	Count
Shared PASS Variants	16,090
Unique to Lab-Evolved	18,026
Unique to Field-Isolated	37,314

Visuals are saved in: results/visualizations/

---

How to Run the Pipeline

Step 1: Create and Activate the Environment

conda env create -f environment.yaml
conda activate snakemake_env

Step 2: Run the Snakemake Workflow

snakemake --cores 4

Step 3: Generate Visualizations

python scripts/visualize_variants.py

Ensure paths in config.yaml are correctly set before running the pipeline.

---

Tools and Dependencies

• Snakemake
• BWA
• SAMtools
• GATK
• bcftools
• Python 3
  • matplotlib
  • matplotlib-venn

---

Project Summary

This pipeline was applied to analyze and compare two E. coli strains: one lab-evolved (from the Lenski long-term evolution experiment) and one field-isolated sample submitted by the FDA. After processing and variant filtration, the field-isolated strain exhibited a significantly higher number of unique variants (37,314) compared to the lab-evolved strain (18,026). The two samples shared 16,090 high-confidence variants. These observations suggest a greater genomic diversity in the field isolate, potentially due to environmental exposure and adaptation.

This project demonstrates how modular, open-source tools can be chained together in a reproducible and interpretable way to explore biological questions around microbial evolution and genomic diversity.