Nextflow SNP calling workflow

May 22, 2026 · 3 min read
projects

nf_snp_calling is an automated, scalable Nextflow pipeline engineered for high-throughput short-read variant calling, joint regenotyping, and cohort-scale population genomics.

While standard GATK Best Practices offer sensitive joint genotyping across large cohorts, their Java-based architectures often incur heavy compute runtimes and high memory overhead. To solve this, nf_snp_calling implements an end-to-end workflow centered around BWA-MEM2, SAMtools, and BCFtools—mimicking the two-pass rigor of GATK joint genotyping while maintaining the speed and low resource footprint of native C/C++ utilities.

flowchart LR subgraph " " subgraph params v9["meta"] v0["ref_fasta"] v4["cohort_id"] end v3{ } v7([BWAMEM2_INDEX]) v11([FASTP]) v12([BWAMEM2_MEM]) v13([SAMTOOLS_FIXMATE]) v14([SAMTOOLS_SORT]) v15([SAMTOOLS_MARKDUP]) v17([SNPCALLING]) v19([BCFTOOLS_MERGE]) v21([FILTER]) v23([CHUNK_VCF]) v27([JOINT_REGENOTYPE]) v29([BCFTOOLS_CONCAT]) v31([FILTER_FINAL]) v33([BCFTOOLS_STATS]) subgraph s1[" "] v4["cohort_id"] v5([SAMTOOLS_FAIDX]) end v0 --> v3 v0 --> v5 v4 --> v5 v0 --> v7 v4 --> v7 v9 --> v11 v0 --> v12 v7 --> v12 v11 --> v12 v12 --> v13 v0 --> v14 v4 --> v14 v13 --> v14 v0 --> v15 v4 --> v15 v14 --> v15 v0 --> v17 v4 --> v17 v15 --> v17 v0 --> v19 v17 --> v19 v4 --> v19 v19 --> v21 v21 --> v23 v0 --> v27 v4 --> v27 v23 --> v27 v15 --> v27 v27 --> v29 v29 --> v31 v31 --> v33 v3 --> s1 end

Steps:

  1. Preprocessing & Alignment: Raw paired-end reads undergo quality filtering and adapter trimming via fastp, followed by indexing and alignment to the reference genome with bwa-mem2.
  2. Alignment Sanitization: Mappings are processed through samtools fixmate (retaining properly paired reads with -m), coordinate-sorted via samtools sort, and marked for duplicates using samtools markdup.
  3. Pass 1 — Per-Sample Discovery & Cohort Merge: Individual candidate variant sites are called across each sample with bcftools mpileup | bcftools call (haploid/custom ploidy supported). The resulting sample VCFs are merged across the cohort (bcftools merge) and filtered (QUAL >= 30, AC > 0) to build an initial candidate site catalog.
  4. VCF Chunking & Parallel Regenotyping: The merged candidate VCF is split into manageable genomic chunks (e.g., 50 kb intervals). All cohort BAMs are then re-evaluated in parallel against each chunk using bcftools mpileup with target site forcing and bcftools call -C alleles --insert-missed to recover marginal coverage, invariant reference alleles, and missing genotype states across every individual.
  5. Concatenation & Final Quality Filtering: The regenotyped chunks are assembled with bcftools concat --naive and subjected to final strict filtering (bcftools view -m2 -M2 -v snps -i 'AC>0 && QUAL>=30') to yield a high-confidence, biallelic SNP matrix.
  6. Downstream Population Analytics: The pipeline generates comprehensive summary metrics via bcftools stats and produces indexed VCFs ready for direct integration into linkage disequilibrium (LD) decay assays, population structure profiling (PCA/Admixture), and demographic modeling.

Highlights

  • Two-Pass Joint Genotyping Logic: Combines cohort-wide candidate discovery with chunked, multi-sample regenotyping (--insert-missed) to resolve reference calls and low-depth alleles without GATK JVM overhead.
  • Streamlined SAMtools/BCFtools Core: Minimizes I/O bottlenecks and intermediate format conversions by chaining native command-line utilities.
  • SLURM-Optimized Parallelism: Features chunked parallel scatter/gather steps and pre-configured execution profiles (e.g., Rhodotorula cohorts) tuned for high-performance computing clusters.
  • Reproducible & Modular: Implemented in Nextflow DSL2 using standard nf-core modules, Conda environment integration, and automated dynamic resource allocation (task.attempt).
Workflow Nextflow