trust1/gatk repository overview

Test GATK image from Broad with few extras added:
BWA,
VCFtools,
ANGSD,
Stacks,
ts (task spooler)

The perl written pipeline for BAM or VCF files generation:

    ~/callingpipe_GATK.pl -R /downloads/input/Some_Species/reference.fasta --GVCF ON --TRDS 2 --SMPL 5 

Three major steps required before starting this pipeline:

[ 1. ] Get (GATK modified) image from Rust's repo (in Ubuntu):

                             sudo docker pull trust1/gatk:version0.4.1

[ 2. ] Fire docker image. For example on Ubuntu it would be such command:

 
                     sudo docker run -v /your/local/folder:/downloads --name GATKvX -it trust1/gatk:version0.4.1       
                                        ====================                                           

Where ' /your/local/folder' -- the path to the local host folder with one subdirectory /your/local/folder//input where raw fast deposited is required:

1. /mnt/rust/temp/input -- source with demultiplexed FASTQs and the reference *.fasta in subdirectory, this sub folder name will be used as species name for the reference ids.
   /your/local/folder/input/fishyfish/ncbiID01.fasta -- species name is 'fishyfish'

2. /your/local/folder/output -- folder with temporary and produced result files. That is the place where files will be collected.

Expecting input FASTQ file names format like these ones:
SP20_indexing10_FCBARCODE_L001_R1_001.fastq.gz
SP20_indexing10_FCBARCODE_L001_R2_001.fastq.gz
"SP20_indexing10_FCBARCODE" will be used as individual sample identifier

[ 3. ] Start the pipeline with following options:

       perl ~/callingpipe_GATK.pl  -R /downloads/input/Some_Species/reference.fasta --GVCF ON --TRDS 2 --SMPL 5 

-R /downloads/input/Some_Species/reference.fasta : provide path on mounted folder. This path to the reference genome fasta inside the container.

--GVCF OFF : skimmed pipeline or full fat run [ --GVCF ON ]. Whole genome sequencing mammalian dataset with 10x coverage of 50 samples will require ~10TB drive and ~100GB RAM on >20CPUs to run [ --GVCF ON ] flag. Time-wise such project could take up to one month on 'minimalistic' system.

--SMPL 5 : 5 is default number of samples in in paralleled jobs (optional parameter)

--TRDS 2 : 2 is default number of threads in each sample (optional parameter). This threads number is most important for BWA and samtools parallelization. Do not use more threads multiplied by samples than number of core CPUs your system have. Default settings would expect 10 CPUs (5x2). GATK Java machine is not sensitive to this number although it used in JAVA garbage collector and GATK "--native-pair-hmm-threads".

⁠Quick test protocol:

This test small single end sequence dataset in less than 10 minutes for quick tools checks and training purpouse. After you launch your container run those commands inside the container:

# copy test data and reference fasta to working folder:
tar  -xzvf quicktestdata.tar.gz -C /downloads 

# pipeline progress output to the text fuile with nohup
nohup perl ~/callingpipe_GATK.pl -R /downloads/input/refgenome/scaffold_889.fasta --GVCF ON --TRDS 4 --SMPL 3 > /downloads/nohup_gatk.txt&

# follow pipeline progress in text file real time by tail command
tail -f /downloads/nohup_gatk.txt

It should take less than ten minutes to reach the end of the pipeline. If you system has less than 12 co-processors reduce --TRDS 4 to smaller number in way combined --TRDS and --SMPL are not exceeding your hardware specs. In example above 4x3=12 CPUs are required.

⁠ Major output files (test data example)

allsamples.FINAL.vcf -- full fat pipeline product unfiltered vcf generated with haplotype caller from allsamples.g.vcf. This is final raw VCF file (no filters applied) produced by haplotype caller from combined g.vcf dataset.

AllSamplesCombined.1.vcf.gz -- Skim pipeline product generated from merging individual (single sample VCF) please note individual non calls at the common sites filled as reference.
This could be too brave for some dataset. Use this file with caution.

AllSamplesCombined.g.vcf -- g.vcf file in case in future you may want combine different datasets/populations. Merging different g.vcf could save lot of computer time rather than reruning pipeline from scratch on combined datasets.

1968527trim_scaffold889
1968527_scaffold889
1968519_scaffold889 -- folders with samples individual temporary files (sam, bam, vcf, flagstaf).

⁠Temporary nonessential output in /downloads/output/ on one sample (1968519_scaffold889) example:

1968519_scaffold889/1968519_scaffold889-pe.sam -- BWA mapped file uncompressed
1968519_scaffold889/1968519_scaffold889-pe.bam -- compressed and sorted alignment
1968519_scaffold889/1968519_scaffold889.1.bam -- duplicates marked bam (ready for most of the callers)
1968519_scaffold889/1968519_scaffold889.2.bam -- base quality recalibrated bam GATK haplotypecaller grade
1968519_scaffold889/1968519_scaffold889.1.vcf.gz -- 1st pass variants in single sample mode (skim pipeline) only non reference sites considered.

⁠Versions comments:

version0.1 -- no parallel jobs for simultaneous sample processing. Sample going trough process one by one all parallelization done trough multithreading control however most GATK walkers are not efficient with multi threading. This version is not suitable for for whole genome sequencing data with [-GVCF ON] parameter.

version0.2 -- added task spooler for sample queuing and sample parallel calculation. Default settings will run 5 samples with 2 threads each. This means 10 CPUs will be required with about 6GB of memory per CPU is needed (>30GB is expected on default settings). Note: the GATK haplotype caller is not very responsive to threads number. The major productivity parameter is '-SMPL ' however keep in mind memory amount on your hardware. Non parallel Java instances for gatk reserving 16GB. This is amount per single pipeline command: the perl instance.

version0.2.1 -- added stacks toolkit for RAD data: Stacks version-2.41 July 8, 2019

version0.3.0 -- added latest GATK v.4.1.4.1 (Nov 28, 2019) into pipeline increase Java memory allocation for some paralleled steps from 2GB to 6GB. Jan 29 2020

version0.3.1 -- main perl script was modified to handle better memory allocation. Each sample will need at least 4-5GB of free memory. In case you run over free memory limit pipeline will not proceed. Complain message will be returned. In such case reduce this number [--SMPL 5] to more feasible amount. Memory allocation will be done automatically: script will check total free memory available on machine and subtract 4GB from that number to reserve this amount for non-paralleled GATK steps (merging VCFs and converting gVCF to VCF). Then this subtracted amount of total memory will be divided equally across number of paralleled sample slots with extra 15% reserve for the JAVA virtual machine itself.

version0.3.2 -- VCF Files with more then 1M SNVs and fragmented chromosomes (>2M scaffolds) can not be merged to single one either with GATK or RTG tools however BCFtools dealing with such inputs just fine. To generate /downloads/output/AllSamplesCombined.1.vcf now using BCF tools.

version0.4.0 -- Sample now merged with bcftools the development fork version on 1st pass VCF merging and G.VCF. The multiallelic SNPs are supported for the both (.1.vcf and g.vcf) merge steps. The combined 1st pass VCF file (skimmed pipeline) /downloads/output/AllSamplesCombined.1.vcf is merged in the way when missed genotypes replaced by homozygous reference call. This is quite bold assumption but I believe it is still much better then just leave gaps in merged dataset for homozygous reference. The replacement of GATK tools for combining gvcf files allows to use fragmented genomes. Neither one of CombineGVCFs or GenomicsDBImport tools from GATK works efficiently with fragmented genomes with smaller memory machines. The one way around is use GATK GenomicsDBImport with splitting genome by chromosome or smaller regions (~15min per contig ~50GB RAM) and then reassemble with GatherVCFs. Unfortunately this path could be only efficiently done in parallel otherwise ETA for processing 1.6M contigs on GenomicsDBImport 1,000 days in single threaded mode.

version0.4.1 -- pipeline script for raw DArT data with STACKS reference caller added. Same folder structure is expected for the script. The raw reads initially cleaned with DART specific trimming with cutadapt on both flanks (5' and 3'). The major intentions of having separate DART pipeline is compare GATK and STACKS produced VCF. Although GATK product seems to be more comprehensive the STAKS pipeline is much faster (hrs vs days) and it also have some flexibility with output formats and connectivity to other third party tools like structure, phylip, popgen etc. Launch GATK container as described above and then start DART/Stacks script with following command:

 perl ~/callingpipe_STACKS.pl  -R /downloads/input/Some_Species/reference.fasta  --TRDS 4 --SMPL 2

Input raw fastq files from DArT expected to be named like this:

 1234567.FASTQ.gz   
 1234568.FASTQ.gz   
 1234569.FASTQ.gz   

Collect pipeline products VCF, structure and phylip formatted files in output folder: /downloads/output

1. /downloads/output/DArT_stackscalled.vcf.gz
   STACKS reference based genotype calls with default stacks settings. All samples from raw fastq used for SNP catalogue and treated as single population. All samples included.

/downloads/output/DArT_stacks_1pop.structure
Structure format note all samples labels as single population

/downloads/output/DArT_stacks_var.phylip
Phylip input file for making phylogenetics tree.

1. /downloads/output/tmp
   This is temporary folder with trimmed reads: *.trmd.fq.gz and BAM files.

⁠Methods.

⁠Variant calling.

The variant calling protocol was performed according to best practice recommendation for non-model organisms from developers of Genome Analysis Tool Kit (GATK v4.1.0.0) [1]. Briefly raw reads were mapped to the reference genomes with BWA MEM algorithm (Version: 0.7.17-r1198-dirty) [2]. The produced alignment BAM files were sorted deduplicated with GATK and initially called with single sample mode by GATK Haplotype caller. At this stage quality control estimation for each sample was performed and outlier samples with low reads or quality metrics were excluded. Then produced variants from selected samples were merged to single VCF and used for database to perform base recalibration step proposed by GATK team. Those base quality recalibrated BAM files used for Haplotype caller in ERC mode. This step generates GVCF (genomic variant calling files) which used for cohort genotyping after combining individual gvcf into multisample gvcf. Filtering and statistics for produced variants performed with VCFtools (0.1.15) [3] and GATK. The containerized automated workflow and documentation available on docker repository (trust1/gatk:version0.4.1):
[ https://hub.docker.com/r/trust1/gatk⁠ ].

###References

• [1] Van der Auwera GA, Carneiro MO, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A, Jordan T, Shakir K, Roazen D, Thibault J, Banks E, Garimella KV, Altshuler D, Gabriel S, DePristo MA. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinformatics. 2013;43(1110):11.10.1-11.10.33. doi: 10.1002/0471250953.bi1110s43. PMID: 25431634; PMCID: PMC4243306.
• [2] Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009 Jul 15;25(14):1754-60. doi: - 10.1093/bioinformatics/btp324. Epub 2009 May 18. PMID: 19451168; PMCID: PMC2705234.
• [3] Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, McVean G, Durbin R; 1000 Genomes Project Analysis Group. The variant call format and VCFtools. Bioinformatics. 2011 Aug 1;27(15):2156-8. doi: 10.1093/bioinformatics/btr330. Epub 2011 Jun 7. PMID: 21653522; PMCID: PMC3137218.