GATK艰躺，全稱(chēng)是Genome Anlysis Toolkit究孕，顧名思義属提，是一套用于分析基因組的工具箱。主要功能是尋找變異位點(diǎn)和基因分型瓶佳，但是實(shí)際上功能超多返十，導(dǎo)致初學(xué)者都不知道從何學(xué)習(xí)GATK荣回。
為了幫助聽(tīng)說(shuō)過(guò)GATK的人更好的入門(mén)GATK察迟， 我，作為一個(gè)初學(xué)者掖疮，在這里介紹一下我的經(jīng)驗(yàn)初茶。
雖然GATK的功能超級(jí)多，但是主要可以歸為以下幾個(gè)方面：
- 診斷和質(zhì)量控制工具(Diagnostics and Quality Control Tools)
- 序列數(shù)據(jù)處理工具(Sequence Data Processing Tools)
- 變異位點(diǎn)探索工具(Variant Discovery Tools)
- 變異位點(diǎn)評(píng)估工具(Variant Evaluation Tools)
- 變異位點(diǎn)操作工具(Variant Manipulation Tools)
- 注釋模塊
- 讀段（reads）過(guò)濾
- 資源文件解碼工具
- 參考序列實(shí)用工具
如何快速建立GATK的心理表征
然后這里面的每一項(xiàng)點(diǎn)開(kāi)都有好多內(nèi)容浊闪，我第一次點(diǎn)開(kāi)的時(shí)候恼布，也是一臉茫然，不知道從何入手搁宾。
但是根據(jù)《認(rèn)知學(xué)習(xí)法》折汞，最好的學(xué)習(xí)方式就是“不要慫，直接上”盖腿，找到一個(gè)已有流程爽待，先把代碼敲上去，然后慢慢理解每一行代碼的作用翩腐，建立一個(gè)模糊的心理表征鸟款，然后循序漸進(jìn)，慢慢學(xué)習(xí)其他工具茂卦，最后就能熟練使用GATK了何什。
所以記下來(lái)主要的任務(wù)，就是帶大家建立關(guān)于GATK的模糊概念等龙。
mapping-by-sequencing其中一個(gè)重要環(huán)節(jié)就是“SNP calling”富俄，我最初用的是samtools和bcftools,結(jié)果的variant特別多（估計(jì)很多是假陽(yáng)性）.雖然最后還是找到了causual mutation, 但是為了保證今后causual mutation的準(zhǔn)確性，我發(fā)現(xiàn)了有文章使用了GATK而咆。他給的代碼如下：
1. Add read groups (Picard tools)
AddOrReplaceReadGroups.jar I=sorted.bam_file O=s1.rg.bam RGLB=genome RGPL=ILLUMINA
RGPU=GATKv4 RGSM=sample_name VALIDATION_STRINGENCY=LENIENT
2. Mark duplicates (Picard tools)
MarkDuplicates.jar INPUT=s1.rg.bam OUTPUT=s2.dedup.bam ASSUME_SORTED=TRUE
VALIDATION_STRINGENCY=LENIENT METRICS_FILE=s2.dedup.metrics
3. Index (samtools)
samtools index s2.dedup.bam
4. Realign reads (create intervals first, then do IndelRealigner) (GATK)
GenomeAnalysisTK.jar -I s2.dedup.bam -R ref_file -T RealignerTargetCreator -o s3.intervals
GenomeAnalysisTK.jar -T IndelRealigner -I s2.dedup.bam -R ref_file -targetIntervals s3.intervals -o
s4.realn.bam
5. Unified genotyper (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s4.realn.bam -glm BOTH -o s5.UG1.vcf -mbq 30 -nt
4
6. Base score recalibrator (GATK)
GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s5.UG1.vcf -o s6.recal
7. Print Reads (GATK)
GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s6.recal -o s7.recal.bam
8. Unified Genotype (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s7.recal.bam -glm BOTH -o s8.UG2.vcf -mbq 30
9. Base score recalibrator (GATK)
GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s8.UG1.vcf -o s9.recal
10. Print Reads (GATK)
GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s9.recal -o s10.recal.bam
11. Unified Genotyper (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s10.recal.bam -glm BOTH
第一步： 不管三七二十一直接把代碼自己敲一遍。
第二步： 讀每行代碼的解釋?zhuān)斫馑囊鈭D
- 用了
picard tools的AddOrReplaceReadGroups.jar加了read group幕袱，根據(jù)以前百度GATK的經(jīng)驗(yàn)暴备，了解GATK要求bam文件的header必須包含@RG，所以這一步應(yīng)該是前面比對(duì)時(shí)候们豌，沒(méi)有在參數(shù)中增加相應(yīng)部分涯捻。所以如果我在比對(duì)的時(shí)候增加了這個(gè)參數(shù)，這一步就可以免了望迎。 - 用
picard tools的MarkDuplicates.jar去重障癌。這是因?yàn)槎鷾y(cè)序有一部橋式PCR擴(kuò)增底物的過(guò)程，在100x以上出現(xiàn)reads重復(fù)辩尊，很大可能是PCR擴(kuò)增重復(fù)了涛浙，所以可以直接去掉。 -
samtools的index，建立索引轿亮，提高檢索效率疮薇。 - 先建立indel區(qū)間文件，然后對(duì)該區(qū)域進(jìn)行reads重比對(duì)我注。這一步我沒(méi)有經(jīng)驗(yàn)按咒，然后我在GATK的論壇搜索了這個(gè)工具，發(fā)現(xiàn)原因是indel會(huì)導(dǎo)致附近的錯(cuò)配但骨，所以需要借由這一步降低indel附近的假陽(yáng)性。但是最新的
HaplotypeCaller or MuTect2已經(jīng)不需要了奔缠，這些工具的haplotype組裝步奏效果類(lèi)似掠抬。這里的UnifiedGenotyper or the original MuTect.還是要的。 - 使用
UnifiedGenotyper尋找variant添坊。論壇搜索發(fā)現(xiàn)剿另，現(xiàn)在推薦用HaplotypeCaller，效果更好贬蛙。 - 根據(jù)上一步得到vcf文件對(duì)第四步得到BAM文件進(jìn)行堿基重校準(zhǔn)，得到校準(zhǔn)所需的文件阳准。
- 使用
PrintReads根據(jù)第6步得到的文件氛堕，對(duì)第四步得到的BAM文件進(jìn)行校準(zhǔn) - 根據(jù)重校準(zhǔn)的BAM文件重新找variant野蝇。
- 根據(jù)第8步的文件進(jìn)行校準(zhǔn)
- 輸出第二次校準(zhǔn)后BAM文件
- 第三次尋找variant。
第三步：根據(jù)上述過(guò)程對(duì)擬南芥中使用GATK尋找variant建立大致的認(rèn)識(shí)：
- 常規(guī)步驟： 先比對(duì)赠摇，比對(duì)后把SAM轉(zhuǎn)成BAM然后排序，建立索引浅蚪，去除PCR重復(fù)
- 對(duì)indel附近進(jìn)行重新比對(duì)
- 由于擬南芥沒(méi)有已知的VCF文件讓UnifiedGenotyper學(xué)習(xí)藕帜，所以需要通兩輪堿基校準(zhǔn)和variant calling降低假陽(yáng)性
第四步： 再看看別的流程，不斷修改最初的認(rèn)識(shí)惜傲。我通過(guò)翻文獻(xiàn)又找到了一個(gè)perl腳本
# 第一行洽故，UnifiedGenotyper
# 輸出read-$i.snv-gatk-snp.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T UnifiedGenotyper -R $cfg{genome} -I read-$i.recal.bam -D known_snp.conv -A AlleleBalance -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200 -glm SNP -out_mode EMIT_VARIANTS_ONLY -log log.UnifiedGenotyper-snp.$i -o read-$i.snv-gatk-snp.vcf");
# 第二行，VariantFiltration,
# 過(guò)濾read-$i.snv-gatk-snp.vcf得read-$i.snv-gatk-snp.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T VariantFiltration -R $cfg{genome} -V read-$i.snv-gatk-snp.vcf --clusterWindowSize 10 --filterExpression \"QUAL < 30.0 || QD < 5.0\" --filterName \"HARD_TO_VALIDATE\" -log log.VariantFiltration-snp.$i -o read-$i.snv-gatk-snp.fil.vcf");
# 第三行 UnifiedGenotyper
# 輸出read-$i.snv-gatk-indel.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T UnifiedGenotyper -R $cfg{genome} -I read-$i.recal.bam -D known_snp.conv -A AlleleBalance -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200 -glm INDEL -out_mode EMIT_VARIANTS_ONLY -log log.UnifiedGenotyper-indel.$i -o read-$i.snv-gatk-indel.vcf");
# 第四行 VariantFiltration
# 過(guò)濾read-$i.snv-gatk-indel.vcf得read-$i.snv-gatk-indel.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T VariantFiltration -R $cfg{genome} -V read-$i.snv-gatk-indel.vcf --clusterWindowSize 10 --filterExpression \"QUAL < 10.0\" --filterName \"HARD_TO_VALIDATE_1\" --filterExpression \"MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)\" --filterName \"HARD_TO_VALIDATE_2\" -log log.VariantFiltration-indel.$i -o read-$i.snv-gatk-indel.fil.vcf");
#第五行 CombineVariants
# 結(jié)合read-$i.snv-gatk-snp.fil.vcf和read-$i.snv-gatk-indel.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T CombineVariants -R $cfg{genome} --variant read-$i.snv-gatk-snp.fil.vcf --variant read-$i.snv-gatk-indel.fil.vcf -o read-$i.snv-gatk.vcf.tmp -log log.CombineVariants.$i");
好像和之前建立的模型不太一樣呀盗誊，這里變成先找variant然后過(guò)濾了时甚。
但是隘弊，這些流程所用的工具基本都來(lái)自于
- 序列數(shù)據(jù)處理工具(Sequence Data Processing Tools)
- 變異位點(diǎn)探索工具(Variant Discovery Tools)
- 變異位點(diǎn)評(píng)估工具(Variant Evaluation Tools)
其實(shí)都是先找到variant，然后通過(guò)各種方法降低假陽(yáng)性撞秋。
小結(jié)
總結(jié)一下：
GATK常用套路就是先找variant长捧，然后降低假陽(yáng)性(例如BQSR, VQSR或直接過(guò)濾等)
有了這樣一個(gè)大致印象后，我再去學(xué)習(xí)GATK Best Practices吻贿，就稍微有點(diǎn)思路了串结。
這些就是我這段日子學(xué)習(xí)GATK的感悟，不涉及到具體的操作舅列。希望大家能夠提供更多GATK相關(guān)的pipeline肌割，讓我能夠不斷更新自己的心理表征。
吐槽
本來(lái)我是想麻煩我同學(xué)讓他向師兄師姐要下他們的GATK流程帐要，便于我借鑒模仿學(xué)習(xí)把敞。但是他說(shuō)他們老板認(rèn)為這是機(jī)密，不能隨便透露榨惠。那么問(wèn)題來(lái)了奋早，假設(shè)你不把代碼公開(kāi)，那么別人無(wú)法重復(fù)出你的生物信息分析結(jié)果赠橙，那該怎么辦耽装？當(dāng)然吐槽歸吐槽，我也能理解期揪，畢竟要是被商業(yè)公司拿走掉奄，就可以拿出去掙錢(qián)了。
不過(guò)我相信凤薛，別人只能模仿你的代碼姓建，卻無(wú)法模仿你的思維模式，所以我很樂(lè)意把自己的代碼分享出去缤苫，只有經(jīng)得起檢驗(yàn)的代碼才是好代碼速兔。
