Genome Analysis Toolkit (GATK) . -- GVA2021
- Deatherage, Daniel E
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 and default settings/values may or may not be the best solution for your work. Be sure to understand the options you are using, the options you are not using but could be, and what others doing similar analysis have done in their work.
Objectives
- Install GATK
- Explore a little of how gatk works
- Provide resources for working with GATK tutorials
idev nodes
Note that the entire tutorial can take place inside an idev session as we will be copying data from earlier tutorials rather than the BioITeam. Not only can it take place on an idev node, it MUST take place on an idev node. While preparing this tutorial, accidentally running the gatk commands on the head node resulted in memory errors.
Tutorial: Installing GATK
Not surprisingly, GATK is available through conda. Since it is so involved and complex I would suggest giving it its own environment perhaps "GVA-gatk"
conda create --name GVA-gatk -c bioconda gatk4 conda activate GVA-gatk gatk --version
The Genome Analysis Toolkit (GATK) v4.2.0.0 HTSJDK Version: 2.24.0 Picard Version: 2.25.0
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:
| flag | Description |
|---|---|
| Tool name, what tool are you trying to use | |
| -R | Reference sequence file |
| -I | Input bam file |
Stealing a nice mnemonic devices from a GATK legacy tutorial (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.
Next you need to convert the .sam file to a .bam file. Since we are keeping gatk in its own environment you will need to briefly switch over to a different environment to access the correct program.
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.
The following lines tell us this is not working as intended
***********************************************************************
A USER ERROR has occurred: Fasta dict file file:///scratch/0004/train402/GVA_GATK/NC_012967.1.dict for reference file:///scratch/0004/train402/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.
***********************************************************************
While the provided link surprisingly doesn't work, it is clear that it is failing because it is not finding a .dict file for our reference sequence. Quickly googling "gatk create dict file" takes us to https://gatk.broadinstitute.org/hc/en-us/articles/360036729911-CreateSequenceDictionary-Picard-
gatk CreateSequenceDictionary -R NC_012967.1.fasta -O NC_012967.1.dict
Using the ls command you can see that the dictionary has been created, and the CountReads command will now run correctly. 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:55:48.483 INFO CountReads - ------------------------------------------------------------ 01:55:48.485 INFO CountReads - The Genome Analysis Toolkit (GATK) v4.2.0.0 01:55:48.485 INFO CountReads - For support and documentation go to https://software.broadinstitute.org/gatk/ 01:55:48.494 INFO CountReads - Executing as train402@c455-084.stampede2.tacc.utexas.edu on Linux v3.10.0-957.5.1.el7.x86_64 amd64 01:55:48.494 INFO CountReads - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_152-release-1056-b12 01:55:48.495 INFO CountReads - Start Date/Time: June 17, 2021 1:55:47 AM CDT 01:55:48.495 INFO CountReads - ------------------------------------------------------------ 01:55:48.496 INFO CountReads - ------------------------------------------------------------ 01:55:48.498 INFO CountReads - HTSJDK Version: 2.24.0 01:55:48.498 INFO CountReads - Picard Version: 2.25.0 01:55:48.498 INFO CountReads - Built for Spark Version: 2.4.5 01:55:48.499 INFO CountReads - HTSJDK Defaults.COMPRESSION_LEVEL : 2 01:55:48.499 INFO CountReads - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false 01:55:48.499 INFO CountReads - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true 01:55:48.500 INFO CountReads - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false 01:55:48.500 INFO CountReads - Deflater: IntelDeflater 01:55:48.500 INFO CountReads - Inflater: IntelInflater 01:55:48.501 INFO CountReads - GCS max retries/reopens: 20 01:55:48.501 INFO CountReads - Requester pays: disabled 01:55:48.501 INFO CountReads - Initializing engine 01:55:50.020 INFO CountReads - Done initializing engine 01:55:50.021 INFO ProgressMeter - Starting traversal 01:55:50.023 INFO ProgressMeter - Current Locus Elapsed Minutes Reads Processed Reads/Minute 01:56:35.721 INFO CountReads - 7600360 read(s) filtered by: WellformedReadFilter 01:56:35.725 INFO ProgressMeter - unmapped 0.8 0 0.0 01:56:35.726 INFO ProgressMeter - Traversal complete. Processed 0 total reads in 0.8 minutes. 01:56:35.726 INFO CountReads - CountReads counted 0 total reads 01:56:35.727 INFO CountReads - Shutting down engine [June 17, 2021 1:56:35 AM CDT] org.broadinstitute.hellbender.tools.CountReads done. Elapsed time: 0.80 minutes. Runtime.totalMemory()=2549612544 Tool returned: 0
What in all that are we actually looking for you might ask?
01:56:35.726 INFO CountReads - CountReads counted 0 total reads
But we know this isn't true.
This tells us that there 7600360 total reads which is certainly greater than 0. 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
In the case of the sam file, you can see that the 7600360 exactly matches what we saw from our fastq files, while it is possible there is an error with the sam to bam conversion the file does have a very large 352 M if there are no reads in that file.
So where did our reads actually go in our GATK analysis? A few lines above we see this line:
01:56:35.721 INFO CountReads - 7600360 read(s) filtered by: WellformedReadFilter
This tells us that all the reads were removed by a specific filter set. Accessing the CountReads help information, we see we can actually disable the WellformedReadFilter with the -DF or --disable-read-filter. giving us an updated CountReads command of:
gatk CountReads -R NC_012967.1.fasta -I SRR030257.bam -DF WellformedReadFilter
Now we see a final line of output being returned by the tool of: 7,600,360 which is what we expected!
GATK has a similar tool 'CountBases' which as the name implies, counts the bases in a sam/bam file.
gatk CountBases -R NC_012967.1.fasta -I SRR030257.bam -DF WellformedReadFilter
This tells us that there are 273,612,960 total bases in the file.
As mentioned this is a very small introduction to GATK adapted from one of the broad's tutorials which can be found here. There are numerous other tutorial links for on the GATK website you can begin to work with now that you have verified your installation and begun troubleshooting some of the quirks of the program.
Welcome to the University Wiki Service! Please use your IID (yourEID@eid.utexas.edu) when prompted for your email address during login or click here to enter your EID. If you are experiencing any issues loading content on pages, please try these steps to clear your browser cache.