如何使用 Snakemake 处理多个输出文件名约定

Question

我有一个场景，其中程序的输出文件名不一致，因此我很难定义 snakemake 规则。

我正在从 SRA 中检索数据，但我不知道该程序是 return 1 个文件还是 2 个文件

这里有两种情况：

情况一

Retrieve SE data - this generates 1 output file named SRR5597649_1.fastq.gz

parallel-fastq-dump --sra-id SRR5597649 --threads 16 --outdir out/ --gzip --skip-technical --readids --dumpbase --split-files --clip

情况 2

Retrieve PE data - this generates 2 output files SRR390728_1.fastq.gz  SRR390728_2.fastq.gz

parallel-fastq-dump --sra-id SRR390728 --threads 16 --outdir out/ --gzip --skip-technical --readids --dumpbase --split-files --clip

我想从 Snakemake 社区了解他们将如何解决上述情况以及如何定义个人规则以及 rule:all

注意：我在尝试@dariober 解决方案后向这个问题添加了更多信息

Snakefile

shell.executable("bash")
import pandas as pd
import glob
import os

configfile: "config.json"

data_dir=os.getcwd()

units_table = pd.read_table("Samples.txt")
samples= list(units_table.Samples.unique())

rule all:
    input:
       expand("raw_fastq/{sample}.dump.done", sample=samples),
       expand("results/salmon/{sample}_salmon_quant",sample=samples),
       expand("results/logs/salmon/{sample}.salmon.log",sample=samples)

rule clean:
     shell: "rm -rf .snakemake/"

include: 'rules/download_sra.smk'
include: 'rules/salmon_quant.smk'

download_sra.smk

rule download_sra:
    """
    Download RNA-Seq data from SRA.
    """
    output: touch("raw_fastq/{sample}.dump.done"),
    params:
        outdir = "raw_fastq",
        threads = 16
    priority:85
    shell: "parallel-fastq-dump --sra-id {wildcards.sample} ... "

salmon_quant.smk

rule salmon_quant:
     input:
             touchfile = 'raw_fastq/{sample}.dump.done',
             index = config['salmon_rat_gentrome_Index']
     output:
             directory("results/salmon/{sample}_salmon_quant")
     log:
             'results/logs/salmon/{sample}.salmon.log'
     priority: 85
     threads: 20
     run:
        fq= sorted(glob.glob(os.path.join('raw_fastq',  wildcards.sample + '_[12].fastq.gz')))

        if len(fq) == 1:
            shell("salmon quant -i {input.index} -l A  -r %s  ..." % fq) # Process as SE
        elif len(fq) == 2:
            shell("salmon quant -i {input.index} -l A -1 %s -2 %s ...." % (fq[0], fq[1])) # Process as PE
        else:
            raise Exception('Wrong number of fastq files!')

我面临的问题：当我尝试'dariober'建议的解决方案时，当%s被替换时似乎有问题shell 命令运行。例如，它被替换为

-r ['raw_fastq/SRR5597645_1.fastq.gz']

而它应该 已被替换为

-r raw_fastq/SRR5597645_1.fastq.gz

可能我在某处弄错了 - 如上面的代码所示，%s 导致沿方括号和单引号的标记。

此外，我想在 download_sra 规则中用 temp 指令包装它，以便在作业完成后自动删除 fastq 文件和 touchfiles。

在此先感谢大家。

Answer 1

一种解决方案可能是在 fastq-dump 完成后触摸一个文件，并将该文件用作其他规则的输入。

在您发布的情况下，实际的 fastq 文件的名称可以从 SRR id 重建，您可以根据后缀 _1.fastq.gz 或 [=12 判断您是否有 1 个或 2 个 fastq 文件=].

这是一个例子：

import glob
import os

sra_id= ['SRR1234', 'SRR4567']

rule all:
    input:
        expand('{sra_id}.dump.done', sra_id= sra_id),
        expand('{sra_id}.bam', sra_id= sra_id),

rule sra_dump:
    output:
        touch('{sra_id}.dump.done'),
    shell:
        r"""
        parallel-fastq-dump --sra-id {wildcards.sra_id} ... 
        """

rule align:
    input:
        '{sra_id}.dump.done',
    output:
        '{sra_id}.bam',
    run:
        fq= sorted(glob.glob(os.path.join(wildcards.sra_id, wildcards.sra_id + '_[12].fastq.gz')))

        if len(fq) == 1:
            shell("bwa mem %s ..." % fq) # Process as SE
        elif len(fq) == 2:
            shell("bwa mem %s %s ..." % (fq[0], fq[1])) # Process as PE
        else:
            raise Exception('Wrong number of fastq files!')

如何使用 Snakemake 处理多个输出文件名约定

How to handle mutliple output filename conventions with Snakemake

bioinformatics

snakemake