找寻gvcf失败的原因

Posted on by

最近走我整理和搭建好的：最新版针对RNA-seq数据的GATK找变异流程， 如果样本样品是正常运行，会输出：

920M Nov 9 02:07 SRR2016956_gatk.gvcf
 12M Nov 9 02:07 SRR2016956_gatk.gvcf.idx
5.5M Nov 9 00:25 SRR2016956_recal.bai
9.0G Nov 9 00:25 SRR2016956_recal.bam
840K Nov 9 00:09 SRR2016956_recal.table
6.2M Nov 8 23:48 SRR2016956_rmd_split.bai

因为我只想保留recal后的bam和call出来的gvcf文件，但是发现有些样本根本就走不通这个流程，就需要debug了。

6.6G Nov 15 08:33 SRR2016955_rmd.bam
6.5M Nov 15 08:33 SRR2016955_rmd.bam.bai
 0 Nov 15 08:33 SRR2016955_rmd_split.bai
 0 Nov 15 08:33 SRR2016955_rmd_split.bam
 134 Nov 15 10:32 SRR2016955_rmd_split_add.bam

也就是说错在：

 time $gatk SplitNCigarReads -R $GENOME \
 -I ${sample}_rmd.bam \
 -O ${sample}_rmd_split.bam

这个步骤非常耗费时间，需要至少100min，所以需要一些技巧来调试。首先需要查看日志，因为GATK的日志非常多，所以需要对比查看，然后我发现运行错误的样本，多了下面的信息：

htsjdk.samtools.util.RuntimeIOException: Attempt to add record to closed writer.
 at htsjdk.samtools.util.AbstractAsyncWriter.write(AbstractAsyncWriter.java:57)
 at htsjdk.samtools.AsyncSAMFileWriter.addAlignment(AsyncSAMFileWriter.java:53)
 at org.broadinstitute.hellbender.utils.read.SAMFileGATKReadWriter.addRead(SAMFileGATKReadWriter.java:21)
 at org.broadinstitute.hellbender.tools.walkers.rnaseq.OverhangFixingManager.writeReads(OverhangFixingManager.java:358)
 at org.broadinstitute.hellbender.tools.walkers.rnaseq.OverhangFixingManager.flush(OverhangFixingManager.java:338)
 at org.broadinstitute.hellbender.tools.walkers.rnaseq.SplitNCigarReads.closeTool(SplitNCigarReads.java:196)
 at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1052)
 at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
 at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
 at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
 at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
 at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
 at org.broadinstitute.hellbender.Main.main(Main.java:292)

很简单就谷歌搜索定位到了讨论帖：https://gatkforums.broadinstitute.org/gatk/discussion/11378/gatk4-splitncigarreads-runtimeioexception-attempt-to-add-record-to-closed-writer

非常感谢前辈们愿意花费十几至二十几分钟写自己的疑惑然后跟开发者公开讨论让我们后学者可以快速解决bug

原来是临时目录问题

I figured out the issue (at least for me). It stems from where SplitNCigarReads is writing the temporary files. For me, it's writing them to the cluster which has very limited disk space. When I redirected this using --TMP_DIR /my/scratch/space everything went smoothly.

The part that still confuses me is that I had already set export _JAVA_OPTIONS=-Djava.io.tmpdir=/my/scratch/space. This is not getting picked up by SplitNCigarReads in GATK4 as I would have expected. After much experimenting I started with a clean environment and simply set --TMP_DIR /my/scratch/space only which worked.

This seems a bit "buggy" to me and it would be great if the GATK development team could look into it and pass Djava.io.tmpdir to --TMP_DIR if possible.

这个时候，就必须回过头看看另外一个：虽然不知道为什么但是我可以解决这个bug，但它就是被解决了 因为那个时候不知道，所以就糊里糊涂的解决了问题而已，但是现在遇到了类似的问题，仍然是很头疼！

那个时候同样的也是 —TMP_DIR 捣鬼，所以我才会在运行GATK命令的时候设置这个 java.io.tmpdir ，通过下面的方式：

　java -Djava.io.tmpdir=/path/to/tmpdir

但是现在就是进退两难了，因为加上这个参数，命令提交到服务器就被kill，如果不加上这个参数，我们的临时目录又没有空间让我运行这个gatk命令。

首先需要简单了解一些java

java -version:查看JDK版本 
whereis java
which java （java执行路径）
echo $JAVA_HOME
echo $PATH

其中java可以操作的参数非常多，这次我们遇到的就是 java.io.tmpdir， 默认的临时文件路径！

java.version Java运行时环境版本
java.vendor Java运行时环境供应商
java.vendor.url Java供应商的 URL
java.home Java安装目录
java.vm.specification.version Java虚拟机规范版本
java.vm.specification.vendor Java虚拟机规范供应商
java.vm.specification.name Java虚拟机规范名称
java.vm.version Java虚拟机实现版本
java.vm.vendor Java虚拟机实现供应商
java.vm.name Java虚拟机实现名称
java.specification.version Java运行时环境规范版本
java.specification.vendor Java运行时环境规范供应商
java.specification.name Java运行时环境规范名称
java.class.version Java类格式版本号
java.class.path Java类路径
java.library.path 加载库时搜索的路径列表
java.io.tmpdir 默认的临时文件路径
java.compiler 要使用的 JIT 编译器的名称
java.ext.dirs 一个或多个扩展目录的路径
os.name 操作系统的名称
os.arch 操作系统的架构
os.version 操作系统的版本
file.separator 文件分隔符（在 UNIX 系统中是“/”）
path.separator 路径分隔符（在 UNIX 系统中是“:”）
line.separator 行分隔符（在 UNIX 系统中是“/n”）
user.name 用户的账户名称
user.home 用户的主目录
user.dir 用户的当前工作目录

其实临时文件并不是一直报错

因为HPC是一百多人一起使用，虽然说临时文件夹有几百个G，但是也架不住人多。

过两天再运行同样的程序居然就ok了，让我也很郁闷。

最后友情宣传生信技能树

​

Comments are closed.