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

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 函数，提供一个通配符字典，其中填充了来自输出语句的已解析部分的值。

---

回复评论的补充:

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

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

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

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

整个过程涉及大约 14 个步骤，可能多达 14 个规则。我实际上不想使用 14 条规则，而是更喜欢在一个规则中完成所有操作。

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

为了充实这个我认为是合理的讨论，让我们假设一个可能出现的简单情况。假设对于许多人中的每一个人，你有 N (>=2) 个肿瘤 VCF 文件，就像我一样，并且你想编写一个规则来生成 N+1 个输出文件，每个肿瘤一个输出文件加上一个与此人关联的更多文件。对人员 ID 使用通配符 X，对人员 X 中的肿瘤 ID 使用通配符 Y。假设该操作是将所有肿瘤 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:
    ...