snakemake - Snakemake,RNA-seq:如何根据所分析样本的特征执行管道的一个子部分或另一个子部分?
问题描述
我正在使用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 语句,但我还没有设法解决这些问题。
您是否有任何提示/建议/方法可以使“转换”发生?
解决方案
如果您事先知道布局,那么最简单的方法是将其存储在某个变量中,如下所示(或者您从配置文件中将其读取到字典中):
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
编辑:
- 当您说“合并”规则时,您是指规则 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
不能使用该指令。