Snakemake 可变数量的文件

Question

我处于一种情况，我想将我的工作流程分散到可变数量的块中，我事先不知道。也许通过具体来解释问题最容易：

有人递给我使用bcl2fastq该no-lane-splitting选项解复用的 FASTQ 文件。我想根据车道拆分这些文件，单独映射每个车道，然后最后再次收集所有内容。但是，我事先不知道车道数。

理想情况下，我想要这样的解决方案，

rule split_fastq_file: (...)  # results in N FASTQ files
rule map_fastq_file: (...)  # do this N times
rule merge_bam_files: (...)  # merge the N BAM files

但我不确定这是可能的。该expand功能要求我知道车道数，并且无法看到如何为此使用通配符。

我应该说我对 Snakemake 比较陌生，而且我可能完全误解了 Snakemake 的工作原理。我花了一些时间来习惯于通过关注输出文件然后倒退来思考“颠倒”的事情。

Answer 1

一种选择是checkpoint在拆分 fastqs 时使用，以便您可以在稍后动态重新评估 DAG 以获得结果车道。

这是一个 MWE 一步一步：

• 设置并制作一个示例 fastq 文件。

# Requires Python 3.6+ for f-strings, Snakemake 5.4+ for checkpoints
import pathlib
import random

random.seed(1)

rule make_fastq:
    output:
        fastq = touch("input/{sample}.fastq")

• 创建 1 到 9 之间的随机数通道，每个通道具有 1 到 9 的随机标识符。请注意，我们将其声明为checkpoint，而不是规则，以便我们稍后可以访问结果。此外，我们将此处的输出声明为特定于示例的目录，以便我们稍后可以将其放入其中以获取创建的通道。

checkpoint split_fastq:
    input:
        fastq = rules.make_fastq.output.fastq
    output:
        lane_dir = directory("temp/split_fastq/{sample}")
    run:
        pathlib.Path(output.lane_dir).mkdir(exist_ok=True)
        n_lanes = random.randrange(1, 10)-
        lane_numbers = random.sample(range(1, 10), k = n_lanes)
        for lane_number in lane_numbers:
            path = pathlib.Path(output.lane_dir) / f"L00{lane_number}.fastq"
            path.touch()

• 做一些中间处理。

rule map_fastq:
    input:
        fastq = "temp/split_fastq/{sample}/L00{lane_number}.fastq"
    output:
        bam = "temp/map_fastq/{sample}/L00{lane_number}.bam"
    run:
        bam = pathlib.Path(output.bam)
        bam.parent.mkdir(exist_ok=True)
        bam.touch()

• 为了合并所有处理过的文件，我们使用一个输入函数来访问在 中创建的通道split_fastq，以便我们可以expand对这些进行动态处理。我们expand在中间处理步骤链中的最后一条规则上执行，在本例中为map_fastq，以便我们要求正确的输入。

def get_bams(wildcards):
    lane_dir = checkpoints.split_fastq.get(**wildcards).output[0]
    lane_numbers = glob_wildcards(f"{lane_dir}/L00{{lane_number}}.fastq").lane_number
    bams = expand(rules.map_fastq.output.bam, **wildcards, lane_number=lane_numbers)
    return bams

• 这个输入函数现在让我们可以轻松访问我们希望合并的 bam 文件，无论有多少，也不管它们可能被调用。

rule merge_bam:
    input:
        get_bams
    output:
        bam = "temp/merge_bam/{sample}.bam"
    shell:
        "cat {input} > {output.bam}"

本示例运行，并与random.seed(1)发生创建三个车道（l001，l002，和l005）。

如果您不想使用checkpoint，我认为您可以通过创建一个输入函数来实现类似的功能，merge_bam该函数打开原始输入 fastq，扫描读取名称以获取车道信息，并预测输入文件应该是什么。然而，这似乎不太稳健。