带有混合通配符的输出文件

Question

我不知道如何用 Snakemake 做到这一点。

首先，说我的“全部”规则是：

rule all:
    input: "SA.txt", "SA_T1.txt", "SA_T2.txt",
        "SB.txt", "SB_T1.txt", "SB_T2.txt", "SB_T3.txt"

注意SA 有两个_T# 文件，而SB 有三个这样的文件，这是其中的关键元素。

现在我想写一个这样的规则来生成这些文件：

rule X:
    output: N="S{X}.txt", T="S{X}_T{Y}.txt"
    (etc.)

但 SnakeMake 要求两个输出模板具有相同的通配符，而这些模板没有。此外，即使 SnakeMake 可以处理多个通配符，它​​也可能希望找到与 S{X}_T{Y}.txt 模板匹配的单个文件名，但我希望它匹配 {X} 与第一个模板的 {X} 匹配的所有文件，即我希望output.T 是一个列表，而不是单个文件。所以看起来这样做的方法是：

def myExpand(T, wildcards):
    T = T.replace("{X}", wildcards.X)
    T = [T.replace("{Y}", S) for S in theXYs[wildcards.X]]
    return T

rule X:
    output: N="S{X}.txt", T=lambda wildcards: myExpand("S{X}_T{Y}.txt", wildcards)
    (etc.)

但我不能这样做，因为 lambda 函数不能用于输出部分。

我该怎么做？

在我看来，这支持在输出语句上支持 lambda 函数，提供一个通配符字典，其中填充来自输出语句的已解析部分的值。

---

响应 cmets 的添加：

需要通配符 Y 的值，因为其他规则对具有通配符 Y 的文件有输入。

我的规则知道 Y（和 X）的不同值，它需要处理从数据库读取到 Python 字典的数据。

X 有很多值，每个 X 值有 2 到 6 个 Y 值。我认为对 X 的每个值使用单独的规则是没有意义的。但是，我可能错了，因为我最近了解到可以将规则放入循环中，并创建多个规则。

有关工作流程的更多信息：我正在将来自一个人的多个肿瘤样本的体细胞变异 VCF 文件合并到一个 VCF 文件中，并这样做使得对于任何一个肿瘤中的每个被调用的变异，所有肿瘤都不会调用该变异分析以确定变体的读取深度，该变量包含在合并的 VCF 文件中。

整个过程涉及大约 14 个步骤，可能多达 14 条规则。我其实并不想使用 14 条规则，但更愿意在一条规则中完成所有操作。

但是，我现在认为解决方案确实是使用许多单独的规则。我之所以避免这样做，部分原因是因为有大量不必要的中间文件，但实际上，这些临时文件无论如何都存在，在单个大规则内。使用多个规则，我可以将它们标记为 temp()，这样 Snakemake 将在最后删除它们。

为了充实这个我认为是合法的讨论，让我们假设一个可能出现的简单情况。假设对于许多人中的每一个人，您有 N (>=2) 个肿瘤 VCF 文件，就像我一样，并且您想编写一个规则来生成 N+1 个输出文件，每个肿瘤一个输出文件加一个与此人关联的更多文件。使用通配符 X 表示人员 ID，使用通配符 Y 表示人员 X 中的肿瘤 ID。假设操作是将所有肿瘤 VCF 文件中存在的所有变体放入人员输出 VCF 文件中，并将所有其他变体放入相应的肿瘤输出文件中它们出现的肿瘤。 假设一个程序从 N 个输入文件生成所有 N+1 个文件。你是怎么写规则的？

你想要这个：

rule:
    output: COMMON="{X}.common.vcf", INDIV="{X}.{Y}.indiv.vcf"
    input: "{X}.{Y}.vcf"
    shell: """
        getCommonAndIndividualVariants --inputs {input} \
        --common {output.COMMON} --indiv {output.INDIV}
        """

但这违反了输出通配符的规则。

我这样做的方式是使用两个规则，虽然不太令人满意，但很有效，第一个规则的输出模板具有更多通配符，第二个模板具有更少的通配符通配符，并让第二条规则创建临时输出文件，这些文件由第一条规则重命名为最终名称：

rule A:
    output: "{X}.{Y}.indiv.vcf"
    input: "{X}.common.vcf"
    run: "for infile in {input}: os.system('mv '+infile+'.tmp'+' '+infile)"

rule B:
    output: "{X}.common.vcf"
    input: lambda wildcards: \
        expand("{X}.{Y}.vcf", **wildcards, Y=getYfromDB(wildcards["X"]))
    params: OUT=lambda wildcards: \
        expand("{X}.{Y}.indiv.vcf.tmp", Y=getYfromDB(wildcards["X"]))
    shell: """
        getCommonAndIndividualVariants --inputs {input} \
        --common {output} --indiv {params.OUT}
        """

Answer 1

我不知道您的工作流程的其余部分如何，最好的解决方案是取决于上下文。

1

如何将规则分成两部分，一个创建"SA.txt", "SA_T1.txt", "SA_T2.txt"，另一个创建"SB.txt", "SB_T1.txt", "SB_T2.txt", "SB_T3.txt"？

2

另一种可能性是在输出指令中只包含 {X} 个文件，但让规则创建其他文件，即使它们不在输出指令中。如果 {Y} 文件是 DAG 的一部分，则此方法不起作用。

3（最佳解决方案？）

第三种可能是最好的解决方案可能是在 X 规则和需要 X 输出的规则中使用聚合通配符。

那么解决办法就是

rule X:
    output: N="S{X_prime}.txt", T="S{Y_prime}.txt"

需要这些文件的规则如下所示：

rule all:
    input:
        expand("S{X_prime}", X_prime="A_T1 A_T2".split()),
        expand("S{Y_prime}", Y_prime="B_T1 B_T2 B_T3".split())

如果这不符合您的要求，我们可以进一步讨论:)

附言。您可能需要使用 wildcard_constraints 来消除规则 X 的输出歧义。

list_of_all_valid_X_prime_values = "A_T1 A_T2".split()
list_of_all_valid_Y_prime_values = "B_T1 B_T2 B_T3".split()

wildcard_constraints:
    X_prime = "({})".format("|".join(list_of_all_valid_X_prime_values))
    Y_prime = "({})".format("|".join(list_of_all_valid_Y_prime_values))

rule all:
    ...

Answer 2

我的理解是，snakemake 的工作步骤如下：

• 查看要求生成的文件的名称。
• 查找在其输出中具有匹配模式的规则。
• 根据上述匹配推断通配符的值。
• 使用通配符确定所选规则的输入名称应该是什么。

您的规则无需输入即可生成其输出，因此推断通配符Y 的值的问题并不明显。

您的规则如何知道它需要使用多少个不同的 Y 值？

如果您找到一种方法来确定Y 的值，只知道X 的值以及预定义的“python 级”函数和变量，那么可能有一种方法可以将Y 作为内部变量您的规则，而不是通配符。

在这种情况下，工作流只能由文件 S{X}.txt 驱动。 S{X}_T{Y}.txt 只是规则执行的副产品，而不是其显式输出。