流程构建 Snakemake使用初步

呦呦呦

在学习生信的过程中，总是听说流程啊，管道啊，pipeline啊，到底这些意味着什么？你能把测序数据从ncbi下载下来，这不叫流程；你可以继续把数据做完质控，然后搞个比对，再做做什么差异分析啊、富集分析啊、各种类型数据的联合分析啊，这也不叫流程，谁知道你中间因为某个包安装用了多久的时间，谁晓得你这个过程是不是 reproducible 的？要是你中间出了问题怎么办，可控性如何？监控性如何？要是这些问题都能解答，我想就应该可以叫做一个流程了吧，流程化、自动化、可控化、高复用化。作为一个新人，要独自解决这些功能还是有些困难，肯定需要“假于物”的，恰好那天看到有个朋友说到用snakemake来些流程，哈哈，这不就有了，下面就说一说我初步学习使用的理解吧。

snakemake 实现流程化

在构建所谓的生信分析pipeline的时候，首先要实现pipeline的流程化，那么snakemake是怎么实现的呢？那就要说到snakemake 的rule语法了。如下所示。

1
2
3
4
5
6
7
8
# in python
rule step_name:
	input:
		"something"
	output:
		"something"
	shell:
		"some commands do with input, result is output"

写好rule后，保存名字为 Snakefile。我们可以使用这样的方式，来构建每一个分析过程，这样做的好处就是直接明了，这一个rule是干什么的，需要什么输入，会产生什么输出，怎么产生都描述完全。我通常把do what作为rule的step_name，这样一来把分析过程先进行步骤分解，写成rule的形式，这就先完成了分析过程的步骤化，然后，因为流程是一个步骤接一个步骤的，有个次序的问题，如何实现？通过input和output来实现，以rule1的output作为rule2的input，这样就可以把步骤串联起来了。对于这个问题，snakemake提供了一个feature，叫做inheritance rule。我们看下面的示例，step_one对测序reads进行比对，得到bam文件，step_two要对bam文件进行统计。对于整个构建过程，我们完全可以采用inheritance rule的方式来进行，这样对于最终结果是怎么来的，可以按图索骥。

1
2
3
4
5
6
7
8
9
10
11
12
13
rule mapping:
	input: "rawdata.fastq"
	output: "rawdata.bam"

# 做法一
rule stats:
	input: "rawdata.bam"
	output: "mapping.stats.txt"

# 做法二
rule stats:
	input: rules.mapping.output
	output: "mapping.stats.txt"

除了上述的 rule 形式，snakemake 还提供了一个 rule all 的语法。它是你这个 pipeline 的最终结果集合，对于你需要的结果，你应该全部显式地列在rule all中，它位于所有rule的最开端。

1
2
3
4
5
rule all:
	input:
		"step_1_result",
		"step_2_result",
		...

有人要问了，为什么这么干？原因我认为有二，第一脚本文件的解析过程是由前至后、由上至下的。那么 rule all 放在前面，相当于你的pipeline有个概览。其二，正如前面所述，step-by-step的input和output是串联的，但是整个流程的 step-by-step 不是单一的，而是有多个step-by-step，它们是并联的。而不论串联并联，都有个结果，那就是 rule all。结合实际运行我们会发现，snakemake 的运行方式就是按图(rule all)索骥，我（rule all)需要什么结果 (input)，就去找生成这个结果的rule，然后运行这个rule；对于这个rule，重复这个依据 input 寻找和执行上游 rule 的过程，对于这个寻找过程，我们可以使用 -rn 参数打印出来。如果你的文件名字是 “Snakefile”, 后面的 -s Snakefile 可以省略。

1
snakemake -rn -s Snakefile

snakemake 实现自动化

虽然我尚未正式做过任何项目，但是依据文献提供的原始数据，一个正常的测序项目，样品数目大多数在 10 个及以上。我们不可能把这些样品都一个个写成一个个 rule，就比如比对过程，你执行一次比对命令执行完成一个样品的比对，通常的做法是使用for循环来实现。但这里是 snakemake，它会帮你的，那就是使用通配符 (wildcards)。一个 wildcards 就是一个匹配的文件名的简单正则。对于每个匹配的文件，snakemake 都会以文件作为 input 执行一次 rule， 看下面的示例。

1
2
3
4
5
6
7
8
# samples => A1.fastq A2.fastq, A3.fastq
rule mapping:
	input:
		"{sample}.fastq"
	output:
		"{sample}.bam"
	shell:
		"bwa mem .... {input} -o {output}"

这样，对于每个 sample 都会进行一次比对。这样就毋须我们自己写 for 循环了。有时候，我们有很多样品，但是，我们仅仅想使用其中某几个作为输入，该如何？ snakemake 提供了constraint wildcards，来对 wildcards 匹配的样品进行限制。

1
2
3
4
5
# samples => A1.fastq A2.fastq, A3.fastq, A14.fastq, A15.fastq
rule mapping:
	input:
		"{sample, \w\d}.fastq"
	...

在我们加入“\w\d”这个正则之后，sample将不会匹配 A14.fastq 和 A15.fastq，这是因为 snakemake 在解析的时候，默认会把 “{sample}.fastq” 解析成 “*.fastq”；但是如果我们使用 constraint wildcards的话，那 “{smaple, \w\d}.fastq” 就会被解析成 “\w\d.fastq”, 这就自然不会匹配 A14.fastq 这样的样品了。除了在 rule 中使用 constraint wildcards，你还可以进行全局声明，不过你得先声明后使用，放在所有的 rule 前面。

1
2
wildcard_constaints:
    sample = "\w\d"

对于上面的那个比对示例，有同学会问了，你这个看起来只是 single end 的测序结果作为输入，那要是 paired end 的测序结果，该怎么输入呢？欸，这里就有两种方式啦，一种是，我再提供一个咯；另一种就是使用 snakemake 提供的 expand 函数。但是，这里问题就来了。

1
2
3
4
5
6
7
8
9
# 普通方式：
rule mapping:
	input:
		"A1_read1.fastq",
		"A1_read2.fastq"
# expand
rule mapping:
	input:
		expand("A1_{read}.fastq", read=["read1", "read2"])

普通的expand是不支持 wildcards 的，看它的使用方式就知道了，与 str.format(read=[]) 形式是一致的。但是，如果要使用 wildcards 怎么办呢？我们需要定义一个 helper，以它作为 input, 这个 helper 你可以正常函数定义，也可以使用匿名函数（如果仅使用一次的话）, 需要注意的是，这个 helper 的参数只能有一个，那就是 wildcards。

1
2
3
rule mapping:
	input:
		lambda wildcards: expand("{sample}_{read}.fastq", sample=wildcards.sample, read=["read1", "read2"])

snakemake 实现流程的可控化

在分析的过程中，分析的步骤辣么多，时间辣么长，如何保证分析过程的正常进行？出现错误如何发现？是不是一步错，步步错，全局失败？是不是只能一水儿的执行下去？问题这么多，snakemake 说看我的。为了实现流程的可控化，snakemake 提供了多个 feature。要控制整个流程，需要控制什么？没错，那就是控制每一个步骤，控制每一个 rule。如何控制？记录运行日志，写配置。尽管 snakemake 自己会生成日志文件，对于每一个 rule，你也可以指定日志生成，这个通过 rule.log 实现。每个 rule 运行命令的控制可以通过 rule.params 来提供命令的参数选项进行控制。下面看示例。

  • 命令参数的可控
1
2
3
4
5
6
7
8
9
10
rule mapping:
	input:
		lambda wildcards: expand("{sample}_{read}.fastq", sample=wildcards.sample, read=["read1", "read2"])
	log:
		"mapping/{sample}.log"
	threads: 8
	params:
		extra = ""
	shell:
		"bwa mem {params.extra} -t {threads} -o {output} {input} > {log}"

这里我们指定了 “bwa mem” 的运行线程数为 8，日志输出到 “mapping/{sample}.log”(还记得 wildcards 吧？) 中，额外的参数你可以写入到 params.extra 中。对于 线程数，虽然可以通过 params 来指定，但是呢，没有在 threads 声明的话，即使指定了，snakemake 仍旧会把它降级为一核的，你还要在运行 snakemake 的时候指定 “–cores 10”之类的参数。到这里你可能会问，可控性仅此吗？当然不是了。我们继续看。

  • 执行 rule 的可控

如前所述，运行过程不会那么一帆风顺的，肯定会出错。根据我们的日志文件，我们定位出错的 rule。在 debug 之后，那肯定要测试一下。这时候以下参数就十分有用了，它可以让我们仅仅执行那些发生变动的 rule。

1
2
3
4
5
6
snakemake -n 									# 打印出运行流程，检验 input、output 
snakemake -p 									# 打印出解析后所有要运行的命令
snakemake -f mapping 							# 仅仅执行 rule mapping及其所有依赖 rule
snakemake -f `snakemake --list-code-changes` 	# 执行所有代码发生变动的 rule
snakemake -f `snakemake --list-input-changes` 	# 执行所有输入发升变动的 rule
snakemake -f `snakemake --list-params-changes` 	# 执行所有参数设定发生变动的 rule
  • 执行 rule 环境的可控

做生信最头疼的东西是什么？配置环境，安装包啊，大批的新手倒在这里，白骨累累。来来来，snakemake 助你披襟斩棘。snakemake 提供的conda 和 wrappers 给你飞一般的感觉。对于前者，我们可以提供一个环境配置文件如 “env.baw.yaml”, 如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
# env.bwa.yaml
channels:
 - bioconda
 - conda-forge
dependencies:
 - bwa = 2.7.1
 - samtools = 2.5.8

# your rule
rule mapping:
	# input, output, ...
	conda:
		"env.bwa.yaml"

然后在启动 snakemake 时候，添加 “–use-conda” 参数，snakemake 将会自动生成 conda 环境，并基于 “env.bwa.yaml” 的配置安装依赖。你以为这就完了吗？NO NO NO，你还有wrappers。对于一些常用的软件，snakemake 把它整合成wrapper，按照 wrapper 的使用指南配置参数之后，snakemake 会为他自动配置运行环境，你什么都不用管, 不足之处就是支持的软件有点少。

1
2
3
4
5
6
7
8
9
10
11
12
13
rule mapping:
	# input, output, ...
	params:
		index = "genome",
		extra = r"-R '@RG\tID:{sample}\tSM:{sample}'",
		sort = "samtools", # can be 'none', 'picard'
		sort_order = "coordinate", # can be 'queryname'
		sort_extra = "" # extra args for samtools/picard
	threads: 8
	log:
		"logs/mapping/{sample}_bwa_mem.log"
	wrapper:
		"0.27.1/bio/bwa/mem"
  • 执行 rule 脚本的可控性

在不使用 wrapper 的时候，我们需要使用 rule.shell 来指定运行的命令，但是分析过程中，我们不仅会用到 bash，还会用到 python，甚至 R，snakemake 也统统让你使用。因为 snakemake 是基于 python 开发的，你的 Snakefile 其中的 python 代码会被 正常解析。这是其一，其二，如果你要在 rule 中使用 python，可以通过 rule.run 来实现。

1
2
3
4
5
6
7
8
rule mapping:
	# input, output,...
	run:
		print(input)
		print(output)
		print(params)
		...
		shell("some commands")

rule.run 可以直接引导 rule 中的其他参数，如 output，input等等，shell是 snakemake引入的 (自己引入：from snakemake.shell import shell)。除此之外，如果你要使用 R 呢？怎么办？snakemake 提供了 rule.script 来满足你。它可以助你引入外部脚本，并执行它，对于 python 脚本，你可以通过 snakemake.input 的方式引用 rule 的各个选项，对于 R 脚本，你可以使用 snakemake@input，如下所示。

1
2
3
4
5
6
7
8
9
10
# Snakemake
rule mapping:
	# input, output
	script:
		"mapping.R"

# mapping.R
in_file = snakemake@input # as a list
out_file = snakemake@output # as a list
dosomething()

summary

在写好这些 rule 之后，你一运行 snakemake -n，应该会出现一些 error （哈哈）。MissingInputError 表明你的 rule.all.input 中有些结果你没有生成哦，需要看看是那个输出结果你写了，但是没有 rule 生成它。Wildcards Error，表明你的通配符有问题，最直观的就是通配符要匹配的值你没有给。这个值得说一下。看下方的示例。

1
2
3
4
5
6
7
8
9
rule all:
	input:
		"{sample}.bam"

rule mapping:
	input:
		"{sample}.fastq"
	output:
		"{sample}.bam"

初看之下，rule.all 中我们需要 “{sample}.bam” 文件，所以需要运行 rule.mapping。然而，snakemake 会给你报错的，说是无法指定确定的 input。为什么呢？通配符不能无中生有，它通配的是什么？通配的是 pipeline 最后结果，这个最后结果是我们指定的。所以 rule.all.input 不能含有不确定的东西，即必须确定。所以首先，rule.all.input应该改成：

1
2
rule all:
	input: expand("{sample}.bam", sample=["A1", "A2", "A3", "A4"])

如此一来，就可以了。但是 wildcards 的使用还有个问题，那就是 rule.input, rule.output, rule.log 这三个的通配符要一致。就是说，当 input 有通配符时，表明该 rule 会执行多次，那么 output 也应该是有多个，log 也应该有多个，不能覆盖，所以都需要一致的通配符 wildcards。但是，有时候，我们不想把 “{sample}.bam” 作为 rule.all 的 input （即最终结果），但仍旧想使用 wildcards，上面的修改方式就还会出现问题。例如：

1
2
3
4
5
6
7
8
9
10
11
12
rule all:
	input: "mapping.stats.txt"

rule mapping:
	input:
		"{sample}.fastq"
	output:
		"{sample}.bam"

rule stats:
	input: "{sample}.bam"
	output: "mapping.stats.txt"

这里还是一样，无法确定通配符要通配的对象。这时候，就要像 rule.all.input 一样，在某一处，确定通配符通配对象，结合前面所说，input/output/log 三种的通配情况应该一致，我们如果要在 rule.mapping 中使用通配，那它就不能发生变化，所以应该改为：

1
2
3
4
# ... 
rule stats:
	input: expand("{sample}.bam", sample=["A1", "A2", "A3", "A4"])
	# ...

如果 rule.mapping 中不使用 wildcards 的话，就可以把 rule.mapping.input/output/log 全部改成 expand 的形式。但是请注意，上述代码仅仅是作为解决通配符问题的例子，如果是比对的话，最好使用通配符，因为比对软件的输入输出是一对一的，不像 fastqc，一个命令输入多个样品，生成多个样品的结果。

最后，snakemake 还有许多 features，比如 temp/protected, cluster, derectory, touch等特性，值得再了解了解。

参考

Snakemake

Comments

▲