GATK

A Nextflow wrapper for the Genome Analysis Toolkit (GATK), modified from the pipeline described in the Bioinformatic Workbook: GATK Best Practices Workflow for DNA-Seq.

This Nextflow wrapper extends GATK's functionality from DNAseq data, to enabling variant calling from RNA-Seq and long-read sequencing data. It provides a comprehensive solution for diverse sequencing platforms, enhancing variant analysis capabilities.

Installation

You will need a working version of nextflow, see here on how to install nextflow. Nextflow modules are avialable on some of the HPC computing resources.

nextflow run isugifNF/GATK --help

Singularity Container

Tools required for the workflow are included in the container aseetharam/gatk:latest and should be automatically pulled by nextflow. (Will only need to run singularity pull if website connection is unstable.)

Test Dataset

A simple DNAseq test dataset (test-data) is available on ISU Box. This dataset contains a small genome (portion of chr1, B73v5 ), and Illumina short reads for 26 NAM lines (including B73) and B73Ab10 line (27 lines total). Only the reads that map to the region of the v5 genome is included, so that this can be tested quickly. There are examples of multiple files belonging to same NAM line as well as single file per NAM line to make sure both conditions works correctly. The end VCF file should have exactly 27 individuals (lines) in them.

wget https://iastate.box.com/shared/static/wt85l6s4nw4kycm2bo0gpgjq752osatu.gz
tar -xf wt85l6s4nw4kycm2bo0gpgjq752osatu.gz

ls -1 test-data/
#> fastq          # <= folder of paired reads
#> read-group.txt
#> ref            # <= folder containing one genome reference

Running the Pipeline

This pipeline can process DNAseq, RNAseq and PacBio long read data, which require different arguments and different input files. Please jump to the section based off of your input below.

GATK variant calling for DNAseq data

Because this pipeline was initially developed to process DNAseq data, we have provided a test DNAseq dataset. You can fetch the pipeline and the test-data folder.

# Fetch pipeline
nextflow run isugifNF/GATK --help

# Fetch the test-data folder from ISU box
wget https://iastate.box.com/shared/static/wt85l6s4nw4kycm2bo0gpgjq752osatu.gz
tar -xf wt85l6s4nw4kycm2bo0gpgjq752osatu.gz

The general format of a run with the pipeline is to provide a genome file (--genome) and Illumina Paired-End Reads files (--reads or --reads_file).

nextflow run isugifNF/GATK \
  --genome test-data/ref/b73_chr1_150000001-151000000.fasta \
  --reads "test-data/fastq/*_{R1,R2}.fastq.gz" \
  -profile slurm,singularity \
  -resume

or

nextflow run isugifNF/GATK \
  --genome test-data/ref/b73_chr1_150000001-151000000.fasta \
  --reads_file read-path.txt \
  -profile slurm,singularity \
  -resume

If you are on a HPC (Nova/Condo/Ceres/Atlas), it is highly recommend to use the submit_nf.slurm script. The # === Load Modules section will need to be modified to get nextflow and singulariy running.

Using the --reads_file

Instead of using a pattern to specify paired reads ("test-data/fastq/*_{R1,R2}.fastq.gz"), we can use a tab-delimited textfile to specify the path to left and right files. The textfile should contain three columns: readname, left read path, right read path.

In our case, for the test-data run the following one-liner to generate a read-path.txt.

for f in test-data/fastq/*_R1.fastq.gz; do echo -e "$(basename $f |cut -f 1 -d "_")\t$(realpath $f)\t$(realpath $f | sed 's/_R1.fastq.gz/_R2.fastq.gz/g')"; done > read-path.txt

DNAseq final output

The final output will be in a results folder. SNPs will be in the VCF file, probably the file with the longest name (e.g. first-round_merged_snps-only_snp-only.pass-only.vcf).

results/
  |_ 01_MarkAdapters/            #<= folders contain intermediate files
  |_ 02_MapReads/
  |_ 03_PrepGATK/
  |_ 04_GATK/
  |_ 05_FilterSNPs/
  |  |_ first-round_merged_snps-only_sorted_snp-only.pass-only.vcf     #<= final SNP file
  |
  |_ report.html
  |_ timeline.html               # <= runtime information for all processes

GATK variant calling for RNAseq

In order to process RNAseq, we use the STAR aligner. To use the RNAseq variant calling pipeline, make the following parameter changes.

nextflow run isugifNF/GATK \
  --seq "rna" \
  --genome "genome/genome.fa" \
  --reads "data/*.{r1,r2}.fq.gz" \
  --gtf "genome/genome.gtf" \
  -profile slurm,singularity \
  -resume

RNAseq final output

The final output will be in a results folder. SNPs will be in the VCF file, probably the file with the longest name (e.g. first-round_merged_snps-only_snp-only.pass-only.vcf).

results/
  |_ 01_MarkAdapters/            #<= folders contain intermediate files
  |_ 02_MapReads/
  |_ 03_PrepGATK/
  |_ 04_GATK/
  |_ 05_FilterSNPs/
  |  |_ first-round_merged_snps-only_sorted_snp-only.pass-only.vcf     #<= final SNP file
  |
  |_ report.html
  |_ timeline.html               # <= runtime information for all processes

GATK variant calling for long read data

In order to process long read data obtained from the PacBio platform, we use the pbmm2 aligner. To use the long read variant calling pipeline, make the following parameter changes.

nextflow run isugifNF/GATK \
  --seq "longread" \
  --genome "genome/genome.fa" \
  --reads "data/long-reads.fastq" \
  -profile slurm,singularity \
  -resume

Long read final output

The final output will be in a results folder. SNPs will be in the VCF file, probably the file with the longest name (e.g. output_snp-only.pass-only.vcf).

results/
  |_ 02_MapReads/
  |_ 03_PrepGATK/
  |_ 04_GATK/
  |_ 05_FilterSNPs/
  |  |_ output_snp-only.pass-only.vcf     #<= final SNP file
  |
  |_ report.html
  |_ timeline.html               # <= runtime information for all processes

Example Runs

Some example runs provided to show nextflow output.

Example runs for DNAseq data

Example run for RNAseq data

Example run for long read data