小白手册：使用Snakemake来轻松搞定生信分析流程

常常做生物信息分析的同学，应该对pipeline不陌生，一般大家都会通过Shell或者Python来搭建自己的一套分析流程。但当样本较多的时候，人为控制和任务调度耗时就比较多了；而且，如果流程中断出错，排查bug也很麻烦，费时费力。

使用一个好的流程工具，可以让你事半功倍，也有利于别人重复自己的数据分析结果。这里为大家介绍一个生信流程搭建神器：Snakemake。

1.Snakemake的优势

Snakemake是Bioconda的作者——算法工程师Johannes Köster为了数据分析能够更有效的重复，用Python3写的工作流管理工具（Workflow Management System）。这个工具的特点是可以根据文件的依赖关系来智能管理命令的并行和串行。

Snakemake的优势：

• 工作流通过基于Python的语言进行描述，容易理解。前后生成文件逻辑清晰，方便维护。
• Snakemake工作流程可以无缝扩展到工作站、集群和云环境而无需修改工作流定义。作业调度方便控制，比如CPU内核、内存等。
• Snakemake可以使用Conda或Singularity自动部署工作流所需的软件依赖性。
• 如果任务中断，当前的文件状态保留，已成功生成的文件不会被覆盖，已运行完的命令不会重复运行。

2.Snakemake的使用教程

（1）安装Snakemake

Snakemake是一个Python3的包，因此需要提前安装Python3。然后用conda或者pip3安装Snakemake。

1. #conda安装

2. conda install -c bioconda -c conda-forge snakemake

4. #pip3安装

5. pip3 install snakemake

（2）准备数据

用如下代码建一个工作目录snakemake-example，并准备虚拟的空文件数据用于测试。

cd $HOME
# 创建一个工作目录
mkdir snakemake-example
cd snakemake-example
# 虚拟基因组数据
touch genome.fa
# 虚拟fastq数据
mkdir fastq
touch fastq/Sample1.R1.fastq.gz fastq/Sample1.R2.fastq.gz
touch fastq/Sample2.R1.fastq.gz fastq/Sample2.R2.fastq.gz

（左右滑动查看完整代码）

当前环境下的所有文件

1. ├── fastq

2. │   ├── Sample1.R1.fastq.gz

3. │   ├── Sample1.R2.fastq.gz

4. │   ├── Sample2.R1.fastq.gz

5. │   └── Sample2.R2.fastq.gz

6. └── genome.fa

（3）创建Snakefile

要创建一个Snakemake工作流，就需要一个Snakefile，通过定义一系列的规则，来描述如何从输入文件得到输出文件。Snakefile本质上是代码扩展的Python脚本。

在snakemake-example目录下创建Snakefile，内容如下。

SAMPLES = ['Sample1', 'Sample2']
rule all:
    input:
        expand('{sample}.txt', sample=SAMPLES)
rule quantify_genes:
    input:
        genome = 'genome.fa',
        r1 = 'fastq/{sample}.R1.fastq.gz',
        r2 = 'fastq/{sample}.R2.fastq.gz'
    output:
'{sample}.txt'
    shell:
'echo {input.genome} {input.r1} {input.r2} > {output}'

一个规则一般由这几部分组成：

• Inputs：输入文件，依赖文件。

• Outputs：输出文件。

• An action：运行的命令。可以是Bash命令，Python脚本，R脚本，R markdown等。

（4）Snakefile逐行解析

1. SAMPLES = ['Sample1', 'Sample2']

SAMPLES是一个字符列表，存入我们的样本名称。

rule all:
    input:
        expand('{sample}.txt', sample=SAMPLES)

rule all的input是流程最终的目标文件，类似于GNU Make，从顶部指定目标。

expand是Snakemake的特有函数，类似列表推导式。

expand('{sample}.txt',sample=SAMPLES)相当于['{sample}.txt'.format(sample=sample) for sample in SAMPLES]。

即我们目标文件：Sample1.txt和Sample2.txt。

1. rule quantify_genes:

2. input:

3. genome = 'genome.fa',

4. r1 = 'fastq/{sample}.R1.fastq.gz',

5. r2 = 'fastq/{sample}.R2.fastq.gz'

6. output:

7. '{sample}.txt'

8. shell:

9. 'echo {input.genome} {input.r1} {input.r2} > {output}'

对于每个目标文件和中间文件，我们需要创建规则，来定义如何从输入文件创建它们。Snakemake通过匹配文件名来确定规则依赖性。输入和输出文件可以包含多个命名通配符（wildcards）。

rule quantify_genes定义了如何得到我们的目标文件Sample1.txt和Sample2.txt。

规则里的output写的{sample}.txt，其中{sample}就是通配符。

Snakemake的工作流程是这样的：

1、先读rule all，了解目标文件是Sample1.txt和Sample2.txt。

2、Snakemake去找哪个rule的output包含目标文件。

3、当Snakemake发现，rule quantify_genes的 output， {sample}.txt可以匹配目标文件 Sample1.txt和 Sample2.txt。它会把 {sample}分别替换成 Sample1和Sample2，对应两个样本，分别生成两个子流程。每个子流程中，出现 {sample}的地方替换回相应的样本名。

（5）snakemake-np

用snakemake-np打印出实际命令来看一下。

参数-n为dry run，即不实际运行；-p为打印命令；-np即为只打印全部命令而不运行。

$ snakemake -np
Building DAG of jobs...
Job counts:
    count   jobs
1   all
2   quantify_genes
3
[TueFeb1117:46:072020]
rule quantify_genes:
    input: fastq/Sample1.R1.fastq.gz, genome.fa, fastq/Sample1.R2.fastq.gz
    output: Sample1.txt
    jobid: 1
    wildcards: sample=Sample1
echo genome.fa fastq/Sample1.R1.fastq.gz fastq/Sample1.R2.fastq.gz > Sample1.txt
[TueFeb1117:46:072020]
rule quantify_genes:
    input: fastq/Sample2.R1.fastq.gz, genome.fa, fastq/Sample2.R2.fastq.gz
    output: Sample2.txt
    jobid: 2
    wildcards: sample=Sample2
echo genome.fa fastq/Sample2.R1.fastq.gz fastq/Sample2.R2.fastq.gz > Sample2.txt
[TueFeb1117:46:072020]
localrule all:
    input: Sample1.txt, Sample2.txt
    jobid: 0
Job counts:
    count   jobs
1   all
2   quantify_genes
3
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

可以看到两个quantify_genes的子流程，分别得到 Sample1.txt和 Sample2.txt。wildcards分别显示为Sample1和Sample2。

（6）执行命令snakemake

1. $ snakemake

2. Building DAG of jobs...

3. Using shell: /usr/bin/bash

4. Provided cores: 1(use --cores to define parallelism)

5. Rules claiming more threads will be scaled down.

6. Job counts:

7. count jobs

8. 1 all

9. 2 quantify_genes

10. 3

12. [TueFeb1118:01:042020]

13. rule quantify_genes:

14. input: genome.fa, fastq/Sample2.R1.fastq.gz, fastq/Sample2.R2.fastq.gz

15. output: Sample2.txt

16. jobid: 1

17. wildcards: sample=Sample2

19. [TueFeb1118:01:052020]

20. Finished job 1.

21. 1 of 3 steps (33%) done

23. [TueFeb1118:01:052020]

24. rule quantify_genes:

25. input: genome.fa, fastq/Sample1.R1.fastq.gz, fastq/Sample1.R2.fastq.gz

26. output: Sample1.txt

27. jobid: 2

28. wildcards: sample=Sample1

30. [TueFeb1118:01:052020]

31. Finished job 2.

32. 2 of 3 steps (67%) done

34. [TueFeb1118:01:052020]

35. localrule all:

36. input: Sample1.txt, Sample2.txt

37. jobid: 0

39. [TueFeb1118:01:052020]

40. Finished job 0.

41. 3 of 3 steps (100%) done

查看我们的目标文件

$cat Sample1.txt
genome.fa fastq/Sample1.R1.fastq.gz fastq/Sample1.R2.fastq.gz
$cat Sample2.txt
genome.fa fastq/Sample2.R1.fastq.gz fastq/Sample2.R2.fastq.gz

我们还可以通过Snakemake绘制一个流程图。

1. snakemake --forceall --dag | dot -Tpng> dag.png

（7）运行成功后重新运行

$ snakemake
Building DAG of jobs...
Nothing to be done.

Snakemake在跑完文件后，会给文件加时间戳。因此当你重新跑，已经跑完的规则就不会重跑。

（8）给流程进一步添加规则

下面我们再加一个规则，让流程稍微复杂一点。

1. SAMPLES = ['Sample1', 'Sample2']

2. rule all:

3. input:

4. 'test.txt'

6. rule quantify_genes:

7. input:

8. genome = 'genome.fa',

9. r1 = 'fastq/{sample}.R1.fastq.gz',

10. r2 = 'fastq/{sample}.R2.fastq.gz'

11. output:

12. '{sample}.txt'

13. shell:

14. 'echo {input.genome} {input.r1} {input.r2} > {output}'

16. rule collate_outputs:

17. input:

18. expand('{sample}.txt', sample=SAMPLES)

19. output:

20. 'test.txt'

21. run:

22. with open(output[0], 'w') as out:

23. for i in input:

24. sample = i.split('.')[0]

25. for line in open(i):

26. out.write(sample + ' '+ line)

代码中， rule all的目标文件变更为 test.txt。

rule collate_outputs的输出文件是目标文件，而输入文件是 rule quantify_genes的输出文件，因此Snakemake流程顺序就是 rule all->rule collate_outputs->rule quantify_genes。

rule collate_outputs的功能是将Sample1.txt和 Sample2.txt的每行加上样本名后，写入 test.txt中。

rule collate_outputs是用的python脚本，因此用run而不是shell。run执行python脚本，shell执行Bash脚本。还可以用script来执行外部脚本。比如 script:'scripts/script.py'或 script:'scripts/script.R'。

执行snakemake

$snakemake
Provided cores: 1
Rules claiming more threads will be scaled down.
Job counts:
        count   jobs
1       all
1       collate_outputs
2
rule collate_outputs:
        input: Sample1.txt, Sample2.txt
        output: test.txt
1 of 2 steps (50%) done
localrule all:
        input: test.txt
2 of 2 steps (100%) done

查看目标文件

1. $cat test.txt

2. Sample1 genome.fa ./fastq/Sample1.R1.fastq.gz ./fastq/Sample1.R2.fastq.gz

3. Sample2 genome.fa ./fastq/Sample2.R1.fastq.gz ./fastq/Sample2.R2.fastq.gz

从流程图可以直观的看出运行命令的先后顺序。

snakemake --forceall --dag | dot -Tpng> dag.png

（9）集群投递

前面我们都是直接运行的，如果想通过集群投递任务运行，可以加上cluster和jobs参数。

1. snakemake --jobs 8 --cluster 'qsub -cwd -q normal -l vf=5g,p=1'

这样可以在normal节点，最多一次提交8个任务，最大内存5g，每个任务使用1个cpu。

如果你的任务依赖于其他任务的结果文件，Snakemake会等待直到这个任务的所有依赖文件都生成后再提交。

3.总结

以上只是一个简单的示例。在我们的日常分析中，为了流程的重复使用和方便后续修改维护，会把参数设置、依赖软件的路径、样本信息写在单独的配置文件里。实现不同功能的规则，也会分开放在不同的脚本（smk后缀）里，然后通过include语句将其包含在主Snakefile中。

├── config.yaml
├── softwares.yaml
├── samples.json
├── scripts
│   ├── script1.py
│   └── script2.R
├── rules
│   ├── 01qc.smk
│   ├── 02mapping.smk
│   └── 03plot.smk
└── Snakefile

如果大家对Snakemake感兴趣，也可以参考以下资源做进一步了解。