最近在做blood tumor paired外显子分析流程，看了一下GATK的论坛，在关于用对照样本作过滤的问题上，曾经的设置质量参数作硬性过滤条件的方法已经out，现在GATK best practice流程里是用大量正常对照样本训练一个正常模型，应该是用机器学习和算法吧(反正我不懂~).
那么基于我的需求也应该走Mutect2的CreateSomaticPanelOfNormals,GATK版本是4.1.4.0，然而搜了一圈没发现适用于最新版的中文教程，试了一下官方的流程好多坑，郁闷了两天还好最后跑通了，分享一下流程,非常感谢生信技能树健明大神的教导和鼓励，在这一步卡住好几天马上要放弃的时候，并且再一次帮助提升这篇小笔记的视觉效果，哈哈~~
官方流程链接，注意是最新版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这个参数对于双端测序也要加。
第二步
必须用GenomicsDBImport 创建文件夹genomicsdb-workspace-path，而不像之前版本生成一个正常vcf合集就结束，这个文件夹里有json文件需要下一步调用，在没跑通的时候我也想用以前版本combinedgvcf替代,事实证明不可以，就要乖乖按流程来并且流程里还有坑…..
先把上一步生成的vcf文件准备好，多参数输入vcf :
samples=$(find .|sed 's/.\///' | grep -E 'vcf.gz$' | sed 's/^/-V /')最坑的是-L 参数，一开始用了靶向bed文件，跑了一晚上三分之一都不到还生成了好几个T的文件，吓的赶紧停掉删除，论坛上搜最后发现这个GenomicsDBImport工具很矫情，貌似不能接受靶向捕获这种很多非连续区域bed文件，最后找到了方法，参数要分别接受每条染色体最后合成一个变量interval
制作interval内容如下：
-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，生成存放结果文件夹 ponall
gatk GenomicsDBImport --java-options "-Xmx32G -Djava.io.tmpdir=./" \ -R ${ref} ${interval} \ --genomicsdb-workspace-path ponall ${samples} 顺利运行了十几分钟就完成，结果文件也不大……..
诡异的是日志文件里只显示了某一条chr的一个位点最后显示success
用SelectVariants工具检查刚才生成的结果里面是否全面覆盖所有染色体
gatk SelectVariants -R ${ref} -V gendb://ponall -O outputcheck.vcf看过之后放心了，全部染色体都在，位点也很多。
第三步
终于回到pon了，用到上一步生成的ponall文件夹内容
gatk CreateSomaticPanelOfNormals --java-options "-Xmx32G -Djava.io.tmpdir=./" \ -V gendb://ponall \ -R ${ref} -O ponall.vcf.gz终于顺利完成~~~
网友评论