GATK使用簡介 Part1/2
GATK(全稱The Genome Analysis Toolkit)是Broad Institute開發的用於二代重測序數據分析的一款軟件，里面包含了很多有用的工具。
網址：http://www.broadinstitute.org/gsa/wiki/index.php/Home_Page
前段時間剛發布了2.x版本的，最近幾天都在不斷更新，網址也搬遷到
http://www.broadinstitute.org/gatk/
不過目前為止原先的文檔還沒有全部轉移過去。2.0版本有一些比較大的更新，同時增加了一些新的功能，具體內容可以參考http://gatk.vanillaforums.com/discussion/67/release-notes-for-gatk-version-2-0，由於2.0目前依然出於比較初步的階段，下面還是以1.x版本進行介紹(后面涉及到的大部分應該變動都不大)。
================================== 下載 =================================
程序本身可以通過主頁上的下載鏈接進去下載：http://www.broadinstitute.org/gsa/wiki/index.php/Downloading_the_GATK
或者可以訪問ftp地址進行下載：ftp://ftp.broadinstitute.org/pub/gsa/GenomeAnalysisTK/
想要得到更新提醒的話可以添加郵件提醒(頁面現在貌似刪掉了，可能是換地址了，麻煩需要的自己找下吧，或者下次我再補上= =)
從以上兩個地址下載到的是已經編譯好的版本，可以直接運行，如果需要下載源程序，可以去github上去找https://github.com/broadgsa/gatk/commits/master，或者linux底下可以通過git工具來下載，下載好的源程序再通過ant進行編譯即可，此處不再贅述。
================================== 運行 =================================
后面以GATK 1.x的最后一個版本1.6-13進行介紹，GATK是一個java程序，解開編譯好的壓縮包，我們可以看到下面有幾個文件：
$ ls
AnalyzeCovariates.jar GenomeAnalysisTK.jar resources
兩個后綴名為jar的文件和一個名為resources的文件夾，這兩個.jar文件就是后面我們要用到java程序，主要要用的就是GenomeAnalysisTK.jar這個程序，resources文件夾里面有一些example的數據，可以用來進行一些測試，另外GATK網站上面還提供了跟human基因組相關的很多數據下載(http://www.broadinstitute.org/gsa/wiki/index.php/GATK_resource_bundle，例如dbsnp的數據等等)，在做跟human相關的一些分析的時候可以直接從上面下載。
Java程序可以在任何安裝有Java平台的系統上運行，運行的時候語法如下：
java -jar <program.jar>
-jar這個參數是必須有的，后面跟你的java程序，例如我們這邊就是
java -jar /home/biosoft/GenomeAnalysisTK-1.6-13/GenomeAnalysisTK.jar …
這邊在引用程序位置時推薦使用絕對路徑，然后除了-jar這個參數以外，還有一個經常用到的參數就是用來設定內存上限的-Xmx參數，用法如下：
java -Xmx<size> -jar …
這邊的size可以用m、g等等做單位，例如
java -Xmx4g -jar /home/biosoft/GenomeAnalysisTK-1.6-13/GenomeAnalysisTK.jar …
就表示把內存的最大使用量限制在4G，對於處理一些比較大的文件，可以適當的把這個值調高一點，來提高運算效率或者防止內存不足程序無法運行。
================================== 全局參數 =================================
GenomeAnalysisTK.jar這一個程序里面包含了很多小的子程序，可以通過
java -jar /home/biosoft/GenomeAnalysisTK-1.6-13/GenomeAnalysisTK.jar --help
來查看具體內容，更詳細的說明請參照
http://www.broadinstitute.org/gsa/gatkdocs/release/
然后簡單的介紹一下幾個全局的參數，這些參數是基本上都會用到的幾個參數
--analysis_type / -T (required String)
這個參數是必須的，用來設定選用的子程序(GATK的參數一般都可以通過一個長參數名或者一個短參數名來進行設定，在每個用法頁面的Argument details里面都會用/給出兩種名稱)，所以一個最基本的GATK命令行的樣子應該是
java -jar/home/biosoft/GenomeAnalysisTK-1.6-13/GenomeAnalysisTK.jar \
-T <String> \
……
--reference_sequence / -R (File)
用來指定fasta格式的參考基因組序列，GATK需要reference序列是經過index的，而且需要兩個index文件，一個是后綴名為.fai的，另外一個是后綴名稱為.dict的，缺少這些文件，或者兩個文件中的內容不一致都可能導致程序報錯，一般情況下沒有這兩個文件GATK在運行時會先生成，如果無法生成或者生成的有問題，也可以通過samtools和picard來手動生成這兩個文件，例如我們的reference序列名稱叫ref.fasta，則用下面的命令行可以生成我們需要的這兩個文件：
Bash語言: 高亮代碼由發芽網提供
## 第一步先用samtools進行index，會在當前目錄下生成一個ref.fasta.fai的文件
samtools faidx ref.fasta
## 在得到.fai文件的基礎上，第二步可以用picard (http://picard.sourceforge.net/)
## 里面的CreateSequenceDictionary這個工具進行index，可以生成另外一個后綴名為.dict
## 的index文件，這一步也可以不做，缺少這個文件的時候運行GATK會自動先生成這個文件，
## 而生成這個文件的過程實際上調用的就是picard里面的這個工具
java -jar /home/biosoft/picard-tools-1.70/CreateSequenceDictionary.jar \
REFERENCE=ref.fasta \
OUTPUT=ref.dict
--input_file / -I (List[String] with default value [])
用來設定輸入的SAM或BAM文件, GATK要求BAM文件也是經過index的，所以如果缺少index文件的話則需要先通過samtools對BAM文件進行index：
samtools index sample.bam
在當前目錄下則會生成sample.bam.bai這個index文件
--intervals / -L (List[IntervalBinding[Feature]])
用來指定某些區域，所有操作只對這些區域進行而不是整個基因組
上面兩個參數都可以設多次，例如在同時處理多個BAM文件的時候，需要通過
-I Sampe01.bam -I Sample02.bam….來進行設定。
--logging_level / -l (String with default value INFO)
用來控制屏幕上顯示的記錄信息的多少，總共有3個等級，INFO, ERROR和FATAL，設定成INFO就表示輸出所有的記錄到屏幕，ERROR就表示只輸出ERROR或FATAL的記錄，FATAL表示只輸出FATAL的記錄，默認等級是INFO，當屏幕上輸出內容太多的時候可以改成ERROR (例如在運行VariantFiltration的時候，默認情況下會輸出大量的WARNING信息)
--num_threads / -nt (Integer with default value 1)
設定線程數，部分(應該很少)GATK子程序是支持多線程的，可以設置這個參數，如果不支持這個參數但是設置了的話運行時會報錯。
--gatk_key / -K (File)
GATK為了不斷改進希望從用戶那邊得到比較多的反饋信息，所以每次程序運行完以后都會生成一些報告發送回服務器，對於有些缺少網絡環境或者不希望發送報告的，可以去申請這個key文件http://www.broadinstitute.org/gsa/wiki/index.php/Phone_home，設置了這個參數以后可以跳過最后的這個步驟。(無法訪問網絡之類的情況下也可以無視，INFO等級下最后完成后會報HttpMethodDirector - I/O exception的錯誤，持續幾十秒，如果確定已經成功完成的話可以直接掐斷退出
其他的全局參數可以參考http://www.broadinstitute.org/gsa/gatkdocs/release/org_broadinstitute_sting_gatk_CommandLineGATK.html，有的用到的不太多，這邊就暫時跳過了。
注：命令行里的斜杠(\)是linux底下用來斷行顯示的，通過\把一行比較長的命令行斷成幾行寫可以方便代碼的閱讀，執行的時候這個地方不會被當作行尾而實現把后面的都連起來的目的，注意\后面不能再有空格或其他的字符，否則就起不到斷行的效果；另外Windows下無效<img title="GATK使用簡介 Part1/2" alt="GATK使用簡介 Part1/2" src="http://www.sinaimg.cn/uc/myshow/blog/misc/gif/E___6726EN00SIGG.gif" type="face" real_src="http://www.sinaimg.cn/uc/myshow/blog/misc/gif/E___6726EN00SIGG.gif" />
下面簡單的介紹一下GATK里面一些常用的子程序(工具)。
==============================UnifiedGenotyper============================
這個可以說是GATK里面用到最多的一個工具，這個工具就是我們call SNP (單核苷酸多態性)或者Indel (插入缺失)的時候用的，最基本的用法如下：
java -jar GenomeAnalysisTK.jar \
-R resources/Homo_sapiens_assembly18.fasta \
-T UnifiedGenotyper \
-I sample1.bam \
-o snps.raw.vcf
前幾個參數之前都講過了，然后簡要解釋一下其他的一些參數：
--out / -o (VCFWriter with default value stdout)
輸出文件，輸出文件的格式為vcf格式，vcf格式的具體信息可以參閱http://www.1000genomes.org/wiki/Analysis/Variant Call Format/vcf-variant-call-format-version-41，這個參數是必須提供的
--output_mode / -out_mode (OUTPUT_MODE with default value EMIT_VARIANTS_ONLY)
輸出模式，總共有3種
EMIT_VARIANTS_ONLY 只輸出variants位點(SNPs和INDELs)
EMIT_ALL_CONFIDENT_SITES 除了variants位點以外與reference一致的可靠位點也輸出
EMIT_ALL_SITES 只要是能call的到的位點都輸出(應該是除了沒覆蓋到的位點，其余位點都輸出)，這個模式下會得到最多的結果
-stand_emit_conf / --standard_min_confidence_threshold_for_emitting (double with default value 30.0)
只有QUAL值大於這個值的才會輸出
-stand_call_conf / --standard_min_confidence_threshold_for_calling (double with default value 30.0)
最低可信值，QUAL值低於這個值而大於上面stand_emit_conf這個值的在結果里面顯示的時候會被標記為LowQual
這兩個參數默認的都是30，就是只有QUAL值大於30的才會輸出到結果中，一般認為對於高覆蓋度的數據(>10x)，QUAL值在30以上就可以認為比較可信
--genotype_likelihoods_model / -glm (Model with default value SNP)
calling模式，默認值為SNP，只輸出SNP的結果，INDEL模式用來輸出INDEL結果，BOTH則會把SNP和INDEL結果都輸出出來
--comp / -comp (List[RodBinding[VariantContext]] with default value [])
用這個參數可以提供一個(或多個，設多個時同樣是設這個參數多次，每次一個文件)已有的vcf結果文件，用於與結果進行比較，如果提供的vcf里面的variants結果與call出來的結果一致，則會在后面加上comp1，[comp2…]，或者用戶可以自己指定每個文件用的標識名稱，例如:
-comp:sample1 sample1.vcf \
-comp:sample2 sample2.vcf \
…
--dbsnp / -D (RodBinding[VariantContext] with default value none)
dbSNP文件，這個參數實際跟上面-comp的參數差不多，不過dbSNP數據在GATK里面有特殊的用途，后面會再提到
=========================== VariantFiltration==========================
用法：
Bash語言: 高亮代碼由發芽網提供
java -Xmx2g -jar GenomeAnalysisTK.jar \
-R ref.fasta \
-T VariantFiltration \
-o output.vcf \
--variant input.vcf \
--filterExpression "AB < 0.2 || MQ0 > 50"\
--filterName "Nov09filters"\
--mask mask.vcf \
--maskName InDel
根據用戶自定義的標准來篩選VCF文件，常用參數如下：
--variant / -V (required RodBinding[VariantContext])
輸入的VCF文件，GATK對文件格式的完整性要求高一點，所以VCF文件格式不符合要求的可能會導致報錯退出。
--out / -o (VCFWriter with default value stdout)
輸出的文件名，結果輸出還是VCF格式的，然后篩選掉的結果不會直接被刪掉，而是用標識給標記出來，就是VCF文件里面FILTER這一列的信息，而VCF開頭也會給出這些標識用的是哪些篩選標准，而所有通過篩選的位點這一列會以PASS標出。
--clusterSize / -cluster (Integer with default value 3)
--clusterWindowSize / -window (Integer with default value 0)
這兩個參數一起用可以設定指定長度內的variants數量的最大值，例如-window設成10，-cluster設成3，表示10bp以內的variants的數量不應當超過3個，如果超過了則會用SnpCluster標示出來。一般這種聚集的SNP可以被認為是不可靠的(10bp里面超過3個(Bowen et al., 2011))。
--filterExpression / -filter (ArrayList[String] with default value [])
--filterName / -filterName (ArrayList[String] with default value [])
這兩個參數也是一起使用的，可以一個第一個參數設置篩選用的條件表達式，第二個參數設置這個篩選在結果里面標識的名稱，兩個參數設置的時候必須是成對的，然后可以設置多次，例如上面的例子中：
--filterExpression "AB < 0.2 || MQ0 > 50" \
--filterName "Nov09filters" \
也可以寫成
--filterExpression "AB < 0.2" \
--filterName "ABfilter" \
--filterExpression "MQ0 > 50" \
--filterName "MQ0filter" \
--invalidatePreviousFilters (boolean with default value false)
如果輸入的VCF文件是已經篩選過的，而且里面篩選掉的內容是已經不要的了，可以通過設置這個參數來去掉這些已經篩選掉的位點；否則這些位點會被再次篩選一遍，然后保留原來的篩選標識不變，后面可能會加上新的篩選標識。
--mask (RodBinding[Feature] with default value none)
--maskName / -maskName (String with default value Mask)
--mask參數用來篩選與某個指定文件中記錄一致的位點，例如篩選repeat區域的位點的時候就可以使用這個參數，這邊給的是一個文件，支持的格式很多，所以適用的范圍其實很廣，這邊暫時先簡單的介紹一下；--maskName同樣是來設置被FILTER掉的位點的標識的。
--missingValuesInExpressionsShouldEvaluateAsFailing (Boolean with default value false)
最后在講一下這個參數，當篩選標准比較多的時候，可能有一些位點沒有篩選條件當中的一條或幾條，例如下面的這個表達式
QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || HaplotypeScore > 13.0
並不一定所有位點都有這些信息，這種情況下GATK運行的時候會報很多WARNING信息，用這個參數可以把這些缺少某些FLAG的位點也給標記成沒有通過篩選的。
=========================== DepthOfCoverage==========================
這個也是我比較喜歡用的一個工具，可以對mapping結果進行一些列的統計得到覆蓋度相關的一些信息，用法：
Bash語言: 高亮代碼由發芽網提供
java -Xmx2g -jar GenomeAnalysisTK.jar \
-R ref.fasta \
-T DepthOfCoverage \
-o file_name_base \
-I input_bams.list
[-geneList refSeq.sorted.txt]\
[-pt readgroup]\
[-ct 4 -ct 6 -ct 10]\
[-L my_capture_genes.interval_list]
*[]里面的表示可選的參數
因為一般粗略統計的時候，只需要_summary和_statistics這兩個文件就能獲得比較有用的一些信息，而對於特定區段或者gene區間等等的，只要再次基礎上指定某段區域或者提供gene位置的一些信息即可，下面就僅僅介紹和整個基因組統計相關的一些參數：
--out / -o (Map[DoCOutputType,PrintStream] with default value None)
指定輸出文件的文件名，默認情況，這個程序運行完成以后除了生成這個文件以外，還會生成一些列其他的文件，前面部分的名稱還是這個，后面會加上_summary, _statistics等等來表示對不同信息的統計結果，每個文件的內容說明如下：
- no suffix:這個就是沒有加任何后綴的文件，每個位點的覆蓋度，這個文件一般會很大，而且很多時候可能根本不需要這么詳細的結果，可以通過后面的參數設置取消輸出這個結果
- _summary:總的統計信息，包括總的鹼基數，基因組的平均覆蓋度，大於某個quality值的鹼基數所占的比例等等
- _statistics:某個覆蓋度范圍之間鹼基的數目，這邊是一個直方圖的統計結果，可以根據這個做出一張覆蓋度的分布圖，很直觀的顯示出覆蓋度的信息。
- _interval_summary:如果設置了特定區段而不是整個基因組的話，相關結果會在這個文件里面
- _interval_statistics:同上，也是特定區段的統計結果
- _gene_summary: gene區間的統計結果
- _gene_statistics:同上，這邊的幾個文件可以方便對某些比較關心的區段，或者是gene區間進行一些統計
- _cumulative_coverage_counts:大於某個覆蓋度值的
- _cumulative_coverage_proportions:同上，一些比例信息
--partitionType / -pt (Set[Partition] with default value [sample])
DepthOfCoverage可以對多個BAM文件進行統計，這個時候可以將多個文件按照不同的方式進行組合統計，所有的組合以BAM文件頭里面的ReadGroup這一行即開頭是@RG的這一行設定的值為准，可以以sample, readgroup, library…等等字段進行組合，例如-pt設成sample的話，sample名稱相同的BAM文件會被當作一個sample來處理，結果統計的就是一個sample的信息，但是如果這幾個文件的readgroup值不一樣的話，根據readgroup進行統計得到的結果就是對每個readgroup而言的。
--includeRefNSites (boolean with default value false)
mapping所用的reference基因組序列可能包含N，如果這個參數設成true的話，它會根據N附近reads的mapping情況來對這段區域進行一個統計
--nBins (int with default value 499)
DepthOfCoverage進行統計的時候是按照不重疊的window (這邊稱為一個bin)來進行統計的，這邊就是指這個bin的個數，也就是最后得到的直方圖里面看到的柱形的個數。
--start (int with default value 1)
--stop (int with default value 500)
這兩個值用來設定開頭和結尾的bin里面應當包含覆蓋度是多少的鹼基，--start設成5就表示第一個bin統計的是所有depth≤5的鹼基的個數，而—stop設成100就表示最后一個bin統計的是所有depth≥100的鹼基的個數，注意這邊stop的值和start的值之間的間隔應該剛好等於上面nBins設置的值，即stop – start = nBins
--summaryCoverageThreshold / -ct (int[] with default value [15])
這個值表示統計depth≥這個值的鹼基所占的比例，可以設定多個值，例如：
-ct 0 -ct 1 -ct 5 -ct 10 -ct 15，則會分別統計depth≥0, 1, 5, 10, 15的鹼基所占的百分比，結果輸出在_summary文件里
--omitDepthOutputAtEachBase / -omitBaseOutput (boolean with default value false)
--omitIntervalStatistics / -omitIntervals (boolean with default value false)
--omitLocusTable / -omitLocusTable (boolean with default value false)
實際用的時候很多太詳細的結果用處可能不大(例如對每個base的統計DepthOutputAtEachBase)，可以通過這幾個參數來取消，同時可以加快程序的運行速度
所以一個快速統計mapping結果的例子大致如下：
Bash語言: 高亮代碼由發芽網提供
java -Xmx4g -jar GenomeAnalysisTK.jar \
-R ref_genome.fasta \
-T DepthOfCoverage \
--omitDepthOutputAtEachBase \
--omitIntervalStatistics \
--omitLocusTable \
-ct 0 -ct 1 -ct 5 -ct 10 -ct 15 -ct 20 -ct 25 -ct 30 \
--nBins 99 --start 1 --stop 100 \
-I sample.bam \
-o sample.depth