Skip to main content

A Nextflow Pipeline for Whole-Genome Variant Calling

Bioinformatics
Data engineering
Nextflow
A modular Nextflow DSL2 pipeline that takes paired-end whole-genome sequencing reads to joint-called, filtered variants with BWA-MEM, SAMtools and GATK, built to run on a PBS cluster and to compare two reference genomes.
Published

2024

NoteAt a glance
  • What: a modular Nextflow DSL2 pipeline that turns raw paired-end whole-genome sequencing (WGS) reads into a joint-called, quality-filtered set of variants
  • Tools: Nextflow, BWA-MEM, SAMtools and GATK4, with FastQC and MultiQC for read quality control
  • Where: built for my BSc Honours project at Stellenbosch University (2024) and run on the Centre for High Performance Computing (CHPC) cluster
  • Scope of this page: how the pipeline is designed and how it works. It contains no study data and no results.

View the code on GitHub

Why I built it

My honours project compared variant calls made against the standard human reference genome, GRCh38 (Schneider et al., 2017), with calls made against a population-specific reference genome.

That comparison set three requirements:

  1. Reproducible: every sample goes through exactly the same steps, with the same tools and settings.
  2. Reference-agnostic: the same pipeline runs against either reference with a single switch, keeping the outputs apart.
  3. HPC-friendly: whole genomes produce very large intermediate files, so the pipeline has to fit a shared cluster’s storage, memory and wall-time limits.

Nextflow (Di Tommaso et al., 2017) fits all three: each step is an isolated process, runs are resumable, and the same code can run on a laptop or a PBS cluster.

How it works

flowchart LR
    A[Paired-end<br/>FASTQ] --> B[Align<br/>BWA-MEM]
    B --> C[SAM → BAM<br/>sort, index<br/>SAMtools]
    C --> D[Base quality<br/>recalibration<br/>GATK BQSR]
    D --> E[Per-sample<br/>GVCF<br/>HaplotypeCaller]
    E --> F[Joint genotyping<br/>GenomicsDBImport<br/>GenotypeGVCFs]
    F --> G[Variant filtering<br/>VQSR<br/>SNPs + indels]
    G --> H[Filtered<br/>cohort VCF]

Read quality was checked with FastQC (Andrews, 2010) and summarised across samples with MultiQC (Ewels et al., 2016) before reads entered the pipeline. The pipeline then follows GATK’s best-practice workflow for germline short variants (Van der Auwera and O’Connor, 2020):

Step Tool (version) What it does
Align BWA-MEM 0.7.7 (Li and Durbin, 2009; Li, 2013) Maps each sample’s read pairs to the chosen reference, adding a read group (sample, library, MGI platform) that GATK needs downstream
Convert, sort, index SAMtools 1.9 (Danecek et al., 2021) Converts the text SAM output to compressed BAM, sorts reads by position and indexes them
Recalibrate base qualities GATK 4.6 BaseRecalibrator and ApplyBQSR (McKenna et al., 2010; DePristo et al., 2011) Corrects systematic errors in the sequencer’s base quality scores, using dbSNP (Sherry, Ward and Sirotkin, 1999) as known variant sites
Call variants per sample GATK 4.6 HaplotypeCaller, GVCF mode (Poplin et al., 2017) Calls variants by local re-assembly and records every position, so samples can be genotyped together later
Joint genotyping GATK 4.6 GenomicsDBImport and GenotypeGVCFs (Poplin et al., 2017) Combines all samples’ GVCFs, chromosome by chromosome (1–22, X, Y), and genotypes them as one cohort
Filter variants GATK 4.2.6 VariantRecalibrator and ApplyVQSR (DePristo et al., 2011) Learns what true variants look like from trusted resources and filters the rest at 99% truth sensitivity, with separate models for SNPs and indels

The variant filtering step trains on well-established resources: HapMap 3.3 (International HapMap 3 Consortium, 2010) and the 1000 Genomes Omni 2.5 genotypes (The 1000 Genomes Project Consortium, 2015) for SNPs, the Mills and 1000 Genomes gold-standard indels (Mills et al., 2011) for indels, and dbSNP (Sherry, Ward and Sirotkin, 1999) as known sites. The filtered SNPs and indels are then merged into a single indexed cohort VCF.

Engineering details

One module per step

Each step is a separate DSL2 module (modules/local/<step>/main.nf), and the main workflow only wires them together. That keeps each step small, testable and easy to swap out:

script.nf (excerpt)
workflow {
    aligned_sams = ALIGN(samples_ch, reference_ch)
    aligned_bams = SAM_TO_BAM(aligned_sams)
    sorted_bams  = SORT_BAM(aligned_bams)
    indexed_bams = INDEX_BAM(sorted_bams.sorted_bam)

    bqsr_bams = BQSR(indexed_bams, reference_ch, reference_index_ch,
                     reference_dict_ch, dbsnp_ch, dbsnp_index_ch)

    gvcfs = HAPLOTYPECALLER(bqsr_bams, reference_ch,
                            reference_index_ch, reference_dict_ch)
    // ... joint genotyping and variant filtering
}

One switch for the reference genome

A single parameter chooses the reference. Outputs are named and stored separately per reference, so a GRCh38 run and a population-reference run never overwrite each other:

modules/local/haplotypecaller/main.nf (excerpt)
publishDir "${params.outdir}/haplotypecaller/${params.use_sac_ref ? 'sac' : 'grch38'}",
           mode: 'symlink'

Dry runs without data

Every module has a stub: block that creates empty placeholder outputs. Running nextflow run script.nf -stub checks that all the channels connect correctly in seconds, without real data and without using cluster time:

modules/local/bqsr/main.nf (excerpt)
stub:
"""
touch ${sampleId}.recal.bam
touch ${sampleId}.recal.bai
"""

Fitting a shared cluster

  • Storage: an aligned whole genome in SAM format can be tens of gigabytes. Intermediate outputs are published as symbolic links (mode: 'symlink') rather than copies, and only the final filtered VCF is copied.
  • Resources: memory and CPUs are set per process, so heavy steps such as alignment get what they need without over-allocating light ones.
  • Software: each process loads its own pinned tool version through the cluster’s environment modules.
nextflow.config (excerpt, cluster details replaced)
process {
    queue          = 'normal'
    clusterOptions = '-P <project-code>'

    withName: 'ALIGN' {
        beforeScript = 'module load BWA/0.7.7'
        memory = '38 GB'
    }
    withName: 'HAPLOTYPECALLER' {
        beforeScript = 'module load gatk/4.6.0.0'
        memory = '8 GB'
        cpus   = 5
    }
}

Joint calling that scales

Calling each sample to a GVCF first, then genotyping the cohort together, means that adding samples later only needs the new samples to be called. It also improves genotype accuracy at sites with low coverage in some samples (Poplin et al., 2017).

What I’d do differently

Looking back, these are the changes I would make:

  • Mark duplicates before recalibration. GATK’s best practices include MarkDuplicates between alignment and BQSR (Van der Auwera and O’Connor, 2020); the pipeline goes straight from sorting to BQSR.
  • Stream instead of writing SAM files. Piping bwa mem directly into samtools sort would skip the large intermediate SAM files entirely.
  • Make quality control part of the pipeline. FastQC and MultiQC were configured but run separately; wiring them in would make each run fully self-contained.
  • Use containers and profiles. Docker or Singularity containers and Nextflow configuration profiles, instead of cluster-specific modules and paths, would let the pipeline run unchanged on any system.
  • Use one GATK version throughout. Variant filtering ran on GATK 4.2.6 while the other GATK steps used 4.6.
  • Compare references against a truth set. Counting shared and unique variants between two references shows that the call sets differ, but not which is more accurate: a reference that already contains a population’s common alleles will report fewer variants at those sites by design. A fair comparison needs a truth set, such as array genotypes for the same samples.

The code

The pipeline code is on GitHub as it was used for the project, with cluster-specific paths and the project code replaced by placeholders. No sample sheet or data is included. Its README lists how to run it and the known issues in this snapshot.

Acknowledgements

This work was done for my BSc Honours in Bioinformatics and Computational Biology. I’m grateful to the people and institutions that made it possible:

My supervisors

Dr C. Uren, Prof M. Möller, Prof D.C. Petersen and Dr Y. Swart, for their guidance throughout the project

Stellenbosch University

Host institution for my BSc Honours

Faculty of Medicine and Health Sciences

Stellenbosch University, for supporting the project

Division of Molecular Biology and Human Genetics

Stellenbosch University, where the project was based

Centre for Bioinformatics and Computational Biology (CBCB)

Stellenbosch University, for bioinformatics support

South African Medical Research Council (SAMRC)

The SAMRC Genomics Platform, for sequencing the data

Centre for High Performance Computing (CHPC)

Computations were performed using facilities provided by the CHPC, South Africa

References

Andrews, S., 2010. FastQC: A quality control tool for high throughput sequence data. Available at: <https://www.bioinformatics.babraham.ac.uk/projects/fastqc/>.
Danecek, P., Bonfield, J.K., Liddle, J., Marshall, J., Ohan, V., Pollard, M.O., Whitwham, A., Keane, T., McCarthy, S.A., Davies, R.M. and Li, H., 2021. Twelve years of SAMtools and BCFtools. GigaScience, [online] 10(2). https://doi.org/10.1093/gigascience/giab008.
DePristo, M.A., Banks, E., Poplin, R., Garimella, K.V., Maguire, J.R., Hartl, C., Philippakis, A.A., Angel, G. del, Rivas, M.A., Hanna, M., McKenna, A., Fennell, T.J., Kernytsky, A.M., Sivachenko, A.Y., Cibulskis, K., Gabriel, S.B., Altshuler, D. and Daly, M.J., 2011. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, [online] 43(5), pp.491–498. https://doi.org/10.1038/ng.806.
Di Tommaso, P., Chatzou, M., Floden, E.W., Barja, P.P., Palumbo, E. and Notredame, C., 2017. Nextflow enables reproducible computational workflows. Nature Biotechnology, [online] 35(4), pp.316–319. https://doi.org/10.1038/nbt.3820.
Ewels, P., Magnusson, M., Lundin, S. and Käller, M., 2016. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics, [online] 32(19), pp.3047–3048. https://doi.org/10.1093/bioinformatics/btw354.
International HapMap 3 Consortium, 2010. Integrating common and rare genetic variation in diverse human populations. Nature, [online] 467(7311), pp.52–58. https://doi.org/10.1038/nature09298.
Li, H., 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. https://doi.org/10.48550/ARXIV.1303.3997.
Li, H. and Durbin, R., 2009. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, [online] 25(14), pp.1754–1760. https://doi.org/10.1093/bioinformatics/btp324.
McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., Garimella, K., Altshuler, D., Gabriel, S., Daly, M. and DePristo, M.A., 2010. The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, [online] 20(9), pp.1297–1303. https://doi.org/10.1101/gr.107524.110.
Mills, R.E., Pittard, W.S., Mullaney, J.M., Farooq, U., Creasy, T.H., Mahurkar, A.A., Kemeza, D.M., Strassler, D.S., Ponting, C.P., Webber, C. and Devine, S.E., 2011. Natural genetic variation caused by small insertions and deletions in the human genome. Genome Research, [online] 21(6), pp.830–839. https://doi.org/10.1101/gr.115907.110.
Poplin, R., Ruano-Rubio, V., DePristo, M.A., Fennell, T.J., Carneiro, M.O., Van der Auwera, G.A., Kling, D.E., Gauthier, L.D., Levy-Moonshine, A., Roazen, D., Shakir, K., Thibault, J., Chandran, S., Whelan, C., Lek, M., Gabriel, S., Daly, M.J., Neale, B., MacArthur, D.G. and Banks, E., 2017. Scaling accurate genetic variant discovery to tens of thousands of samples. [online] https://doi.org/10.1101/201178.
Schneider, V.A., Graves-Lindsay, T., Howe, K., Bouk, N., Chen, H.-C., Kitts, P.A., Murphy, T.D., Pruitt, K.D., Thibaud-Nissen, F., Albracht, D., Fulton, R.S., Kremitzki, M., Magrini, V., Markovic, C., McGrath, S., Steinberg, K.M., Auger, K., Chow, W., Collins, J., Harden, G., Hubbard, T., Pelan, S., Simpson, J.T., Threadgold, G., Torrance, J., Wood, J.M., Clarke, L., Koren, S., Boitano, M., Peluso, P., Li, H., Chin, C.-S., Phillippy, A.M., Durbin, R., Wilson, R.K., Flicek, P., Eichler, E.E. and Church, D.M., 2017. Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly. Genome Research, [online] 27(5), pp.849–864. https://doi.org/10.1101/gr.213611.116.
Sherry, S.T., Ward, M. and Sirotkin, K., 1999. dbSNP—database for single nucleotide polymorphisms and other classes of minor genetic variation. Genome Research, [online] 9(8), pp.677–679. https://doi.org/10.1101/gr.9.8.677.
The 1000 Genomes Project Consortium, 2015. A global reference for human genetic variation. Nature, [online] 526(7571), pp.68–74. https://doi.org/10.1038/nature15393.
Van der Auwera, G.A. and O’Connor, B.D., 2020. Genomics in the cloud: Using Docker, GATK, and WDL in Terra. 1st ed. O’Reilly Media.