Gatk VariantRecalibrator 上的 Snakemake

Question

我是使用 snakemake 的新手，我在对 snakemake 执行 step 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') 打印堆栈跟踪。 [2019 年 5 月 30 日星期四 08:05:30] 规则 vqsr 错误：职位编号：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 是错误的，但我不知道是什么错误。请指教。

Answer 1

我建议您使用 --printshellcmds 参数运行 Snakemake。这会给你确切的命令，它运行s 在你的 shell 中（不是你在来自 python 胆量的异常中得到的那个）。您可以手动复制该命令和运行。

让我澄清一下：Snakemake 运行不是您在 shell 部分中指定的命令，而是运行一个子进程，添加了额外的参数并设置了环境变量（两者因 shell 和 OS 而异），然后它使用 python 设施来运行应用程序。运行使用 --printshellcmds 参数的 Snakemake 可以让您看到没有格式问题的命令（当所有通配符都已解析时），但没有像 set -euo pipefail;.

这样的细节开销

Gatk VariantRecalibrator 上的 Snakemake

Snakemake on gatk VariantRecalibrator

python

bioinformatics

snakemake

vcf-variant-call-format