snakemake: пропущено последнее правило, изменив все командные строки

ниже мои коды snakemake, если я не закомментирую код line28,29, который является правилом all-> input-> 1-я, 2-я командные строки, тогда я не могу получить последнее правило varscan_somatic , то есть результат пробного прогона выглядит так:

Job counts:
        count   jobs
        1       all
        25      mpileup_analysis
        25      normal_mpileup
        25      tumor_mpileup
        76
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

но если я закомментирую код line28,29, который является правилом all-> input-> 1-я, 2-я командные строки, то я могу получить последнее правило varscan_somatic, то есть сухой запустить вывод вроде этого:

Job counts:
        count   jobs
        1       all
        25      mpileup_analysis
        25      normal_mpileup
        25      tumor_mpileup
        25      varscan_somatic
        101
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

Я не знаю, почему это происходит? Кто-нибудь может дать мне совет? Большое спасибо за любую помощь.

import re
import os

mydict = dict()
with open("config.txt") as HD:
    for line in HD:
        line = line.rstrip()
        if line.startswith("#"):
            continue
        value,field = re.split("\s*=\s*",line)
        mydict[value] = field

VarScan    = mydict['VarScan']
SAMtools   = mydict['SAMtools']
REFERENCE  = mydict['REFERENCE']
PL_MPPANA  = mydict['MPPANA']

chrlist = [i for i in os.listdir("call_region") if i.endswith(".region.bed")]
def replace(str1):
    str2 = str1.replace(".region.bed","")
    return str2
chrlist = map(replace,chrlist)

configfile:"paired_test.yaml"

rule all:
    input:
#        expand("varscan_somatic/{sample}/{chrid}.normal.mpileup.analysis",sample=config['samples'],chrid=chrlist),
#        expand("varscan_somatic/{sample}/{chrid}.tumor.mpileup.analysis",sample=config['samples'],chrid=chrlist),
        expand("varscan_somatic/{sample}/{chrid}.snp.vcf",sample=config['samples'],chrid=chrlist),
        expand("varscan_somatic/{sample}/{chrid}.indel.vcf",sample=config['samples'],chrid=chrlist)

rule normal_mpileup:
    input:
        bam=lambda wc:config['samples'][wc.sample][3],
        bed="call_region/{chrid}.region.bed"
    output:
        "varscan_somatic/{sample}/{chrid}.normal.mpileup"
    log:
        "log/varscan_somatic/{sample}/{chrid}.normal.mpileup.log"
    shell:
        "{SAMtools} mpileup -f {REFERENCE} -l {input.bed} "
        "{input.bam} -d1000 -Q10 -q10 -o {output} "
        "1>{log} 2>&1"

rule tumor_mpileup:
    input:
        bam=lambda wc:config['samples'][wc.sample][1],
        bed="call_region/{chrid}.region.bed"
    output:
        "varscan_somatic/{sample}/{chrid}.tumor.mpileup"
    log:
        "log/varscan_somatic/{sample}/{chrid}.tumor.mpileup.log"
    shell:
        "{SAMtools} mpileup -f {REFERENCE} -l {input.bed} "
        "{input.bam} -d1000 -Q10 -q10 -o {output} "
        "1>{log} 2>&1"

rule mpileup_analysis:
    input:
        tumor="varscan_somatic/{sample}/{chrid}.tumor.mpileup",
        norml="varscan_somatic/{sample}/{chrid}.normal.mpileup"
    output:
        tumor="varscan_somatic/{sample}/{chrid}.tumor.mpileup.analysis",
        norml="varscan_somatic/{sample}/{chrid}.normal.mpileup.analysis"
    log:
        "log/varscan_somatic/{sample}/{chrid}.mpileup_analysis.log"
    shell:
        "{PL_MPPANA} {input.tumor} {output.tumor} 1>{log} 2>&1 "
        "&& {PL_MPPANA} {input.norml} {output.norml} 1>>{log} 2>&1"

rule varscan_somatic:
    input:
        tumor="varscan_somatic/{sample}/{chrid}.tumor.mpileup",
        norml="varscan_somatic/{sample}/{chrid}.normal.mpileup",
        temp1="varscan_somatic/{sample}/{chrid}.tumor.mpileup.analysis",
        temp2="varscan_somatic/{sample}/{chrid}.normal.mpileup.analysis"
    output:
        "varscan_somatic/{sample}/{chrid}.snp",
        "varscan_somatic/{sample}/{chrid}.indel"
    log:
        "log/varscan_somatic/{sample}/{chrid}.varscan_somatic.log"
    params:
        "varscan_somatic/{sample}/{chrid}",
        "--validation 1 --output-vcf 1"
    shell:
        "{VarScan} somatic {input.tumor} {input.norml} {params} 1>{log} 2>&1"

>>>config.txt
VarScan = /path/my/varscan
REFERENCE = /path/my/hg19.fa
SAMtools = /path/my/samtools
MPPANA = /path/my/mppana.pl

>>> paired.yaml
samples:
    S01:['S01','S01.bqsr.bam','S02','S02.bqsr.bam']

ниже: сводка snakemake без комментариев

snakemake -s bin/varscan_somatic_paired.py -np --forceall --summary
Building DAG of jobs...
output_file     date    rule    version log-file(s)     status  plan
varscan_somatic/S001/chr21.tumor.mpileup.analysis       Thu Sep 26 02:34:56 2019        mpileup_analysis        -       log/varscan_somatic/S001/chr21.mpileup_analysis.log ok      update pending
varscan_somatic/S001/chr21.normal.mpileup.analysis      Thu Sep 26 02:34:56 2019        mpileup_analysis        -       log/varscan_somatic/S001/chr21.mpileup_analysis.log ok      update pending
varscan_somatic/S001/chr10.tumor.mpileup.analysis       Wed Sep 25 22:56:59 2019        mpileup_analysis        -       log/varscan_somatic/S001/chr10.mpileup_analysis.log ok      update pending

ниже: сводка snakemake с комментариями

Building DAG of jobs...
output_file     date    rule    version log-file(s)     status  plan
varscan_somatic/S001/chr19.snp.vcf      Wed Sep 25 22:14:13 2019        varscan_somatic -       log/varscan_somatic/S001/chr19.varscan_somatic.log ok       update pending
varscan_somatic/S001/chr19.indel.vcf    Wed Sep 25 22:14:13 2019        varscan_somatic -       log/varscan_somatic/S001/chr19.varscan_somatic.log ok       update pending
varscan_somatic/S001/chr14.snp.vcf      Thu Sep 26 01:22:17 2019        varscan_somatic -       log/varscan_somatic/S001/chr14.varscan_somatic.log ok       update pending

введите здесь описание изображения

введите здесь описание изображения


snakemake
person Danielle    schedule 23.09.2019    source источник
comment
Я не вижу причин, по которым выполнение varscan_somatic должно зависеть от этих закомментированных строк. Попробуйте добавить параметр --summary в свою команду snakemake с комментариями и без них и опубликуйте результат, если он не слишком большой.   -  person dariober    schedule 23.09.2019
comment
Вы пробовали --forceall?   -  person Dmitry Kuzminov    schedule 23.09.2019
comment
@DmitryKuzminov Я пробовал --forceall, но ничего не менять .$ snakemake -s bin/varscan_somatic_paired.py -np --forceall Building DAG of jobs... Job counts: count jobs 1 all 25 mpileup_analysis 25 normal_mpileup 25 tumor_mpileup 76   -  person Danielle    schedule 26.09.2019
comment
@dariober Я попробовал --summary, и результат очень большой. Поэтому я выделю заголовок 3 строчки в вопросе. Я не знаю, могу ли я передать вам файлы в Stack Overflow?   -  person Danielle    schedule 26.09.2019
comment
@Danielle Здесь вы можете представить картину зависимостей. Попробуйте snakemake --dag.   -  person Dmitry Kuzminov    schedule 26.09.2019
comment
@DmitryKuzminov Я загрузил картинки, но вижу, что они не загружаются. Думаю, может быть, из-за брандмауэра я не могу загружать картинки в stackoverflow.   -  person Danielle    schedule 29.09.2019
comment
Я вижу, что изображение DAG соответствует результатам пробного прогона или сумме. Я думаю, что это ошибка змейки.   -  person Danielle    schedule 29.09.2019
comment
@Danielle Я не уверен, что предоставленный вами код является действительным кодом, который вы запускаете. По крайней мере, идентификация комментариев недействительна. Попробуйте упростить пример. Сначала оставьте единственную пару образец / chrid в config. Затем замените shellsection простым touch {output}. Приведите минимальный пример, который не работает. Я попробую это на своей машине и подтвердю, если это не сработает. PS Я не верю, что это ошибка в Snakemake.   -  person Dmitry Kuzminov    schedule 29.09.2019

Ответы (1)


arrow_upward
0
arrow_downward

Этот код неверен: до исправления:

def replace(str1):
    str2 = str1.replace(".region.bed","")
    return str2
chrlist = map(replace,chrlist)

после исправления:

def replace(str1):
    str2 = str1.replace(".region.bed","")
    return str2
chrlist = list(map(replace,chrlist))

тогда все ок.

person Danielle    schedule 29.09.2019