Snakemake，RNA-seq：如何根据所分析样本的特征执行管道的一个子部分或另一个子部分？

Posted 2023-02-25

【中文标题】Snakemake，RNA-seq：如何根据所分析样本的特征执行管道的一个子部分或另一个子部分？

【英文标题】：Snakemake, RNA-seq : How can I execute one subpart of a pipeline or another subpart based on the characteristics of the sample that is analysed?

【发布时间】：2020-08-09 07:31:51

【问题描述】：

我正在使用 snakemake 设计一个 RNAseq 数据分析管道。虽然我已经设法做到了，但我想让我的管道尽可能地具有适应性，并使其能够在同一运行分析中处理单读取 (SE) 数据或双端 (PE) 数据，而不是一次分析 SE 数据，另一次分析 PE 数据。

我的管道应该是这样设计的：

提供 1 个文件（SE 数据）或 2 个文件（PE 数据）的数据集下载 --> 针对 1 个文件的规则 A 集

OR

针对 2 个文件的规则 B 集 --> 采用 1 或 2 个输入文件并将其合并的规则 进入单个输出--> 最终规则集。

注意：A 的所有规则都有 1 个输入和 1 个输出，B 的所有规则都有 2 个输入和 2 个输出，它们各自的命令如下所示：

1 输入：somecommand -i input -o output 2 个输入：somecommand -i1 input1 -i2 input2 -o1 output1 -o2 output2

注2：除了输入/输出不同外，集合A和B的所有规则都具有相同的命令、参数/等...

换句话说，我希望我的管道能够根据示例在规则集 A 或规则集 B 的执行之间切换，方法是在开始时在配置文件中提供有关示例的信息 (样本 1 是 SE，样本 2 是 PE……这是事先知道的）或要求 snakemake 在数据集下载后计算文件数，以便为每个样本选择合适的下一组规则。如果您看到其他方法可以做到这一点，欢迎您告诉我们。

我考虑过使用检查点、输入函数和 if/else 语句，但我没有设法解决这些问题。

您有任何提示/建议/方法来实现“转换”吗？

【参考方案1】：

如果您事先知道布局，那么最简单的方法是将其存储在某个变量中，如下所示（或者您从配置文件中将其读取到字典中）：

layouts = "sample1": "paired", "sample2": "single", ... etc

然后你可以做的是像这样“合并”你的规则（我猜你是在谈论修剪和对齐，所以这就是我的例子）：

ruleorder: B > A

rule A:
    input:
        sample.fastq.gz
    output:
        trimmed_sample.fastq.gz
    shell:
        "somecommand -i input -o output"

rule B:
    input:
        input1=sample_R1.fastq.gz,
        input2=sample_R2.fastq.gz
    output:
        output1=trimmed_sample_R1.fastq.gz,
        output2=trimmed_sample_R2.fastq.gz
    shell:
        "somecommand -i1 input.input1 -i2 input.input2 -o1 output.output1 -o2 output.output2"


def get_fastqs(wildcards):
    output = dict()
    if layouts[wildcards.sample] == "single":
        output["input"] = "trimmed_sample2.fastq.gz"
    elif layouts[wildcards.sample] == "paired":
        output["input1"] = "trimmed_sample1_R1.fastq.gz"
        output["input2"] = "trimmed_sample1_R2.fastq.gz"
    return output


rule alignment:  
    def input:
        unpack(get_fastqs)
    def output:
        somepath/sample.bam
    shell:
        ...

这里发生了很多事情。

首先您需要一个规则顺序，这样蛇形才能知道如何处理模棱两可的情况 规则 A 和 B 都必须存在（除非您对输出文件进行了一些修改）。 对齐规则需要一个输入函数来确定它需要哪个输入。

一些自我推销：我做了一个蛇形管道，它可以做很多事情，包括 RNA-seq 和在线下载样本并自动确定它们的布局（单端与双端）。请看看它是否解决了您的问题：https://vanheeringen-lab.github.io/seq2science/content/workflows/rna_seq.html

---

编辑：

当您说“合并”规则时，您是指规则 A、B 和对齐方式吗？

我的措辞不清楚。合并我的意思是“合并 单端和双端和双端逻辑在一起，因此您可以继续使用单个规则（例如计数表，您可以命名它）。

规则顺序：为什么选择 B > A？确保配对样本不会最终以单端规则运行？

没错！当一条规则需要 trimmed_sample1_R1.fastq.gz 时，Snakemake 怎么知道你的样本名称？样品的名称是 sample1，还是 sample1_R1？两者都可以，这让snakemake抱怨它不知道如何解决这个问题。添加规则顺序时，您告诉 Snakemake，如果不清楚，请按此顺序解决。

对齐规则中的命令需要 1 或 2 个输入。我打算在 params 指令中使用 if/else 来选择输入。我这样想对吗？ （我认为您在管道中也这样做了）

是的，我们就是这样解决的。我们这样做是因为我们希望每个规则都有自己的环境。如果您不使用单独的 conda 环境进行对齐，那么您可以这样做更干净/更漂亮，就像这样

rule alignment:  
    input:
        unpack(get_fastqs)
    output:
        somepath/sample.bam
    run:
        if layouts[wildcards.sample] == "single":
            shell("single-end command")
        if layouts[wildcards.sample] == "paired":
            shell("paired-end command")

我觉得这个选项比我们在 seq2science 管道中所做的要清晰得多。但是在 seq2science 管道中，我们支持许多不同的对齐器，并且它们都有不同的 conda 环境，因此不能使用 run 指令。

【讨论】：

感谢您的回答。今天我会仔细看看这一切，我会及时通知你，但这似乎很有希望！ （你猜对了修剪、对齐和诸如此类的东西！）。

好吧，所以我看了你的答案（还有你的管道，真的很棒，这对我有很大帮助！）我需要更多解释：1.当你说“合并”规则，您是指规则 A、B 和对齐吗？或者只是规则 A 和 B ？还是更多规则？ 2.规则顺序：为什么选择B>A？确保配对样本最终不会在单端规则中运行？ 3.对齐规则中的命令需要1或2个输入。我打算在 params 指令中使用 if/else 来选择输入。我这样想对吗？ （我认为您在管道中也这样做了）。

如果是这样，您能否对您的答案进行简短编辑以提及它？这样，你的回答就完美了，我可以关闭这篇文章了！！！非常感谢您的帮助！

@athiebaut 我在编辑中回答了您的问题。乐于助人！

@bli 那些“def”不应该在那里！我删除了它们。解包，“解包”字典，参见：snakemake.readthedocs.io/en/stable/snakefiles/…