python - gatk VariantRecalibrator 上的 Snakemake
问题描述
我是使用snakemake 的新手,在snakemake 上执行步骤gatk VariantRecalibrator 时遇到问题,它产生错误,但脚本在非snakemake 格式时可以正常运行。
import snakemake.io
import os
REF="/data/data/reference/refs/ucsc.hg19.fasta"
HM="/data/data/variant_call/hapmap_3.3.hg19.sites.vcf"
OMNI="/data/data/variant_call/1000G_omni2.5.hg19.sites.vcf"
SNPS="/data/data/variant_call/1000G_phase1.snps.high_confidence.hg19.sites.vcf"
DBSNP="/data/data/variant_call/dbsnp_138.hg19.vcf"
NAME="CHS"
rule all:
input: "VCFs/{name}.recal.vcf".format(name=NAME),
"VCFs/{name}.output.tranches".format(name=NAME)
rule vqsr:
input: vcf="VCFs/SRS008640.raw.vcf",
ref=REF,
hm=HM,
omni=OMNI,
snps=SNPS,
dbsnp=DBSNP
output: recal="VCFs/{name}.recal.vcf".format(name=NAME),
tranches="VCFs/{name}.output.tranches".format(name=NAME),
rscript="VCFs/{name}.output.plots.R".format(name=NAME)
params: java_opts="-Xmx16g"
shell: "gatk --java-options -Xmx16g VariantRecalibrator \
-R {input.ref} \
-V {input.vcf} \
--resource:hapmap,known=false,training=true,truth=true,prior=15.0 {input.hm} \
--resource:omni,known=false,training=true,truth=false,prior=12.0 {input.omni} \
--resource:1000G,known=false,training=true,truth=false,prior=10.0 {input.snps} \
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 {input.dbsnp} \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode SNP \
-O {output.recal} \
--tranches-file {output.tranches} \
--rscript-file {output.rscript}"
错误:设置系统属性 GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') 以打印堆栈跟踪。[Thu May 30 08:05:30 2019] 规则 vqsr 中的错误:jobid:1 输出:VCFs/CHS.recal.vcf、VCFs/CHS.output.tranches、VCFs/CHS.output.plots.R
RuleException:
CalledProcessError in line 28 of /data/data/Samples/snakemake-example/WGS-test/step7.smk:
Command ' set -euo pipefail; gatk --java-options -Xmx16g VariantRecalibrator -R /data/data/reference/refs/ucsc.hg19.fas ta -V VCFs/SRS008640.raw.vcf --resource:hapmap,known=false,training=true,truth=true,prior=15.0 /data/data/variant_call /hapmap_3.3.hg19.sites.vcf --resource:omni,known=false,training=true,truth=false,prior=12.0 /data/data/variant_call/1000 G_omni2.5.hg19.sites.vcf --resource:1000G,known=false,training=true,truth=false,prior=10.0 /data/data/variant_call/1000G _phase1.snps.high_confidence.hg19.sites.vcf --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /data/data/ variant_call/dbsnp_138.hg19.vcf -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode SNP -O VCFs/CHS. recal.vcf --tranches-file VCFs/CHS.output.tranches --rscript-file VCFs/CHS.output.plots.R ' returned non-zero exit sta tus 2.
File "/data/data/Samples/snakemake-example/WGS-test/step7.smk", line 28, in __rule_vqsr
File "/root/miniconda3/envs/bioinfo/lib/python3.6/concurrent/futures/thread.py", line 56, in run
Removing output files of failed job vqsr since they might be corrupted:
VCFs/CHS.recal.vcf, VCFs/CHS.output.tranches, VCFs/CHS.output.plots.R
Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
Complete log: /data/data/Samples/snakemake-example/WGS-test/.snakemake/log/2019-05-30T065011.676785.snakemake.log
如果我使用相同的代码,我可以运行以创建 recal 文件和 tranches,并且可以转到下一步 applyvqsr,但是如果我将它放入 snakemake 它有错误并且第 27 行是 gatk --java-options -Xmx16g VariantRecalibrator是错误,但我不知道它是什么错误。请指教。
解决方案
我建议您使用--printshellcmds
参数运行 Snakemake。这将为您提供它在您的 shell 中运行的确切命令(而不是您在来自 python 胆量的异常中获得的命令)。您可以复制该命令并手动运行。
让我澄清一下:Snakemake 不是按照您在本shell
节中指定的那样运行命令,而是运行一个子进程,添加其他参数并设置环境变量(两者都因 shell 和操作系统而异),然后它使用 python 工具运行一个应用程序。使用参数运行 Snakemake--printshellcmds
将允许您查看没有格式问题的命令(当所有通配符都被解析时),但没有像set -euo pipefail;
.
推荐阅读
- google-sheets - 如何使用谷歌表格条件格式选项将列与单元格 id 而不是值进行比较?
- angular - 如何在 Angular 6 中使用 ng-for 处理 JSON 非结构化数据
- react-redux - 我了解 Redux 状态可用于容器,但动作创建者如何通过道具对容器可用,
- android-studio - 属性 application@appComponentFactory value=(android.support.v4.app.CoreComponentFactory)
- ios - 可以在 iOS 中配置/编程 VPN Killswitch 吗?
- amazon-web-services - 该区域和这些 S3 配置可能有什么问题?
- c# - 在 asp.net mvc 中预期条件的上下文中指定的非布尔类型的表达式
- java - 从 Spring AOP 方面返回自定义 HttpServletResponse 不起作用
- php - 如何在php中包含函数后检索变量
- kubernetes - Kubernetes 中两个或多个节点的 Nginx 请求