GATK使用方法詳解-plob最詳盡說明書

-.z.GATK使用方法詳解一、使用GATK前須知事項：〔1〕對GATK的測試主要使用的是人類全基因組和外顯子組的測序數(shù)據(jù)，而且全部是基于illumina數(shù)據(jù)格式，目前還沒有提供其他格式文件〔如IonTorrent〕或者實驗設(shè)計〔RNA-Seq〕的分析方法?！?〕GATK是一個應(yīng)用于前沿科學(xué)研究的軟件，不斷在更新和修正，因此，在使用GATK進展變異檢測時，最好是下載最新的版本，目前的版本是2.8.1〔2014-02-25〕。下載：./gatk/download?！?〕在GATK使用過程中〔見下面圖〕，有些步驟需要用到變異信息，對于這些變異，GATK只提供了人類的變異信息，可以在GATK的FTP站點下載〔GATKresourcebundle〕。如果要研究的不是人類基因組，需要自行構(gòu)建變異，GATK提供了詳細的構(gòu)建方法?！?〕GATK在進展BQSR和VQSR的過程中會使用到R軟件繪制一些圖，因此，在運行GATK之前最好先檢查一下是否正確安裝了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果畫圖時出現(xiàn)錯誤，會提示需要安裝的包的名稱。二、GATK的使用流程GATK最正確使用方案：共3大步驟，即:原始數(shù)據(jù)的處理-->變異檢測-->初步分析。原始數(shù)據(jù)的處理1.對原始下機fastq文件進展過濾和比對〔mapping〕對于Illumina下機數(shù)據(jù)推薦使用bwa進展mapping。Bwa比對步驟大致如下：〔1〕對參考基因組構(gòu)建索引：例子：bwainde*-abwtswhg19.fa。構(gòu)建索引時需要注意的問題：bwa構(gòu)建索引有兩種算法，兩種算法都是基于BWT的，這兩種算法通過參數(shù)-ais和-abwtsw進展選擇。其中-abwtsw對于短的參考序列是不工作的，必須要大于等于10Mb；-ais是默認(rèn)參數(shù)，這個參數(shù)不適用于大的參考序列，必須要小于等于2G?！?〕尋找輸入reads文件的SA坐標(biāo)。對于pairend數(shù)據(jù)，每個reads文件單獨做運算，singleend數(shù)據(jù)就不用說了，只有一個文件。pairend：singleend：主要參數(shù)說明：-oint：允許出現(xiàn)的最大gap數(shù)。-eint：每個gap允許的最大長度。-dint：不允許在3’端出現(xiàn)大于多少bp的deletion。-iint：不允許在reads兩端出現(xiàn)大于多少bp的indel。-lint：Read前多少個堿基作為seed，如果設(shè)置的seed大于read長度，將無法繼續(xù)，最好設(shè)置在25-35，與-k2配合使用。-kint：在seed中的最大編輯距離，使用默認(rèn)2，與-l配合使用。-tint：要使用的線程數(shù)。-Rint：此參數(shù)只應(yīng)用于pairend中，當(dāng)沒有出現(xiàn)大于此值的最正確比對結(jié)果時，將會降低標(biāo)準(zhǔn)再次進展比對。增加這個值可以提高配比照對的準(zhǔn)確率，但是同時會消耗更長的時間，默認(rèn)是32。-Iint：表示輸入的文件格式為Illumina1.3+數(shù)據(jù)格式。-Bint：設(shè)置標(biāo)記序列。從5’端開場多少個堿基作為標(biāo)記序列，當(dāng)-B為正值時，在比對之前會將每個read的標(biāo)記序列剪切，并將此標(biāo)記序列表示在BCSAM標(biāo)簽里，對于pairend數(shù)據(jù)，兩端的標(biāo)記序列會被連接。-b：指定輸入格式為bam格式。這是一個很奇怪的功能，就是對其它軟件的bam文件進展重新比對的意思〔3〕生成sam格式的比對文件。如果一條read比對到多個位置，會隨機選擇一種。例子：singleend：參數(shù)：-nint：如果reads比對次數(shù)超過多少次，就不在*A標(biāo)簽顯示。-rstr：定義頭文件?！甊G\tID:foo\tSM:bar’，如果在此步驟不進展頭文件定義，在后續(xù)GATK分析中還是需要重新增加頭文件。pairend：bwasampe-a500read1.fq.gz.sairead2.fq.gz.sairead1.fq.gzread2.fq.gz>read.sam參數(shù)：-aint：最大插入片段大小。-oint：pairend兩reads中其中之一所允許配對的最大次數(shù)，超過該次數(shù)，將被視為singleend。降低這個參數(shù)，可以加快運算速度，對于少于30bp的read，建議降低-o值。-rstr：定義頭文件。同singleend。-nint：每對reads輸出到結(jié)果中的最多比對數(shù)。對于最后得到的sam文件，將比對上的結(jié)果提取出來〔awk即可處理〕，即可直接用于GATK的分析。注意：由于GATK在下游的snp-calling時，是按染色體進展call-snp的。因此，在準(zhǔn)備原始sam文件時，可以先按染色體將文件分開，這樣會提高運行速度。但是當(dāng)數(shù)據(jù)量缺乏時，可能會影響后續(xù)的VQSR分析，這是需要注意的。2.對sam文件進展進展重新排序〔reorder〕由BWA生成的sam文件時按字典式排序法進展的排序〔le*icographically〕進展排序的〔chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chr*，chrY〕，但是GATK在進展callsnp的時候是按照染色體組型〔karyotypic〕進展的〔chrM，chr1，chr2…chr22，chr*，chrY〕，因此要對原始sam文件進展reorder?？梢允褂胮icard-tools中的ReorderSam完成。eg.java-jarpicard-tools-1.96/ReorderSam.jarI=hg19.samREFERENCE=hg19.fa注意：1)這一步的頭文件可以人工加上，同時要確保頭文件中有的序號在下面序列中也有對應(yīng)的。雖然在GATK上的說明chrM可以在最前也可以在最后，但是當(dāng)把chrM放在最后時可能會出錯。2)在進展排序之前，要先構(gòu)建參考序列的索引。e.g.samtoolsfaid*hg19.fa。最后生成的索引文件：hg19.fa.fai。3)如果在上一步想把大文件切分成小文件的時候，頭文件可以自己手工加上，之后運行這一步就好了。3.將sam文件轉(zhuǎn)換成bam文件〔bam是二進制文件，運算速度快〕這一步可使用samtoolsview完成。e.g.samtoolsview-bShg19.reorder_00.sam-ohg19.sam_01.bam

4.對bam文件進展sort排序處理這一步是將sam文件中同一染色體對應(yīng)的條目按照坐標(biāo)順序從小到大進展排序。可以使用picard-tools中SortSam完成。e.g.java-jarpicard-tools-1.96/SortSam.jarSORT_ORDER=coordinate5.對bam文件進展加頭〔head〕處理GATK2.0以上版本將不再支持無頭文件的變異檢測。加頭這一步可以在BWA比對的時候進展，通過-r參數(shù)的選擇可以完成。如果在BWA比對期間沒有選擇-r參數(shù)，可以增加這一步驟。可使用picard-tools中AddOrReplaceReadGroups完成。e.g.java-jarpicard-tools-1.96/AddOrReplaceReadGroups.jarID=hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19IDstr：輸入reads集ID號；LB：read集文庫名；PL：測序平臺〔illunima或solid〕；PU：測序平臺下級單位名稱〔run的名稱〕；SM：樣本名稱。注意：這一步盡量不要手動加頭，本人嘗試過屢次手工加頭，雖然看起來與軟件加的頭是一樣的，但是程序卻無法運行。6.Merge如果一個樣本分為多個lane進展測序，則在進展下一步之前可以將每個lane的bam文件合并。e.g.java-jarpicard-tools-1.70/MergeSamFiles.jarINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bam……INPUT=lane8.bamOUTPUT=sample.bam7.DuplicatesMarking在制備文庫的過程中，由于PCR擴增過程中會存在一些偏差，也就是說有的序列會被過量擴增。這樣，在比對的時候，這些過量擴增出來的完全一樣的序列就會比對到基因組的一樣位置。而這些過量擴增的reads并不是基因組自身固有序列，不能作為變異檢測的證據(jù)，因此，要盡量去除這些由PCR擴增所形成的duplicates，這一步可以使用picard-tools來完成。去重復(fù)的過程是給這些序列設(shè)置一個flag以標(biāo)志它們，方便GATK的識別。還可以設(shè)置REMOVE_DUPLICATES=true來丟棄duplicated序列。對于是否選擇標(biāo)記或者刪除，對結(jié)果應(yīng)該沒有什么影響，GATK官方流程里面給出的例子是僅做標(biāo)記不刪除。這里定義的重復(fù)序列是這樣的：如果兩條reads具有一樣的長度而且比對到了基因組的同一位置，則就認(rèn)為這樣的reads是由PCR擴增而來，就會被GATK標(biāo)記。e.g.java-jarpicard-tools-1.96/MarkDuplicates.jarREMOVE_DUPLICATES=falseMA*_FILE_HANDLES_FOR_READ_ENDS_MAP=8000注意：dedup這一步只要在library層面上進展就可以了，例如一個sample如果建了多個庫的話，對每個庫進展dedup即可，不需要把所有庫合成一個sample再進展dedup操作。其實并不能準(zhǔn)確的定義被mask的reads到底是不是duplicates，重復(fù)序列的程度與測序深度和文庫類型都有關(guān)系。最主要目的就是盡量減小文庫構(gòu)建時引入文庫的PCRbias。8.要對上一步得到的結(jié)果生成索引文件可以用samtools完成，生成的索引后綴是bai。e.g.9.Localrealignmentaroundindels這一步的目的就是將比對到indel附近的reads進展局部重新比對，將比對的錯誤率降到最低。一般來說，絕大局部需要進展重新比對的基因組區(qū)域，都是因為插入/缺失的存在，因為在indel附近的比對會出現(xiàn)大量的堿基錯配，這些堿基的錯配很容易被誤認(rèn)為SNP。還有，在比對過程中，比對算法對于每一條read的處理都是獨立的，不可能同時把多條reads與參考基因組比對來排錯。因此，即使有一些reads能夠正確的比對到indel，但那些恰恰比對到indel開場或者完畢位置的read也會有很高的比對錯誤率，這都是需要重新比對的。Localrealignment就是將由indel導(dǎo)致錯配的區(qū)域進展重新比對，將indel附近的比對錯誤率降到最低。主要分為兩步：第一步，通過運行RealignerTargetCreator來確定要進展重新比對的區(qū)域。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TRealignerTargetCreator-known-known參數(shù)說明：-R：參考基因組；-T：選擇的GATK工具；-I：輸入上一步所得bam文件；-o：輸出的需要重新比對的基因組區(qū)域結(jié)果；-ma*Interval：允許進展重新比對的基因組區(qū)域的最大值，不能太大，太大消耗會很長時間，默認(rèn)值500；-known：的可靠的indel位點，指定的可靠的indel位點，重比對將主要圍繞這些位點進展，對于人類基因組數(shù)據(jù)而言，可以直接指定GATKresourcebundle里面的indel文件〔必須是vcf文件〕。對于knownsites的選擇很重要，GATK中每一個用到knownsites的工具對于knownsites的使用都是不一樣的，但是所有的都有一個共同目的，那就是分辨真實的變異位點和不可信的變異位點。如果不提供這些knownsites的話，這些統(tǒng)計工具就會產(chǎn)生偏差，最后會嚴(yán)重影響結(jié)果的可信度。在這些需要知道knownsites的工具里面，只有UnifiedGenotyper和HaplotypeCaller對knownsites沒有太嚴(yán)格的要求。如果你所研究的對象是人類基因組的話，那就簡單多了，因為GATK上對如何使用人類基因組的knownsites做出了詳細的說明，具體的選擇方法如下表，這些文件都可以在GATKresourcebundle中下載。TooldbSNP129dbSNP>132Millsindels1KGindelsHapMapOmniRealignerTargetCreator**IndelRealigner**BaseRecalibrator***(UnifiedGenotyper/HaplotypeCaller)*VariantRecalibrator****VariantEval*但是如果你要研究的不是人類基因組的話，那就有點麻煩了，./gatk/guide/article"id=1243，這個上是做非人類基因組時，大家分享的經(jīng)歷，可以參考一下。這個knownsites如果實在沒有的話，也是可以自己構(gòu)建的：首先，先使用沒有經(jīng)過矯正的數(shù)據(jù)進展一輪SNPcalling；然后，挑選最可信的SNP位點進展BQSR分析；最后，在使用這些經(jīng)過BQSR的數(shù)據(jù)進展一次真正的SNPcalling。這幾步可能要重復(fù)好屢次才能得到可靠的結(jié)果。第二步，通過運行IndelRealigner在這些區(qū)域進展重新比對。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TIndelRealigner運行完畢后，生成的hg19.dedup.realn_07.bam即為最后重比對后的文件。注意：1.第一步和第二步中使用的輸入文件〔bam文件〕、參考基因組和indel文件必須是一樣的文件。2.當(dāng)在一樣的基因組區(qū)域發(fā)現(xiàn)多個indel存在時，這個工具會從其中選擇一個最有可能存在比對錯誤的indel進展重新比對，剩余的其他indel不予考慮。3.對于454下機數(shù)據(jù)，本工具不支持。此外，這一步還會忽略bwa比對中質(zhì)量值為0的read以及在CIGAR信息中存在連續(xù)indel的reads。10.Basequalityscorerecalibration這一步是對bam文件里reads的堿基質(zhì)量值進展重新校正，使最后輸出的bam文件中reads中堿基的質(zhì)量值能夠更加接近真實的與參考基因組之間錯配的概率。這一步適用于多種數(shù)據(jù)類型，包括illunima、solid、454、CG等數(shù)據(jù)格式。在GATK2.0以上版本中還可以對indel的質(zhì)量值進展校正，這一步對indelcalling非常有幫助舉例說明，在reads堿基質(zhì)量值被校正之前，我們要保存質(zhì)量值在Q25以上的堿基，但是實際上質(zhì)量值在Q25的這些堿基的錯誤率在1%，也就是說質(zhì)量值只有Q20，這樣就會對后續(xù)的變異檢測的可信度造成影響。還有，在邊合成邊測序的測序過程中，在reads末端堿基的錯誤率往往要比起始部位更高。另外，AC的質(zhì)量值往往要低于TG。BQSR的就是要對這些質(zhì)量值進展校正。BQSR主要有三步：第一步：利用工具BaseRecalibrator，根據(jù)一些knownsites，生成一個校正質(zhì)量值所需要的數(shù)據(jù)文件，GATK以“.grp〞為后綴命名。e.g.java-jarGenomeAnalysisTK.jar-TBaseRecalibrator-Rhg19.fa第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp來生成校正后的數(shù)據(jù)文件，也是以“.grp〞命名，這一步主要是為了與校正之前的數(shù)據(jù)進展比擬，最后生成堿基質(zhì)量值校正前后的比擬圖，如果不想生成最后BQSR比擬圖，這一步可以省略。e.g.java-jarGenomeAnalysisTK.jar-TBaseRecalibrator-Rhg19.fa-oGATK第三步：利用工具PrintReads將經(jīng)過質(zhì)量值校正的數(shù)據(jù)輸出到新的bam文件中，用于后續(xù)的變異檢測。e.g.java-jarGenomeAnalysisTK.jar-TPrintReads-Rhg19.fa主要參數(shù)說明：-bqsrBAQGOP：BQSRBAQgapopen罰值，默認(rèn)值是40，如果是對全基因組數(shù)據(jù)進展BQSR分析，設(shè)置為30會更好。-lqt：在計算過程中，該算法所能考慮的reads兩端的最小質(zhì)量值。如果質(zhì)量值小于該值，計算過程中將不予考慮，默認(rèn)值是2。注意：〔1〕當(dāng)bam文件中的reads數(shù)量過少時，BQSR可能不會正常工作，GATK建議base數(shù)量要大于100M才能得到比擬好的結(jié)果。〔2〕除非你所研究的樣本所得到的reads數(shù)實在太少，或者比對結(jié)果中的mismatch根本上都是實際存在的變異，否則必須要進展BQSR這一步。對于人類基因組，即使有了dbSNP和千人基因組的數(shù)據(jù)，還有很多mismatch是錯誤的。因此，這一步能做一定要做。11.分析和評估BQSR結(jié)果這一步會生成評估前后堿基質(zhì)量值的比擬結(jié)果，可以選擇使用圖片和表格的形式展示。e.g.java-jarGenomeAnalysisTK.jar-TAnalyzeCovariates-Rhg19.fa參數(shù)解釋：-before：基于原始比對結(jié)果生成的第一次校對表格。-after：基于第一次校對表格生成的第二次校對表格。-plots：評估BQSR結(jié)果的報告文件。-csv：生成報告中圖標(biāo)所需要的所有數(shù)據(jù)。

12.Reducebamfile這一步是使用ReduceReads這個工具將bam文件進展壓縮，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有進展變異檢測所需要的信息。這樣不僅能夠節(jié)省存儲空間，也方便后續(xù)變異檢測過程中對數(shù)據(jù)的處理。e.g.java-jarGenomeAnalysisTK.jar-TReduceReads-Rhg19.fa到此為止，GATK流程中的第一大步驟就完畢了，完成了variantscalling所需要的所有準(zhǔn)備工作，生成了用于下一步變異檢測的bam文件。變異檢測1.VariantCallingGATK在這一步里面提供了兩個工具進展變異檢測——UnifiedGenotyper和HaplotypeCaller。其中HaplotypeCaller一直還在開發(fā)之中，包括生成的結(jié)果以及計算模型和命令行參數(shù)一直在變動，因此，目前使用比擬多的還是UnifiedGenotyper。此外，HaplotypeCaller不支持Reduce之后的bam文件，因此，中選擇使用HaplotypeCaller進展變異檢測時，不需要進展Reducereads。UnifiedGenotyper是集合多種變異檢測方法而成的一種VariantsCaller，既可以用于單個樣本的變異檢測，也可以用于群體的變異檢測。UnifiedGenotyper使用貝葉斯最大似然模型，同時估計基因型和基因頻率，最后對每一個樣本的每一個變異位點和基因型都會給出一個準(zhǔn)確的后驗概率。e.g.java-jarGenomeAnalysisTK.jar-glmBOTH-lINFO-Rhg19.fa-TUnifiedGenotyper-stand_call_conf10-stand_emit_conf30上述命令將對輸入的bam文件中的所有樣本進展變異檢測，最后生成一個vcf文件，vcf文件中會包含所有樣本的變異位點和基因型信息。但是現(xiàn)在所得到的結(jié)果是最原始的、沒有經(jīng)過任何過濾和校正的Variants集合。這一步產(chǎn)生的變異位點會有很高的假陽性，尤其是indel，因此，必須要進展進一步的篩選過濾。這一步還可以指定對基因組的*一區(qū)域進展變異檢測，只需要增加一個參數(shù)-L：target_interval.list，格式是bed格式文件。主要參數(shù)解釋：-A：指定一個或者多個注釋信息，最后輸出到vcf文件中。-*A：指定不做哪些注釋，最后不會輸出到vcf文件中。-D：的snp文件。-glm：選擇檢測變異的類型。SNP表示只進展snp檢測；INDEL表示只對indel進展檢測；BOTH表示同時檢測snp和indel。默認(rèn)值是SNP。-hets：雜合度的值，用于計算先驗概率。默認(rèn)值是0.001。-ma*AltAlleles：容許存在的最大altallele的數(shù)目，默認(rèn)值是6。這個參數(shù)要特別注意，不要輕易修改默認(rèn)值，程序設(shè)置的默認(rèn)值幾乎可以滿足所有的分析，如果修改了可能會導(dǎo)致程序無法運行。-mbq：變異檢測時，堿基的最小質(zhì)量值。如果小于這個值，將不會對其進展變異檢測。這個參數(shù)不適用于indel檢測，默認(rèn)值是17。-minIndelt：在做indelcalling的時候，支持一個indel的最少read數(shù)量。也就是說，如果同時有多少條reads同時支持一個候選indel時，軟件才開場進展indelcalling。降低這個值可以增加indelcalling的敏感度，但是會增加消耗的時間和假陽性。-minIndelFrac：在做indelcalling的時候，支持一個indel的reads數(shù)量占比對到該indel位置的所有reads數(shù)量的百分比。也就是說，只有同時滿足-minIndelt和-minIndelFrac兩個參數(shù)條件時，才會進展indelcalling。-onlyEmitSamples：當(dāng)指定這個參數(shù)時，只有指定的樣本的變異檢測結(jié)果會輸出到vcf文件中。-stand_emit_conf：在變異檢測過程中，所容許的最小質(zhì)量值。只有大于等于這個設(shè)定值的變異位點會被輸出到結(jié)果中。-stand_call_conf：在變異檢測過程中，用于區(qū)分低質(zhì)量變異位點和高質(zhì)量變異位點的閾值。只有質(zhì)量值高于這個閾值的位點才會被視為高質(zhì)量的。低于這個質(zhì)量值的變異位點會在輸出結(jié)果中標(biāo)注LowQual。在千人基因組方案第二階段的變異檢測時，利用35*的數(shù)據(jù)進展snpcalling的時候，當(dāng)設(shè)置成50時，有大概10%的假陽性。-dcov：這個參數(shù)用于控制檢測變異數(shù)據(jù)的coverage(*)，4*的數(shù)據(jù)可以設(shè)置為40，大于30*的數(shù)據(jù)可以設(shè)置為200。注意：GATK進展變異檢測的時候，是按照染色體排序順序進展的〔先callchr1，然后chr2，然后chr3…最后chrY〕，并非多條染色體并行檢測的，因此，如果數(shù)據(jù)量比擬大的話，建議分染色體分別進展，對性染色體的變異檢測可以同常染色體方法。大多數(shù)參數(shù)的默認(rèn)值可以滿足大多數(shù)研究的需求，因此，在做變異檢測過程中，如果對參數(shù)意義不是很明確，不建議修改。2.對原始變異檢測結(jié)果進展過濾〔hardfilterandVQSR〕這一步的目的就是對上一步call出來的變異位點進展過濾，去掉不可信的位點。這一步可以有兩種方法，一種是通過GATK的VariantFiltration，另一種是通過GATK的VQSR〔變異位點質(zhì)量值重新校正〕進展過濾。通過GATK上提供的最正確方案可以看出，GATK是推薦使用VASR的，但使用VQSR數(shù)據(jù)量一定要到達要求，數(shù)據(jù)量太小無法使用高斯模型。還有，在使用VAQR時，indel和snp要分別進展。VQSR原理介紹：這個模型是根據(jù)已有的真實變異位點〔人類基因組一般使用HapMap3中的位點，以及這些位點在Omni2.5MSNP芯片中出現(xiàn)的多態(tài)位點〕來訓(xùn)練，最后得到一個訓(xùn)練好的能夠很好的評估真?zhèn)蔚腻e誤評估模型，可以叫他適應(yīng)性錯誤評估模型。這個適應(yīng)性的錯誤評估模型可以應(yīng)用到call出來的原始變異集合中的變異位點和新發(fā)現(xiàn)的變異位點，進而去評估每一個變異位點發(fā)生錯誤的概率，最終會給出一個得分。這個得分最后會被寫入vcf文件的INFO信息里，叫做VQSLOD，就是在訓(xùn)練好的混合高斯模型下，一個位點是真實的概率比上這個位點可能是假陽性的概率的logoddsratio〔對數(shù)差異比〕，因此，可以定性的認(rèn)為，這個值越大就越好。VQSR主要分兩個步驟，這兩個步驟會使用兩個不同的工具：VariantRecalibrator和ApplyRecalibration。i)VariantRecalibratorVariantRecalibrator：通過大量的高質(zhì)量的變異集合的各個注釋〔包括很多種，后面介紹〕的值來創(chuàng)立一個高斯混合模型，然后用于評估所有的變異位點。這個文件最后將生成一個recalibration文件。原理簡單