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

年4月19日GATK使用方法詳解plob最詳盡說明書文檔僅供參考GATK使用方法詳解一、使用GATK前須知事項(xiàng)：（1）對GATK的測試主要使用的是人類全基因組和外顯子組的測序數(shù)據(jù)，而且全部是基于illumina數(shù)據(jù)格式，當(dāng)前還沒有提供其它格式文件（如IonTorrent）或者實(shí)驗(yàn)設(shè)計(jì)（RNA-Seq）的分析方法。（2）GATK是一個(gè)應(yīng)用于前沿科學(xué)研究的軟件，不斷在更新和修正，因此，在使用GATK進(jìn)行變異檢測時(shí)，最好是下載最新的版本，當(dāng)前的版本是2.8.1（-02-25）。下載網(wǎng)站：。（3）在GATK使用過程中（見下面圖），有些步驟需要用到已知變異信息，對于這些已知變異，GATK只提供了人類的已知變異信息，能夠在GATK的FTP站點(diǎn)下載（GATKresourcebundle）。如果要研究的不是人類基因組，需要自行構(gòu)建已知變異，GATK提供了詳細(xì)的構(gòu)建方法。（4）GATK在進(jìn)行BQSR和VQSR的過程中會(huì)使用到R軟件繪制一些圖，因此，在運(yùn)行GATK之前最好先檢查一下是否正確安裝了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果畫圖時(shí)出現(xiàn)錯(cuò)誤，會(huì)提示需要安裝的包的名稱。二、GATK的使用流程GATK最佳使用方案：共3大步驟，即:\o"GATK使用方法詳解（原始數(shù)據(jù)的處理）"原始數(shù)據(jù)的處理-->\o"GATK使用方法詳解（變異檢測）"變異檢測

-->\o"GATK使用方法詳解（初步分析）"初步分析。原始數(shù)據(jù)的處理1.對原始下機(jī)fastq文件進(jìn)行過濾和比對（mapping）對于Illumina下機(jī)數(shù)據(jù)推薦使用bwa進(jìn)行mapping。Bwa比對步驟大致如下：（1）對參考基因組構(gòu)建索引：例子：bwaindex-abwtswhg19.fa。構(gòu)建索引時(shí)需要注意的問題：bwa構(gòu)建索引有兩種算法，兩種算法都是基于BWT的，這兩種算法經(jīng)過參數(shù)-ais和-abwtsw進(jìn)行選擇。其中-abwtsw對于短的參考序列是不工作的，必須要大于等于10Mb；-ais是默認(rèn)參數(shù)，這個(gè)參數(shù)不適用于大的參考序列，必須要小于等于2G。（2）尋找輸入reads文件的SA坐標(biāo)。對于pairend數(shù)據(jù)，每個(gè)reads文件單獨(dú)做運(yùn)算，singleend數(shù)據(jù)就不用說了，只有一個(gè)文件。pairend：bwaalnhg19.faread1.fq.gz-t4-I>read1.fq.gz.saibwaalnhg19.faread2.fq.gz-t4-I>read2.fq.gz.saisingleend：bwaalnhg19.faread.fq.gz-l30-k2-t4-I>read.fq.gz.sai主要參數(shù)說明：-oint：允許出現(xiàn)的最大gap數(shù)。-eint：每個(gè)gap允許的最大長度。-dint：不允許在3’端出現(xiàn)大于多少bp的deletion。-iint：不允許在reads兩端出現(xiàn)大于多少bp的indel。-lint：Read前多少個(gè)堿基作為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é)果時(shí)，將會(huì)降低標(biāo)準(zhǔn)再次進(jìn)行比對。增加這個(gè)值能夠提高配對比正確準(zhǔn)確率，可是同時(shí)會(huì)消耗更長的時(shí)間，默認(rèn)是32。-Iint：表示輸入的文件格式為Illumina1.3+數(shù)據(jù)格式。-Bint：設(shè)置標(biāo)記序列。從5’端開始多少個(gè)堿基作為標(biāo)記序列，當(dāng)-B為正值時(shí)，在比對之前會(huì)將每個(gè)read的標(biāo)記序列剪切，并將此標(biāo)記序列表示在BCSAM標(biāo)簽里，對于pairend數(shù)據(jù)，兩端的標(biāo)記序列會(huì)被連接。-b：指定輸入格式為bam格式。這是一個(gè)很奇怪的功能，就是對其它軟件的bam文件進(jìn)行重新比正確意思bwaalnhg19.faread.bam>read.fq.gz.sai（3）生成sam格式的比對文件。如果一條read比對到多個(gè)位置，會(huì)隨機(jī)選擇一種。例子：singleend：bwasamsehg19.faread.fq.gz.sairead.fq.gz>read.fq.gz.sam參數(shù)：-nint：如果reads比對次數(shù)超過多少次，就不在XA標(biāo)簽顯示。-rstr：定義頭文件。‘@RG\tID:foo\tSM:bar’，如果在此步驟不進(jìn)行頭文件定義，在后續(xù)\o"ViewallpostsinGATK"GATK分析中還是需要重新增加頭文件。pairend：bwasampe-a500read1.fq.gz.sairead2.fq.gz.sairead1.fq.gzread2.fq.gz>read.sam參數(shù)：-aint：最大插入片段大小。-oint：pairend兩reads中其中之一所允許配正確最大次數(shù)，超過該次數(shù)，將被視為singleend。降低這個(gè)參數(shù)，能夠加快運(yùn)算速度，對于少于30bp的read，建議降低-o值。-rstr：定義頭文件。同singleend。-nint：每對reads輸出到結(jié)果中的最多比對數(shù)。對于最后得到的sam文件，將比對上的結(jié)果提取出來（awk即可處理），即可直接用于\o"ViewallpostsinGATK"GATK的分析。注意：由于\o"ViewallpostsinGATK"GATK在下游的snp-calling時(shí)，是按染色體進(jìn)行call-snp的。因此，在準(zhǔn)備原始sam文件時(shí)，能夠先按染色體將文件分開，這樣會(huì)提高運(yùn)行速度。可是當(dāng)數(shù)據(jù)量不足時(shí)，可能會(huì)影響后續(xù)的VQSR分析，這是需要注意的。2.對sam文件進(jìn)行進(jìn)行重新排序（reorder）由BWA生成的sam文件時(shí)按字典式排序法進(jìn)行的排序（lexicographically）進(jìn)行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），可是\o"ViewallpostsinGATK"GATK在進(jìn)行callsnp的時(shí)候是按照染色體組型（karyotypic）進(jìn)行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要對原始sam文件進(jìn)行reorder。能夠使用picard-tools中的ReorderSam完成。eg.java-jarpicard-tools-1.96/ReorderSam.jarI=hg19.samO=hg19.reorder_00.samREFERENCE=hg19.fa注意：1)這一步的頭文件能夠人工加上，同時(shí)要確保頭文件中有的序號在下面序列中也有對應(yīng)的。雖然在\o"ViewallpostsinGATK"GATK網(wǎng)站上的說明chrM能夠在最前也能夠在最后，可是當(dāng)把chrM放在最后時(shí)可能會(huì)出錯(cuò)。2)在進(jìn)行排序之前，要先構(gòu)建參考序列的索引。e.g.samtoolsfaidxhg19.fa。最后生成的索引文件：hg19.fa.fai。3)如果在上一步想把大文件切分成小文件的時(shí)候，頭文件能夠自己手工加上，之后運(yùn)行這一步就好了。

3.將sam文件轉(zhuǎn)換成bam文件（bam是二進(jìn)制文件，運(yùn)算速度快）這一步可使用samtoolsview完成。e.g.samtoolsview-bShg19.reorder_00.sam-ohg19.sam_01.bam

4.對bam文件進(jìn)行sort排序處理這一步是將sam文件中同一染色體對應(yīng)的條目按照坐標(biāo)順序從小到大進(jìn)行排序。能夠使用picard-tools中SortSam完成。e.g.java-jarpicard-tools-1.96/SortSam.jarINPUT=hg19.sam_01.bamOUTPUT=hg19.sam.sort_02.bamSORT_ORDER=coordinate

5.對bam文件進(jìn)行加頭（head）處理\o"ViewallpostsinGATK"GATK2.0以上版本將不再支持無頭文件的變異檢測。加頭這一步能夠在BWA比正確時(shí)候進(jìn)行，經(jīng)過-r參數(shù)的選擇能夠完成。如果在BWA比對期間沒有選擇-r參數(shù)，能夠增加這一步驟。可使用picard-tools中AddOrReplaceReadGroups完成。e.g.java-jarpicard-tools-1.96/AddOrReplaceReadGroups.jarI=hg19.sam.sort_02.bamO=hg19.reorder.sort.addhead_03.bamID=hg19IDLB=hg19IDPL=illuminePU=hg19PUSM=hg19IDstr：輸入reads集ID號；LB：read集文庫名；PL：測序平臺（illunima或solid）；PU：測序平臺下級單位名稱（run的名稱）；SM：樣本名稱。注意：這一步盡量不要手動(dòng)加頭，本人嘗試過多次手工加頭，雖然看起來與軟件加的頭是一樣的，可是程序卻無法運(yùn)行。

6.Merge如果一個(gè)樣本分為多個(gè)lane進(jìn)行測序，那么在進(jìn)行下一步之前能夠?qū)⒚總€(gè)lane的bam文件合并。e.g.java-jarpicard-tools-1.70/MergeSamFiles.jarINPUT=lane1.bamINPUT=lane2.bamINPUT=lane3.bamINPUT=lane4.bam……INPUT=lane8.bamOUTPUT=sample.bam

7.DuplicatesMarking在制備文庫的過程中，由于PCR擴(kuò)增過程中會(huì)存在一些偏差，也就是說有的序列會(huì)被過量擴(kuò)增。這樣，在比正確時(shí)候，這些過量擴(kuò)增出來的完全相同的序列就會(huì)比對到基因組的相同位置。而這些過量擴(kuò)增的reads并不是基因組自身固有序列，不能作為變異檢測的證據(jù)，因此，要盡量去除這些由PCR擴(kuò)增所形成的duplicates，這一步能夠使用picard-tools來完成。去重復(fù)的過程是給這些序列設(shè)置一個(gè)flag以標(biāo)志它們，方便\o"ViewallpostsinGATK"GATK的識別。還能夠設(shè)置REMOVE_DUPLICATES=true來丟棄duplicated序列。對于是否選擇標(biāo)記或者刪除，對結(jié)果應(yīng)該沒有什么影響，\o"ViewallpostsinGATK"GATK官方流程里面給出的例子是僅做標(biāo)記不刪除。這里定義的重復(fù)序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置，那么就認(rèn)為這樣的reads是由PCR擴(kuò)增而來，就會(huì)被\o"ViewallpostsinGATK"GATK標(biāo)記。e.g.java-jarpicard-tools-1.96/MarkDuplicates.jarREMOVE_DUPLICATES=falseMAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000INPUT=hg19.reorder.sort.addhead_03.bamOUTPUT=hg19.reorder.sort.addhead.dedup_04.bamMETRICS_FILE=hg19.reorder.sort.addhead.dedup_04.metrics注意：dedup這一步只要在library層面上進(jìn)行就能夠了，例如一個(gè)sample如果建了多個(gè)庫的話，對每個(gè)庫進(jìn)行dedup即可，不需要把所有庫合成一個(gè)sample再進(jìn)行dedup操作。其實(shí)并不能準(zhǔn)確的定義被mask的reads到底是不是duplicates，重復(fù)序列的程度與測序深度和文庫類型都有關(guān)系。最主要目的就是盡量減小文庫構(gòu)建時(shí)引入文庫的PCRbias。

8.要對上一步得到的結(jié)果生成索引文件能夠用samtools完成，生成的索引后綴是bai。e.g.samtoolsindexhg19.reorder.sort.addhead.dedup_04.bam9.Localrealignmentaroundindels這一步的目的就是將比對到indel附近的reads進(jìn)行局部重新比對，將比正確錯(cuò)誤率降到最低。一般來說，絕大部分需要進(jìn)行重新比正確基因組區(qū)域，都是因?yàn)椴迦?缺失的存在，因?yàn)樵趇ndel附近的比對會(huì)出現(xiàn)大量的堿基錯(cuò)配，這些堿基的錯(cuò)配很容易被誤認(rèn)為SNP。還有，在比對過程中，比對算法對于每一條read的處理都是獨(dú)立的，不可能同時(shí)把多條reads與參考基因組比對來排錯(cuò)。因此，即使有一些reads能夠正確的比對到indel，但那些恰恰比對到indel開始或者結(jié)束位置的read也會(huì)有很高的比對錯(cuò)誤率，這都是需要重新比正確。Localrealignment就是將由indel導(dǎo)致錯(cuò)配的區(qū)域進(jìn)行重新比對，將indel附近的比對錯(cuò)誤率降到最低。主要分為兩步：第一步，經(jīng)過運(yùn)行RealignerTargetCreator來確定要進(jìn)行重新比正確區(qū)域。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TRealignerTargetCreator-Ihg19.reorder.sort.addhead.dedup_04.bam-ohg19.dedup.realn_06.intervals-knownMills_and_1000G_gold_standard.indels.hg19.vcf-known1000G_phase1.indels.hg19.vcf參數(shù)說明：-R：參考基因組；-T：選擇的\o"ViewallpostsinGATK"GATK工具；-I：輸入上一步所得bam文件；-o：輸出的需要重新比正確基因組區(qū)域結(jié)果；-maxInterval：允許進(jìn)行重新比正確基因組區(qū)域的最大值，不能太大，太大耗費(fèi)會(huì)很長時(shí)間，默認(rèn)值500；-known：已知的可靠的indel位點(diǎn)，指定已知的可靠的indel位點(diǎn)，重比對將主要圍繞這些位點(diǎn)進(jìn)行，對于人類基因組數(shù)據(jù)而言，能夠直接指定\o"ViewallpostsinGATK"GATKresourcebundle里面的indel文件（必須是vcf文件）。對于knownsites的選擇很重要，\o"ViewallpostsinGATK"GATK中每一個(gè)用到knownsites的工具對于knownsites的使用都是不一樣的，可是所有的都有一個(gè)共同目的，那就是分辨真實(shí)的變異位點(diǎn)和不可信的變異位點(diǎn)。如果不提供這些knownsites的話，這些統(tǒng)計(jì)工具就會(huì)產(chǎn)生偏差，最后會(huì)嚴(yán)重影響結(jié)果的可信度。在這些需要知道knownsites的工具里面，只有UnifiedGenotyper和HaplotypeCaller對knownsites沒有太嚴(yán)格的要求。如果你所研究的對象是人類基因組的話，那就簡單多了，因?yàn)閈o"ViewallpostsinGATK"GATK網(wǎng)站上對如何使用人類基因組的knownsites做出了詳細(xì)的說明，具體的選擇方法如下表，這些文件都能夠在\o"ViewallpostsinGATK"GATKresourcebundle中下載。TooldbSNP129dbSNP>132Millsindels1KGindelsHapMapOmniRealignerTargetCreatorXXIndelRealignerXXBaseRecalibratorXXX(UnifiedGenotyper/HaplotypeCaller)XVariantRecalibratorXXXXVariantEvalX

可是如果你要研究的不是人類基因組的話，那就有點(diǎn)麻煩了，，這個(gè)網(wǎng)站上是做非人類基因組時(shí)，大家分享的經(jīng)驗(yàn)，能夠參考一下。這個(gè)knownsites如果實(shí)在沒有的話，也是能夠自己構(gòu)建的：首先，先使用沒有經(jīng)過矯正的數(shù)據(jù)進(jìn)行一輪SNPcalling；然后，挑選最可信的SNP位點(diǎn)進(jìn)行BQSR分析；最后，在使用這些經(jīng)過BQSR的數(shù)據(jù)進(jìn)行一次真正的SNPcalling。這幾步可能要重復(fù)好多次才能得到可靠的結(jié)果。第二步，經(jīng)過運(yùn)行IndelRealigner在這些區(qū)域內(nèi)進(jìn)行重新比對。e.g.java-jarGenomeAnalysisTK.jar-Rhg19.fa-TIndelRealigner-targetIntervalshg19.dedup.realn_06.intervals-Ihg19.reorder.sort.addhead.dedup_04.bam-ohg19.dedup.realn_07.bam-knownMills_and_1000G_gold_standard.indels.hg19.vcf-known1000G_phase1.indels.hg19.vcf運(yùn)行結(jié)束后，生成的hg19.dedup.realn_07.bam即為最后重比對后的文件。注意：1.第一步和第二步中使用的輸入文件（bam文件）、參考基因組和已知indel文件必須是相同的文件。2.當(dāng)在相同的基因組區(qū)域發(fā)現(xiàn)多個(gè)indel存在時(shí)，這個(gè)工具會(huì)從其中選擇一個(gè)最有可能存在比對錯(cuò)誤的indel進(jìn)行重新比對，剩余的其它indel不予考慮。3.對于454下機(jī)數(shù)據(jù)，本工具不支持。另外，這一步還會(huì)忽略bwa比對中質(zhì)量值為0的read以及在CIGAR信息中存在連續(xù)indel的reads。10.Basequalityscorerecalibration這一步是對bam文件里reads的堿基質(zhì)量值進(jìn)行重新校正，使最后輸出的bam文件中reads中堿基的質(zhì)量值能夠更加接近真實(shí)的與參考基因組之間錯(cuò)配的概率。這一步適用于多種數(shù)據(jù)類型，包括illunima、solid、454、CG等數(shù)據(jù)格式。在\o"ViewallpostsinGATK"GATK2.0以上版本中還能夠?qū)ndel的質(zhì)量值進(jìn)行校正，這一步對indelcalling非常有幫助舉例說明，在reads堿基質(zhì)量值被校正之前，我們要保留質(zhì)量值在Q25以上的堿基，可是實(shí)際上質(zhì)量值在Q25的這些堿基的錯(cuò)誤率在1%，也就是說質(zhì)量值只有Q20，這樣就會(huì)對后續(xù)的變異檢測的可信度造成影響。還有，在邊合成邊測序的測序過程中，在reads末端堿基的錯(cuò)誤率往往要比起始部位更高。另外，AC的質(zhì)量值往往要低于TG。BQSR的就是要對這些質(zhì)量值進(jìn)行校正。BQSR主要有三步：第一步：利用工具BaseRecalibrator，根據(jù)一些knownsites，生成一個(gè)校正質(zhì)量值所需要的數(shù)據(jù)文件，\o"ViewallpostsinGATK"GATK網(wǎng)站以“.grp”為后綴命名。e.g.java-jarGenomeAnalysisTK.jar-TBaseRecalibrator-Rhg19.fa-IChrALL.100.sam.dedup.realn.07.bam-knownSitesdbsnp_137.hg19.vcf-knownSitesMills_and_1000G_gold_standard.indels.hg19.vcf-knownSites1000G_phase1.indels.hg19.vcf-oChrALL.100.sam.recal.08-1.grp第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp來生成校正后的數(shù)據(jù)文件，也是以“.grp”命名，這一步主要是為了與校正之前的數(shù)據(jù)進(jìn)行比較，最后生成堿基質(zhì)量值校正前后的比較圖，如果不想生成最后BQSR比較圖，這一步能夠省略。e.g.java-jarGenomeAnalysisTK.jar-TBaseRecalibrator-Rhg19.fa-IChrALL.100.sam.dedup.realn.07.bam-BQSRChrALL.100.sam.recal.08-1.grp-o\o"ViewallpostsinGATK"GATK/hg19.recal.08-2.grp-knownSitesdbsnp_137.hg19.vcf-knownSitesMills_and_1000G_gold_standard.indels.hg19.vcf-knownSites1000G_phase1.indels.hg19.vcf第三步：利用工具PrintReads將經(jīng)過質(zhì)量值校正的數(shù)據(jù)輸出到新的bam文件中，用于后續(xù)的變異檢測。e.g.java-jarGenomeAnalysisTK.jar-TPrintReads-Rhg19.fa-IChrALL.100.sam.dedup.realn.07.bam-BQSRChrALL.100.sam.recal.08-1.grp-oChrALL.100.sam.recal.08-3.grp.bam主要參數(shù)說明：-bqsrBAQGOP：BQSRBAQgapopen罰值，默認(rèn)值是40，如果是對全基因組數(shù)據(jù)進(jìn)行BQSR分析，設(shè)置為30會(huì)更好。-lqt：在計(jì)算過程中，該算法所能考慮的reads兩端的最小質(zhì)量值。如果質(zhì)量值小于該值，計(jì)算過程中將不予考慮，默認(rèn)值是2。注意：（1）當(dāng)bam文件中的reads數(shù)量過少時(shí)，BQSR可能不會(huì)正常工作，\o"ViewallpostsinGATK"GATK網(wǎng)站建議base數(shù)量要大于100M才能得到比較好的結(jié)果。（2）除非你所研究的樣本所得到的reads數(shù)實(shí)在太少，或者比對結(jié)果中的mismatch基本上都是實(shí)際存在的變異，否則必須要進(jìn)行BQSR這一步。對于人類基因組，即使有了dbSNP和千人基因組的數(shù)據(jù)，還有很多mismatch是錯(cuò)誤的。因此，這一步能做一定要做。

11.分析和評估BQSR結(jié)果這一步會(huì)生成評估前后堿基質(zhì)量值的比較結(jié)果，能夠選擇使用圖片和表格的形式展示。e.g.java-jarGenomeAnalysisTK.jar-TAnalyzeCovariates-Rhg19.fa-beforeChrALL.100.sam.recal.08-1.grp-afterChrALL.100.sam.recal.08-2.grp-csvChrALL.100.sam.recal.grp.09.csv-plotsChrALL.100.sam.recal.grp.09.pdf參數(shù)解釋：-before：基于原始比對結(jié)果生成的第一次校對表格。-after：基于第一次校對表格生成的第二次校對表格。-plots：評估BQSR結(jié)果的報(bào)告文件。-csv：生成報(bào)告中圖標(biāo)所需要的所有數(shù)據(jù)。

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