Snakemake

Question

我有大量的输入文件是这样组织的：

data/
├── set1/
│   ├── file1_R1.fq.gz
│   ├── file1_R2.fq.gz
│   ├── file2_R1.fq.gz
│   ├── file2_R2.fq.gz
|   :
│   └── fileX_R2.fq.gz
├── another_set/
│   ├── asdf1_R1.fq.gz
│   ├── asdf1_R2.fq.gz
│   ├── asdf2_R1.fq.gz
│   ├── asdf2_R2.fq.gz
|   :
│   └── asdfX_R2.fq.gz
:   
└── many_more_sets/
    ├── zxcv1_R1.fq.gz
    ├── zxcv1_R2.fq.gz
    :
    └── zxcvX_R2.fq.gz

如果您熟悉生物信息学——这些当然是来自配对末端测序运行的 fastq 文件。我正在尝试生成一个读取所有这些的 snakemake 工作流程，但我已经在第一条规则上失败了。这是我的尝试：

configfile: "config.yaml"

rule all:
    input:
        read1=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R1.fq.gz", output=config["output"]),
        read2=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R2.fq.gz", output=config["output"])

rule clip_and_trim_reads:
    input:
        read1=expand("{data}/{set}/{{sample}}_R1.fq.gz", data=config["raw_data"], set=config["sets"]),
        read2=expand("{data}/{set}/{{sample}}_R2.fq.gz", data=config["raw_data"], set=config["sets"])
    output:
        read1=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R1.fq.gz", output=config["output"]),
        read2=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R2.fq.gz", output=config["output"])
    threads: 16
    shell:
        """
        someTool -o {output.read1} -p {output.read2} \
        {input.read1} {input.read2}
        """

我无法将 clip_and_trim_reads 指定为目标，因为 Target rules may not contain wildcards. 我尝试添加 all 规则，但这给了我另一个错误：

$ snakemake -np
Building DAG of jobs...
WildcardError in line 3 of /work/project/Snakefile:
Wildcards in input files cannot be determined from output files:
'sample'

我也尝试使用 dynamic() 函数作为 all 规则，它确实找到了文件，但仍然给我这个错误：

$ snakemake -np
Dynamic output is deprecated in favor of checkpoints (see docs). It will be removed in Snakemake 6.0.
Building DAG of jobs...
MissingInputException in line 7 of /work/project/ladsie_002/analyses/scripts/2019-08-02_splice_leader_HiC/Snakefile:
Missing input files for rule clip_and_trim_reads:
data/raw_data/set1/__snakemake_dynamic___R1.fq.gz
data/raw_data/set1/__snakemake_dynamic___R2.fq.gz
data/raw_data/set1/__snakemake_dynamic___R2.fq.gz
data/raw_data/set1/__snakemake_dynamic___R1.fq.gz
[...]

我有一百多个不同的文件，所以我非常想避免创建包含所有文件名的列表。有什么想法可以实现吗？

Answer 1

我认为您误解了 snakemake 的工作原理。当你运行 snakemake 你定义你想要的输出时，要么在命令行上，要么将生成 Snakefile 中第一条规则的输入（你的规则全部）。由于你没有指定任何输出（snakemake -np），Snakemake 将尝试生成规则所有的输入。

你的规则的输入基本上都是：

"somepath/clipped_and_trimmed_reads/{sample}_R1.fq.gz"

不幸的是，Snakemake 不知道如何从中生成输出...我们需要告诉 Snakemake 使用哪些文件。我们可以手动执行此操作：

rule all:
    input:
        "somepath/clipped_and_trimmed_reads/file1_R1.fq.gz",
        "somepath/clipped_and_trimmed_reads/asdf1_R1.fq.gz",
        "somepath/clipped_and_trimmed_reads/zxcv1_R1.fq.gz"

但是随着我们获得更多文件，这变得非常麻烦，而且正如您在问题中指定的那样，这不是您想要的。我们需要编写一个小函数来为我们获取所有文件名。

import glob
import re

data=config["raw_data"]
samples = []
locations = {}
for file in glob.glob(data + "/**", recursive=True):
    if '_R1.fq.gz' in file:
        split = re.split('/|_R1', file)
        filename, directory = split[-2], split[-3]
        locations[filename] = directory  # we will need this one later
        samples.append(filename)

我们现在可以将其提供给我们的规则：

rule all:
    input:
        read1=expand("{output}/clipped_and_trimmed_reads/{sample}_R1.fq.gz", output=config["output"], sample=samples),
        read2=expand("{output}/clipped_and_trimmed_reads/{sample}_R2.fq.gz", output=config["output"], sample=samples)

请注意，我们不再将 sample 作为通配符，但我们 'expand' 将其放入我们的 read1 和 read2 中，从而实现输出和样本的所有可能组合。

然而我们只完成了一半.. 如果我们像这样调用 Snakemake，它会确切地知道我们想要哪个输出，以及哪个规则可以生成这个（规则 clip_and_trim_reads）。然而它仍然不知道去哪里寻找这些文件。幸运的是，我们已经有了一个字典来存储这些（存储在位置）。

rule clip_and_trim_reads:
    input:
        read1=lambda wildcards: expand("{data}/{set}/{sample}_R1.fq.gz", data=config["raw_data"], set=locations[wildcards.sample], sample=wildcards.sample),
        read2=lambda wildcards: expand("{data}/{set}/{sample}_R2.fq.gz", data=config["raw_data"], set=locations[wildcards.sample], sample=wildcards.sample)
    output:
        ...

现在一切正常了！！甚至更好；由于规则 clip_and_trim_reads 的所有结果都写入了一个文件夹，因此从这里继续应该容易得多！

p.s。我没有测试过任何代码，所以第一次尝试时可能并非一切正常。但是消息仍然存在。

Snakemake - 从输入文件动态派生目标

Snakemake - dynamically derive the targets from input files

python

bioinformatics

python-3.x