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

Question

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

中分析 PE 数据。

我的流水线应该是这样设计的：

提供1个文件（SE数据）或2个文件（PE数据）的数据集下载-->
一组规则 A 特定于 1 个文件 OR 一组规则 B 特定于 2 个文件 -->
采用 1 或 2 个输入文件并合并的规则 it/them 到单个输出 -->
最终规则集。

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

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

注2：除了inputs/outputs不同外，A组和B组的所有规则命令都是一样的，parameters/etc...

换句话说，我希望我的管道能够根据示例在执行规则集 A 或规则集 B 之间切换，方法是在配置文件中提供有关示例的信息开始（样本 1 是 SE，样本 2 是 PE...这是事先已知的）或要求 snakemake 在数据集下载后计算文件数量，以便为每个样本选择正确的下一组规则。如果您看到另一种方法，欢迎您提出来。

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

你有什么 hints/advice/ways 来实现这种“转变”吗？

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:
        ...

这里发生了很多事情。

首先你需要一个规则顺序，这样 snakemake 就知道如何处理不明确的情况
规则 A 和 B 都必须存在（除非您对输出文件做某事）。
对齐规则需要一个输入函数来确定它需要哪个输入。

Some self-promotion：我制作了一个 snakemake 管道，它可以做很多事情，包括 RNA-seq 和在线下载示例并自动确定它们的布局（single-end vs paired-end ).请看看它是否解决了您的问题：https://vanheeringen-lab.github.io/seq2science/content/workflows/rna_seq.html

编辑:

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

我的措辞不清楚。通过合并，我的意思是“合并 single-end 和 paired-end 以及 paired-end 逻辑在一起，因此您可以继续使用单个规则（例如，计数 table，您可以命名）。

规则顺序：为什么选择 B > A？确保配对样本不会在 single-end 规则中以运行结束？

没错！当规则需要 trimmed_sample1_R1.fastq.gz 时，Snakemake 怎么知道你的样本的名字？样本的名称是 sample1，还是 sample1_R1？它可以是任何一个，这使得 snakemake 抱怨它不知道如何解决这个问题。添加ruleorder的时候告诉Snakemake，不清楚的时候，按照这个顺序解决。

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

是的，我们就是这样解决的。我们这样做是因为我们希望每条规则都有自己的环境。如果你不使用单独的 conda 环境进行对齐，那么你可以这样做 cleaner/prettier，像这样

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 指令。

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?

snakemake

rna-seq