Snakemake: ошибка при попытке создать несколько выходных файлов

Я пишу конвейер snakemake, чтобы взять общедоступные файлы sra, преобразовать их в файлы fastq, а затем запустить их через выравнивание, пиковый вызов и регрессию оценки LD.

У меня проблема в правиле под названием SRA2fastq ниже, в котором я использую parallel-fastq-dump для преобразования файлов SRA в файлы с парным концом fastq. Это правило генерирует два вывода для каждого файла SRA: SRRXXXXXXX_1 и SRRXXXXXXX_2.

Вот мой файл конфигурации:

samples:
    fullard2018_NpfcATAC_1: SRR5367824
    fullard2018_NpfcATAC_2: SRR5367798
    fullard2018_NpfcATAC_3: SRR5367778
    fullard2018_NpfcATAC_4: SRR5367754
    fullard2018_NpfcATAC_5: SRR5367729

И вот несколько первых правил моего Snakefile:

# read config info into this namespace
configfile: "config.yaml"
print (config['samples'])

rule all:
    input:
        expand("fastq_files/{SRA}_{num}.fastq.gz", SRA=[config['samples'][x] for x in config['samples']], num=[1,2]),
        expand("FastQC/{SRA}_{num}_fastqc.html", SRA=[config['samples'][x] for x in config['samples']], num=[1,2]),
        "FastQC/fastq_multiqc.html",
        expand("peak_files/{sample}_peaks.blrm.narrowPeak", sample=config['samples']),
        "peak_files/Fullard2018_peaks.mrgd.blrm.narrowPeak",
        expand("LD_annotation_files/Fullard_2018.{chr}.l2.ldscore.gz", chr=range(1,23))

rule SRA_prefetch:
    params:
        SRA="{SRA}"
    output:
        "/home/c1477909/ncbi/public/sra/{SRA}.sra"
    log:
        "logs/prefetch/{SRA}.log"
    shell:
        "prefetch {params.SRA}"

rule SRA2fastq:
    input:
        "/home/c1477909/ncbi/public/sra/{SRA}.sra"
    output:
        "fastq_files/{SRA}_1.fastq.gz",
        "fastq_files/{SRA}_2.fastq.gz"
    log:
        "logs/SRA2fastq/{SRA}.log"
    shell:
        """
        parallel-fastq-dump --sra-id {input} --threads 8 \
        --outdir fastq_files --split-files --gzip
        """

rule fastqc:
    input:
        rules.SRA2fastq.output
    output:
        # Output needs to end in '_fastqc.html' for multiqc to work
        html="FastQC/{SRA}_{num}_fastqc.html"
    log:
        "logs/FASTQC/{SRA}_{num}.log"
    wrapper:
        "0.27.1/bio/fastqc"

rule multiqc_fastq:
    input:
        lambda wildcards: expand("FastQC/{SRA}_{num}_fastqc.html", SRA=[config['samples'][x] for x in config['samples']], num=[1,2])
    output:
        "FastQC/fastq_multiqc.html"
    wrapper:
        "0.27.1/bio/multiqc"

rule bowtie2:
    input:
        sample=lambda wildcards: expand("fastq_files/{SRA}_{num}.fastq.gz", SRA=config['samples'][wildcards.sample], num=[1,2])
    output:
        "bam_files/{sample}.bam"
    log:
        "logs/bowtie2/{sample}.txt"
    params:
        index=config["index"],  # prefix of reference genome index (built with bowtie2-build),
        extra=""
    threads: 8
    wrapper:
       "0.27.1/bio/bowtie2/align"

Однако, когда я запускаю Snakefile, я получаю следующую ошибку:

Error in job SRA2fastq while creating output files fastq_files/SRR5367754_1.fastq.gz, fastq_files/SRR5367754_2.fastq.gz

Я видел эту ошибку много раз раньше, и обычно она возникает, когда имя выходного файла, сгенерированного программой, не совсем соответствует имени выходного файла, которое вы указываете в соответствующем правиле snakemake. Однако здесь дело обстоит не так, как если бы я запускаю команду snakemake, генерируемую для этого конкретного правила отдельно, файлы создаются, как ожидалось, и имена файлов совпадают. Вот пример одного экземпляра правила, взятого после запуска snakemake -np:

rule SRA2fastq:
    input: /home/c1477909/ncbi/public/sra/SRR5367779.sra
    output: fastq_files/SRR5367779_1.fastq.gz, fastq_files/SRR5367779_2.fastq.gz
    log: logs/SRA2fastq/SRR5367779.log
    jobid: 18
    wildcards: SRA=SRR5367779

    parallel-fastq-dump --sra-id /home/c1477909/ncbi/public/sra/SRR5367779.sra --threads 8 --outdir fastq_files --split-files --gzip

Обратите внимание, что выходные файлы, сгенерированные командой parallel-fastq-dump, запущенной отдельно (т.е. без использования snakemake), имеют имена, указанные в правиле SRA2fastq:

ls fastq_files
SRR5367729_1.fastq.gz  SRR5367729_2.fastq.gz

Я немного озадачен этим, поскольку эту ошибку обычно легко исправить, но я не могу понять, в чем проблема. Я попытался изменить раздел вывода SRA2fastq на:

    output:
        file1="fastq_files/{SRA}_1.fastq.gz",
        file2="fastq_files/{SRA}_2.fastq.gz"

Однако это вызывает ту же ошибку. Я также пробовал просто указать один выходной файл, но позже это повлияет на правило bowtie2, так как я получаю input files missing ошибку.

Есть идеи, что здесь происходит? Что-то мне не хватает при поиске нескольких выходных файлов в одном правиле?

Огромное спасибо


person Darren    schedule 21.08.2018    source источник
comment
Похоже, это должно работать. Попробовать средство отслеживания проблем на bitbucket?   -  person The Unfun Cat    schedule 21.08.2018
comment
Каков статус команды оболочки? Попробуйте запустить свою команду parallel-fastq-dump, а затем echo $?. Возможно, программное обеспечение возвращает код ошибки, даже если оно работает.   -  person The Unfun Cat    schedule 21.08.2018
comment
@TheUnfunCat, спасибо, но точно, что вы имеете в виду во втором заявлении там ... не уверен, как `` проверить статус команды оболочки '' и нужно ли мне это делать при запуске snakemake или при запуске команды в качестве проверки напрямую в командной строке.   -  person Darren    schedule 21.08.2018
comment
Когда вы запускаете команды оболочки, они возвращаются с кодом выхода. Snakemake исследует это, чтобы узнать, успешно ли работала программа. Вы должны сделать это по команде в качестве проверки. Даже если это не причина вашей ошибки, полезно знать: shellscript.sh/exitcodes.html   -  person The Unfun Cat    schedule 22.08.2018