最近在做blood tumor paired外顯子分析流程，看了一下GATK的論壇，在關(guān)于用對照樣本作過濾的問題上览徒，曾經(jīng)的設(shè)置質(zhì)量參數(shù)作硬性過濾條件的方法已經(jīng)out，現(xiàn)在GATK best practice流程里是用大量正常對照樣本訓(xùn)練一個正常模型颂龙，應(yīng)該是用機器學(xué)習(xí)和算法吧(反正我不懂~).

那么基于我的需求也應(yīng)該走Mutect2的CreateSomaticPanelOfNormals,GATK版本是4.1.4.0习蓬，然而搜了一圈沒發(fā)現(xiàn)適用于最新版的中文教程，試了一下官方的流程好多坑措嵌，郁悶了兩天還好最后跑通了躲叼，分享一下流程,非常感謝生信技能樹健明大神的教導(dǎo)和鼓勵，在這一步卡住好幾天馬上要放棄的時候企巢，并且再一次幫助提升這篇小筆記的視覺效果押赊，哈哈~~

官方流程鏈接，注意是最新版GATK 4.1.4.0 MUTECT2的 CreateSomaticPanelOfNormals

<https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_mutect_CreateSomaticPanelOfNormals.php>

第一步

gatk Mutect2 --java-options "-Xmx32G -Djava.io.tmpdir=./" \    -R ${ref} -I ${id}*.applybqsr.bam  \    -max-mnp-distance 0 --independent-mates  \    -L ${targetv5} -O ${dirvcfm2}/${id}.raw.vcf.gz

這里掉了兩個坑包斑，一是-max-mnp-distance 0 一定要加流礁，不然下一步會報錯，再一個是—independent-mates這個參數(shù)對于雙端測序也要加罗丰。

第二步

必須用GenomicsDBImport 創(chuàng)建文件夾genomicsdb-workspace-path神帅，而不像之前版本生成一個正常vcf合集就結(jié)束，這個文件夾里有json文件需要下一步調(diào)用萌抵，在沒跑通的時候我也想用以前版本combinedgvcf替代,事實證明不可以找御，就要乖乖按流程來并且流程里還有坑…..
先把上一步生成的vcf文件準備好元镀，多參數(shù)輸入vcf :

samples=$(find .|sed 's/.\///' | grep -E 'vcf.gz$' | sed 's/^/-V /')

最坑的是-L 參數(shù)，一開始用了靶向bed文件霎桅，跑了一晚上三分之一都不到還生成了好幾個T的文件栖疑，嚇的趕緊停掉刪除，論壇上搜最后發(fā)現(xiàn)這個GenomicsDBImport工具很矯情滔驶，貌似不能接受靶向捕獲這種很多非連續(xù)區(qū)域bed文件遇革，最后找到了方法，參數(shù)要分別接受每條染色體最后合成一個變量interval
制作interval內(nèi)容如下：

-L chr1 -L chr2 -L chr3 -L chr4 -L chr5 -L chr6 -L chr7 -L chr8 -L chr9 -L chr10 -L chr11 -L chr12 -L chr13 -L chr14 -L chr15 -L chr16 -L chr17 -L chr18 -L chr19 -L chr20 -L chr21 -L chr22 -L chrX -L chrY

運行GenomicsDBImport揭糕，生成存放結(jié)果文件夾 ponall

gatk GenomicsDBImport --java-options "-Xmx32G -Djava.io.tmpdir=./" \     -R ${ref} ${interval}  \     --genomicsdb-workspace-path ponall   ${samples}

順利運行了十幾分鐘就完成萝快，結(jié)果文件也不大……..
詭異的是日志文件里只顯示了某一條chr的一個位點最后顯示success
用SelectVariants工具檢查剛才生成的結(jié)果里面是否全面覆蓋所有染色體

gatk SelectVariants -R ${ref} -V gendb://ponall  -O outputcheck.vcf

看過之后放心了，全部染色體都在著角，位點也很多揪漩。

第三步

終于回到pon了，用到上一步生成的ponall文件夾內(nèi)容

gatk CreateSomaticPanelOfNormals --java-options "-Xmx32G -Djava.io.tmpdir=./"  \     -V gendb://ponall  \     -R ${ref}  -O ponall.vcf.gz

終于順利完成~~~

?