Snakemake на gatk
Я новичок в использовании змейки, у меня возникла проблема при выполнении шага gatk VariantRealibrator на змейке, он выдал ошибку, но скрипт может работать без ошибок, если он не в формате змейки.
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'), чтобы распечатать трассировку стека. [Чт, 30 мая, 08:05:30 2019] Ошибка в правиле 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 и транши, и могу перейти к следующему шагу applyvqsr, однако, если я поставлю его в змейку, у него будет ошибка, а строка 27 будет gatk --java-options -Xmx16g VariantRekalator. это ошибка, но я не знаю, что это за ошибка. Пожалуйста посоветуй.
Ответы 1
Я бы посоветовал запускать Snakemake с параметром --printshellcmds. Это даст вам точную команду, которую он запускает в вашей оболочке (а не ту, которую вы получаете в исключении, которое исходит из кишок python). Вы можете скопировать эту команду и запустить вручную.
Поясню: Snakemake запускает команду не так, как вы указали в разделе shell, а запускает подпроцесс, добавляет дополнительные параметры и устанавливает переменные среды (обе различаются в зависимости от оболочки и ОС), затем использует средства Python для запустить приложение. Запуск Snakemake с параметром --printshellcmds позволит вам увидеть команду без проблем с форматированием (когда разрешены все подстановочные знаки), но без лишних деталей, таких как set -euo pipefail;.
он может создавать .recal.vcf, output.plots.R и output.tranches, когда я запускаю командную строку вручную
Вы запускали Snakemake с --printshellcmds? что было на выходе?
Создание DAG заданий... Количество заданий: количество заданий 1 все 1 vqsr 2 [Пн, 3 июня, 03:39:18 2019] правило vqsr: ввод: VCFs/SRS008640.raw.vcf, /data/data/reference/refs/ ucsc.hg19.fasta, /data/data/variant_call/hapm ap_3.3.hg19.sites.vcf, /data/data/variant_call/1000G_omni2.5.hg19.sites.vcf, /data/data/variant_call/1000G _phase1 .snps.high_confidence.hg19.sites.vcf, /data/data/variant_call/dbsnp_138.hg19.vcf вывод: VCFs/CHS.recal.vcf, VCFs/CHS.output.tranches, VCFs/CHS.output.plots.R Идентификатор задания: 1
gatk --java-options -Xmx16g VariantRecalibrator -R /data/data/reference/refs/ucsc.hg19.fasta -V VCFs/S RS008640.raw.vcf --resource:hapmap,known=false, training=true,true =true, prior=15.0 /data/data/variant_call/hapmap_3.3.hg19.sites.vcf --resource:omni,known=false, training=true,true=false, prior=12.0 /data/ data /variant_call/1000G_omni2.5.hg19.sites.vcf --resource:1000G,known=false, training=true, правда=false, prior= 10.0 /data/data/variant_call/1000G_phase1.snps.high_confidence.hg 19.sites.vcf --
ресурс:dbsnp,known=true, training=false,true=false, prior=2.0 /data/data/variant_call/dbsnp_138.hg19.vcf -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -режим SNP -O VCFs/CHS.recal.vcf --tranches-file VCFs/CHS.out put.tranches --rscript-file VCFs/CHS.output.plots.R
[Пн, 3 июня, 03:39:18 2019] localrule all: input: VCFs/CHS.recal.vcf, VCFs/CHS.output.tranches, VCFs/CHS.output.plots.R jobid: 0 Количество заданий: количество заданий 1 all 1 vqsr 2 Это пробный запуск (флаг -n). Порядок заданий не отражает порядок выполнения.
Хорошо, теперь у нас есть точная команда, запущенная в оболочке. Что возвращает эта (точная) команда при запуске из оболочки?
Ошибка: установите системное свойство GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') для печати трассировки стека. [Чт, 30 мая, 08:05:30 2019] Ошибка в правиле vqsr: jobid: 1 вывод: VCFs/CHS.recal.vcf, VCFs/CHS.output.tranches, VCFs/CHS.output.plots.R
RuleException: CalledProcessError в строке 28 файла /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,true =true, prior=15.0 /data/data/variant_call /hapmap_3.3.hg19.sites.vcf --resource:omni,known=false, training=true, правда=false, prior=12.0 /data/data /variant_call/1000 G_omni2.5.hg19.sites.vcf --resource: 1000G, известный = ложь, обучение = правда, правда = ложь, ранее= 10,0
/data/data/variant_call/1000G _phase1.snps.high_confidence.hg19.sites.vcf --resource:dbsnp,known=true, training=false, true=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 VCF/CHS. recal.vcf --tranches-file VCFs/CHS.output.tranches --rscript-file VCFs/CHS.output.plots.R ' вернул ненулевой статус выхода 2. File "/data/data/Samples/snakemake- пример/WGS-test/step7.smk", строка 28, в __rule_vqsr
Файл «/root/miniconda3/envs/bioinfo/lib/python3.6/concurrent/futures/thread.py», строка 56, выполняется удаление выходных файлов неудачного задания vqsr, поскольку они могут быть повреждены: VCFs/CHS. recal.vcf, VCFs/CHS.output.tranches, VCFs/CHS.output.plots.R Выключение может занять некоторое время. Выход из-за сбоя выполнения задания. Посмотрите сообщение об ошибке выше. Полный журнал: /data/data/Samples/snakemake-example/WGS-test/.snakemake/log/2019-05-30T065011.676785.snakemake.log
Runtime.totalMemory()=7386169344 ************************************************* ***************************** Произошла ОШИБКА ПОЛЬЗОВАТЕЛЯ: невозможно выполнить команду Rscript: добавьте каталог Rscript в свою среду ${ПУТЬ} **************************************************** ************************** Установите системное свойство GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true'), чтобы распечатать трассировку стека .