Genome Analysis Toolkit (GATK) . -- GVA2020


Overview

The Genome Analysis Toolkit (GATK) is a set of programs developed by the broad institute with an extensive website. As mentioned in the final presentation, it has the ability to perform much of the analysis required for calling genomic variants as well as many many other things. Why you may ask yourself did this magical tool only appear on the final day of the class? GATK uses read mappers, read aligners, variant callers, and all the other things (or similar things) that you have been introduced to throughout the course so we have actually been going over what you needed to know in smaller more digestible chunks.

This tutorial is quite small and does not showcase but the smallest drop in a bucket of what GATK is capable of doing. This is because the broad itself has developed many many tutorials for all the different things GATK does and extensive forums are available if the tutorials are not enough to get you through what you are trying to do. Finally, as the makers of the software they have put out and maintain what they regard as the best way to use their product in the form of 'best practices'. If you are going to use GATK, its a real real real good idea to make sure you are following their best practices because that is a situation where people will raise a big eyebrow if you say you are going against the flow.

While GATK is great, one stop shops often are often not the best at everything they do, don't be afraid to use other programs. Particularly following what other researchers are doing in your field

Objectives

  1. Load GATK on lonestar
  2. Use the sample data provided by the broad (through TACC) to verify that TACC is working
  3. Explore a little of what is under the hood.

Tutorial: Loading GATK

While you may think based on the overview that GATK is an obvious choice for a module on TACC, you may be surprised to learn that seemingly every other year TACC removes it as a module, and this is a bad year. On the plus side, it means that once we install it for you locally, the only issue will be if you need to update the version, and recent changes to GATK have made it much easier to work with.

Instaling GATK and verifying it is functioning
# set up directories
mkdir $WORK/src
mkdir -p $HOME/local/bin


# download file and extract it
cd $WORK/src   
wget https://github.com/broadinstitute/gatk/releases/download/4.1.7.0/gatk-4.1.7.0.zip
unzip gatk-4.1.7.0.zip

# copy executables to somewhere already in your $PATH variable (remember we set this up on Monday in your .bashrc file)
cp gatk-4.1.7.0/gatk $HOME/local/bin
cp gatk-4.1.7.0/*.jar $HOME/local/bin

# verify correctly installed
cds
gatk -help # if this does not output a large list of colored text, try the following command and if that does not output colored text get my attention
gatk --list

If you see 316 lines of a long scrolling output detailing some copyright information and a bunch of different commands everything is correctly loaded. While individual tools will require different options and the program itself takes many different options only 3 things are ALWAYS required:

flagDescription

Tool name, what tool are you trying to use
-RReference sequence file
-IInput bam file

Stealing a nice mnemonic devices from a GATK toturial (which is condensed below), these 3 arguments don't have to be in this order, but if you learn them in this order, you will be able to remember them if you TRI. Remember, specific tools will require additional arguments.

Getting sample data

Rather than using sample data specifically for this tutorial, we will instead do a small tutorial based on our read mapping tutorial from day 2 of the course. Assuming you completed that tutorial you the following tutorial should work.

You are trying to copy the SRR030257.sam file from the $SCRATCH/GVA_bowtie2_mapping/bowtie2/SRR030257.sam and the NC_012967.1.fasta file from the $SCRATCH directory
mkdir $SCRATCH/GVA_GATK
cd $SCRATCH/GVA_GATK
cp /scratch/01821/ded/GVA_bowtie2_mapping/bowtie2/SRR030257.sam .
cp /scratch/01821/ded/GVA_bowtie2_mapping/NC_012967.1.fasta .

Next you need to convert the .sam file to a .bam file.

Refresher on how to convert .sam files into .bam files
samtools view -S -b SRR030257.sam > SRR030257.bam 


# Do you remember what tutorial we used this command in before?


Tutorial: Use GATK to count the number of reads in a bam file

Using the following information we will use gatk the CountReads tool to count the number of reads in the SRR030257.bam file which was from the NC_012967.fasta reference file. Pay attention to the the words in bold and the table/discussion in the previous tutorial section and see if you can figure out how to do this on your own.

 Check your answer

gatk  CountReads -R NC_012967.1.fasta -I SRR030257.bam

 Did this command do what you expected? click here and continue once you decide if you think it worked as expected

The following line tells us things did not go as intended.

The following lines tell us this is not working as intended

***********************************************************************


A USER ERROR has occurred: Fasta index file file:///scratch/01821/ded/GVA_GATK/NC_012967.1.fasta.fai for reference file:///scratch/01821/ded/GVA_GATK/NC_012967.1.fasta does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.


***********************************************************************

The text clearly tells us that a fasta index file is missing, and from the gatk link, we see that they are using the samtools faidx command to generate their indexes of fasta files.

index command
samtools faidx NC_012967.1.fasta

Now we repeat our CountReads command, and see that the results have indeed changed.

 Did this command do what you expected? click here and continue once you decide if you think it worked as expected

The following line tells us things still did not go as intended.

The following lines tell us this is not working as intended

***********************************************************************


A USER ERROR has occurred: Fasta dict file file:///scratch/01821/ded/GVA_GATK/NC_012967.1.dict for reference file:///scratch/01821/ded/GVA_GATK/NC_012967.1.fasta does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.


***********************************************************************

That link contains a discussion about a .dict file that is required for fasta references that might normally be generated by GATK but as we did not do the read mapping inside of GATK we lack it. Based on the reading on the linked web page I came up with the following solution to create that dictionary file.

Note that I handle creating the dictionary file directly through GATK rather than through isolating a specific .jar file
gatk CreateSequenceDictionary -R NC_012967.1.fasta -O NC_012967.1.dict

Now if we retry our CountReads command we get very different output (note that yours will have subtle differences in things like the names of directories):

INFO: Failed to detect whether we are running on Google Compute Engine.
01:26:37.269 INFO  CountReads - ------------------------------------------------------------
01:26:37.269 INFO  CountReads - The Genome Analysis Toolkit (GATK) v4.1.7.0
01:26:37.269 INFO  CountReads - For support and documentation go to https://software.broadinstitute.org/gatk/
01:26:37.269 INFO  CountReads - Executing as ded@login2 on Linux v4.4.103-6.38-default amd64
01:26:37.269 INFO  CountReads - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_151-b12
01:26:37.269 INFO  CountReads - Start Date/Time: June 24, 2020 1:26:37 AM CDT
01:26:37.269 INFO  CountReads - ------------------------------------------------------------
01:26:37.270 INFO  CountReads - ------------------------------------------------------------
01:26:37.270 INFO  CountReads - HTSJDK Version: 2.21.2
01:26:37.270 INFO  CountReads - Picard Version: 2.21.9
01:26:37.270 INFO  CountReads - HTSJDK Defaults.COMPRESSION_LEVEL : 2
01:26:37.270 INFO  CountReads - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
01:26:37.270 INFO  CountReads - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
01:26:37.270 INFO  CountReads - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
01:26:37.270 INFO  CountReads - Deflater: IntelDeflater
01:26:37.270 INFO  CountReads - Inflater: IntelInflater
01:26:37.270 INFO  CountReads - GCS max retries/reopens: 20
01:26:37.270 INFO  CountReads - Requester pays: disabled
01:26:37.270 INFO  CountReads - Initializing engine
01:26:37.589 INFO  CountReads - Done initializing engine
01:26:37.589 INFO  ProgressMeter - Starting traversal
01:26:37.589 INFO  ProgressMeter -        Current Locus  Elapsed Minutes       Reads Processed     Reads/Minute
01:26:48.143 INFO  CountReads - 7600360 read(s) filtered by: WellformedReadFilter 


01:26:48.146 INFO  ProgressMeter -             unmapped              0.2                     0              0.0
01:26:48.147 INFO  ProgressMeter - Traversal complete. Processed 0 total reads in 0.2 minutes.
01:26:48.147 INFO  CountReads - CountReads counted 0 total reads
01:26:48.147 INFO  CountReads - Shutting down engine
[June 24, 2020 1:26:48 AM CDT] org.broadinstitute.hellbender.tools.CountReads done. Elapsed time: 0.18 minutes.
Runtime.totalMemory()=2804940800
Tool returned:
0


What in all that are we actually looking for you might ask? 

01:26:48.143 INFO  CountReads - 7600360 read(s) filtered by: WellformedReadFilter

01:26:48.147 INFO  ProgressMeter - Traversal complete. Processed 0 total reads in 0.2 minutes.

This tells us that the bam file contains 7600360 total reads and that none were removed by any filtering options. The lack of anything being removed should make sense since we didn't try to filter anything out.

You can check this against the original fastq files SRR030257_1.fastq and SRR030257_1.fastq (check them in your GVA_bowtie2_mapping directory)
grep -c "^+$" $SCRATCH/GVA_bowtie2_mapping/SRR030257_1.fastq
grep -c "^+$" $SCRATCH/GVA_bowtie2_mapping/SRR030257_2.fastq

While we haven't discussed sam format in much detail, each read gets its own line, and if you compare the .sam file and the original fastq file listed above, you see that each line on the sam file seems to start with 'SRR030257' followed by a number related to the read. This gives us a base handle to check with grep

help checking sam file
grep -c "^SRR030257" SRR030257.sam

In the case of the sam file, you can see that the 7600360 exactly matches what we saw on the GATK output, while the fastq files each show 3800180 (or 1/2 the total reads). While this is somewhat trivial and 


As mentioned this is a very small introduction to GATK adapted from one of the broad's tutorials which can be found here http://gatkforums.broadinstitute.org/gatk/discussion/1209/howto-run-the-gatk-for-the-first-time. Feel free to explore that link and the other tutorial links for taking GATK further.

Return to GVA2020 home page.