2

我正在使用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 语句，但我还没有设法解决这些问题。

您是否有任何提示/建议/方法可以使“转换”发生？

4

1 回答 1

5

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

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 都必须存在（除非你对输出文件做 sth hacky）。
  • 对齐规则需要一个输入函数来确定它需要哪个输入。

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

编辑：

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

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

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

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

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

于 2020-08-09T12:02:40.050 回答