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

Question

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

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

• 数据集下载，提供 1 个文件（SE 数据）或 2 个文件（PE 数据）-->
• 一组特定于 1 个文件的规则 A 或一 组特定于 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 语句，但我还没有设法解决我的问题。

你有什么提示/建议/方法可以让“切换”发生？

Answer 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

编辑：

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

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

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

确切地！当规则需要trimmed_sample1_R1.fastq.gz 时，Snakemake 如何知道您的样本名称？是样本的名字，sample1，还是sample1_R1？它可以是两者之一，这让蛇形抱怨它不知道如何解决这个问题。当您添加规则顺序时，您告诉 Snakemake，当不清楚时，按此顺序解决。

3. 对齐规则中的命令需要 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无法使用该指令。