如何在 snakemake 中处理可变数量的重复

Question

我正在处理大量的 fastq 文件。我有很多文件夹，每个样本一个，每个包含多个重复（repl），例如一个文件夹可能包含：

Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R{r}_{repl}.fastq.gz
{r} = 1 or 2
{repl} = [1...n+1] is a three digit integer with leading zeros

每个文件夹都有不同数量的文件，例如例如，n 可能是 13，在这种情况下，将有以下文件：

Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_001.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_002.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_003.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_004.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_005.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_006.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_007.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_008.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_009.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_010.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_011.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_012.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_013.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R1_014.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_001.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_002.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_003.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_004.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_005.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_006.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_007.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_008.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_009.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_010.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_011.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_012.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_013.fastq.gz
Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001_R2_014.fastq.gz

...

我需要在我的 shell 命令 r1 和 r2 中传递给鲑鱼，其中 r1 是文件名中包含 'R1' 的样本的文件列表，同样适用于 R2。

因为我不知道每个文件夹中可能会有多少副本，所以我想检查文件是否存在，return 我自己的函数 repl_expand 中的列表。然而，它会生成所有样本中所有文件的列表，而不是受限于每个样本，因为它是从 expand(["{sample}"], sample=DATASETS).

传递的列表。

似乎我可以将“{sample”} 传递给我的函数，但这只是一个字符串...用实际样本名称替换 {sample} 显然是由 snakemake 在调用我的函数之后，但如果我使用 snakemake 扩展函数，则在调用之前。

需要明确的是，我希望鲑鱼运行三次，每次 r1 都是示例的 'R1' 文件列表，对于 r2 也是如此。

因此，我希望能够在三个调用中的每一个中仅将单个迭代示例名称传递给我的函数，以便为该示例的鲑鱼命令行生成 r1 和 r2。

这是我的 Snakefile，它确实产生了三个鲑鱼运行，但它们都是一样的，包含所有样本的数据。

import os

MAX_REPLICATES = 30
DATASETS = ["Sample_AJL-0704-ME/AJL-0704-ME_GCCAAT_L001",
 "Sample_BG-1058-ME/BG-1058-ME_ACAGTG_L002",
 "Sample_BO-0712-ME/BO-0712-ME_GTGAAA_L002"]
# and many more...

def repl_expand(templates, direction):
    ss = []
    for templ in templates:
        templ = templ + "_R{dir}_{repl:03d}.fastq.gz"
        s = []
        for i in range (1, MAX_REPLICATES):
            name = templ.format(dir=direction, repl=i)
            if os.path.exists(templ.format(dir=direction, repl=i)):
                s.append(name)
        ss.append(s)
    return ss

SALMON = "/path/to/salmon/salmon-1.4.0_linux_x86_64/salmon-latest_linux_x86_64/bin/salmon"

rule all:
    input: expand(["quants/{dataset}.sf"], dataset=DATASETS)


rule salmon_quant:
    input:
        r1 = repl_expand(expand(["{sample}"], sample=DATASETS), 1),
        r2 = repl_expand(expand(["{sample}"], sample=DATASETS), 2),
        index = "path/to/salmon/salmon_sa_index"
    output:
        "quants/{sample}.sf"
    params:
        dir = "quants",
        libtype = 'A',
        threads = 12
    shell:
        "{SALMON} quant \
         -i {input.index} \
         --libtype {params.libtype} \
         --threads {params.threads} \
         --validateMappings \
         --gcBias \
         -o {params.dir} \
         -1 {input.r1} -2 {input.r2}"

我是 snakemake 的新手，所以可能完全错误地处理了这个问题，很高兴能走上正确的道路...

谢谢

Answer 1

首先要注意：您不需要向 expand 函数提供列表。此函数采用模板字符串和 returns 字符串列表作为变量替换的产物。所以规则 all 应该是：

rule all:
    input: expand("quants/{dataset}.sf", dataset=DATASETS)

下一句是关于通配符命名的。但是 Snakemake 不关心你是否在不同的规则中用不同的名字调用同一个实体，我建议你保持一致和明确以避免错误，所以 rule salmon_quant 的输出可能听起来像：

rule salmon_quant:
    output:
        "quants/{dataset}.sf"

有趣的是，您在 rule salmon_quant 中使用名称“{sample}”是在两种不同的上下文中，因此它们的含义完全不同。 output 部分（我重命名的部分）中的 {sample} 是通配符，而 input[=34= 中的 {sample} ] 是 expand 函数的变量。这就是为什么您观察到该替换“异常”的原因。

如果您希望通配符与输出部分中的通配符相匹配，您应该这样做：

rule salmon_quant: input: r1 = lambda wildcards: repl_expand([wildcards.dataset], 1), r2 = lambda wildcards: repl_expand([wildcards.dataset], 2), index = "path/to/salmon/salmon_sa_index" output: "quants/{dataset}.sf"

实际上，如果您简化了 repl_expand 的实现并删除了不需要的循环 for templ in templates:。
，那么您就不需要那个只有一个元素的丑陋列表了

如何在 snakemake 中处理可变数量的重复

How to handle variable numbers of replicates in snakemake

python

snakemake

fastq