这里有一个流程，希望你会喜欢

A review of bioinformatics pipeline framework 的作者对已有的工具进行很好的分类

工具推荐

作者的看法：

1. implicit，也就是Make rule语法更适合用于整合不同执行工具

2. 基于配置的流程更加稳定，也比较适合用于集群分配任务。

最后作者建议是：

• 如果实验室既不是纯粹的生物学试验（不需要workbench这种UI界面），也不需要高性能基于类的流程设计， 不太好选， 主要原则是投入和产出比

• 如果实验室进行的是重复性的研究，那么就需要对数据和软件进行版本控制， 建议是 configuration-based pipelines

• 如果实验室做的是探索性的概念证明类工作（exploratory proofs-of-concept)，那么需要的是 DSL-based pipeline。

• 如果实验室用不到高性能计算机(HPC)，只能用云服务器，就是server-based frameworks.

目前已有的流程可以在awesome-pipeline 进行查找。

就目前来看，pipeline frameworks & library 这部分的框架中 nextflow 是点赞数最多的生物学相关框架。只可惜nextflow在运行时需要创建fifo，而在NTFS文件系统上无法创建，所以我选择 snakemake , 一个基于Python写的DSL流程框架。

环境准备

为了能够顺利完成这部分的教程，请准备一个Linux环境，如果使用Windows，则按照biostarhandbook(一)分析环境和数据可重复部署一个虚拟机，并安装miniconda3。

如下步骤会下载所需数据，并安装所需要的软件，并且启动工作环境。

1. wget https://bitbucket.org/snakemake/snakemake-tutorial/get/v3.11.0.tar.bz2

2. tar -xf v3.11.0.tar.bz2 --strip 1

3. cd snakemake-snakemake-tutorial-623791d7ec6d

4. conda env create --name snakemake-tutorial --file environment.yaml

5. source activate snakemake-tutorial

6. # 退出当前环境

7. source deactivate

当前环境下的所有文件

1. ├── data

2. │   ├── genome.fa

3. │   ├── genome.fa.amb

4. │   ├── genome.fa.ann

5. │   ├── genome.fa.bwt

6. │   ├── genome.fa.fai

7. │   ├── genome.fa.pac

8. │   ├── genome.fa.sa

9. │   └── samples

10. │   ├── A.fastq

11. │   ├── B.fastq

12. │   └── C.fastq

13. ├── environment.yaml

14. └── README.md

基础：一个流程实例

如果你编译过软件，那你应该见过和用过 make, 但是你估计也没有仔细想过make是干嘛用的。Make是最常用的软件构建工具，诞生于1977年，主要用于C语言的项目，是为了处理编译时存在各种依赖关系，尤其是部分文件更新后，Make能够重新生成需要更新的文件以及其对应的文件。

Snakemake和Make功能一致，只不过用Python实现，增加了许多Python的特性，并且和Python一样非常容易阅读。下面将使用Snakemake写一个变异检测流程。

第一步：序列比对

Snakemake非常简单，就是写各种rule来完成不同的任务。我们第一条rule就是将序列比对到参考基因组上。如果在命令行下就是 bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam。 但是按照Snakemake的规则就是下面的写法。

1. # 用你擅长的文本编辑器

2. vim Snakefile

3. # 编辑如下内容

4. rule bwa_map:

5.  input:

6.  'data/genome.fa',

7.  'data/samples/A.fastq'

8.  output:

9.  'mapped_reads/A.bam'

10.  shell:

11.  '''

12.  bwa mem {input} | samtools view -Sb - > {output}

13.  '''

解释一下：这几行定义了一个规则(rule)，在这个规则下，输入(input)有两个，而输出(output)只有一个，在 shell中运行命令，只不过里面的文件都用 {}形式替代。伪执行一下: snakemake -np mapped_reads/A.bam检查一下是否会出错，真实运行情况如下（不带规则，默认执行第一个规则）:

run snakemake

第二步：使用通配推广序列比对规则

如果仅仅是上面这样子处理一个文件，还无法体现 snakemake的用途，毕竟还不如手动敲代码来的方便。 snakemake的一个有点在于它能够使用文件名通配的方式对一类文件进行处理。将上面的 A改成 {sample},就可以将符合 *.fastq的文件处理成 *.bam.

1. rule bwa_map:

2.  input:

3.  'data/genome.fa',

4.  'data/samples/{sample}.fastq'

5.  output:

6.  'mapped_reads/{sample}.bam'

7.  shell:

8.  '''

9.  bwa mem {input} | samtools view -Sb - > {output}

10.  '''

那么，用 snakemake -np mapped_reads/{A,B,C}.bam，就会发现，他非常机智的就比对了 B.fastq和 C.fastq，而不会再比对一遍A.fastq, 也不需要你写一堆的判断语句去手动处理。

规则统配

当然，如果你用 touch data/samples/A.fastq改变A.fastq的时间戳，他就会认位A.fastq文件发生了改变，那么重复之前的命令就会比对A.fastq。

第三步：比对后排序

比对后的文件还需要进一步的排序，才能用于后续分析，那么规则该如何写呢？

1. rule samtools_sort:

2.  input:

3.  'mapped_reads/{sample}.bam'

4.  output:

5.  'sorted_reads/{sample}.bam'

6.  shell:

7.  'samtools sort -T sorted_reads/{wildcards.sample}'

8.  ' -O bam {input} > {output}'

以之前的输出作为输出文件名，输出到另一个文件夹中。和之前的规则基本相同，只不过这里用到了 wildcards.sample来获取通配名用作 -T的临时文件的前缀 sample实际名字。

运行 snakemake -np sorted_reads/B.bam，你就会发现他就会非常智能的先比对再排序。这是因为 snakemake会自动解决依赖关系，并且按照依赖的前后顺序进行执行。

第四步： 建立索引和对任务可视化

这里我们再写一个规则，对之前的排序后的BAM文件建立索引。

1. rule samtools_index:

2.  input:

3.  'sorted_reads/{sample}.bam'

4.  output:

5.  'sorted_reads/{sample}.bam.bai'

6.  shell:

7.  'samtools index {input}'

目前已经写了三个规则，那么这些规则的执行和依赖关系如何呢？ snakemake提供了 --dag选项用于 dot命令进行可视化

1. snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg

运行流程

第五步：基因组变异识别

基因组变异识别需要整合之前所有的BAM文件，你可能会打算这样写

1. rule bcftools_call:

2.  input:

3.  fa='data/genome.fa',

4.  bamA='sorted_reads/A.bam'

5.  bamB='sorted_reads/B.bam'

6.  baiA='sorted_reads/A.bam.bai'

7.  baiB='sorted_reads/B.bam.bai'

8.  output:

9.  'calls/all.vcf'

10.  shell:

11.  'samtools mpileup -g -f {input.fa} {input.bamA} {input.bamB} | '

12.  'bcftools call -mv - > {output}'

这样写的却没有问题，但是以后每多一个样本就需要多写一个输入，太麻烦了。这里就体现出Snakemake和Python所带来的特性了，我们可以用列表推导式的方法搞定。

1. ['sorted_reads/{}.bam'.format(sample) for sample in ['A','B']]

进一步，可以在规则外定义 SAMPLES=['A','B']，则规则内的输入可以写成 bam=['sorted_reads/{}.bam'.format(sample) for sample in SAMPLES]. 由于列表推导式比较常用，但是写起来有点麻烦，snakemake定义了 expand进行简化, 上面可以继续改写成 expand('sorted_reads/{sample}.bam', sample=SAMPLES)

那么最后的规则就是

1. SAMPLES=['A','B']

2. rule bcftools_call:

3.  input:

4.  fa='data/genome.fa',

5.  bam=expand('sorted_reads/{sample}.bam', sample=SAMPLES),

6.  bai=expand('sorted_reads/{sample}.bam.bai', sample=SAMPLES)

7.  output:

8.  'calls/all.vcf'

9.  shell:

10.  'samtools mpileup -g -f {input.fa} {input.bam} | '

11.  'bcftools call -mv - > {output}'

小练习： 请用snakemake生成当前的DAG图。

第六步：编写报告

上面都是在规则里执行shell脚本，snakemake的一个优点就是可以在规则里面写Python脚本，只需要把 shell改成 run，此外还不需要用到引号。

1. rule report:

2.  input:

3.  'calls/all.vcf'

4.  output:

5.  'report.html'

6.  run:

7.  from snakemake.utils import report

8.  with open(input[0]) as vcf:

9.  n_calls = sum(1 for l in vcf if not l.startswith('#'))

11.  report('''

12.  An example variant calling workflow

13.  ===================================

15.  Reads were mapped to the Yeast

16.  reference genome and variants were called jointly with

17.  SAMtools/BCFtools.

19.  This resulted in {n_calls} variants (see Table T1_).

20.  ''', output[0], T1=input[0])

这里还用到了 snakemake的一个函数，report，可以对markdown语法进行渲染生成网页。

第七步：增加目标规则

之前运行snakemake都是用的 snakemake 目标文件名, 除了目标文件名外，snakemake还支持规则名作为目标。通常我们按照习惯定义一个 all规则，来生成结果文件。

1. rule all:

2.  input:

3.  'report.html

总结

总结以下学习过程，知识点如下：

• Snakemake基于规则执行命令，规则一般由 input, output,shell三部分组成。

• Snakemake可以自动确定不同规则的输入输出的依赖关系，根据时间戳来判断文件是否需要重新生成

• Snakemake 以{sample}.fa形式进行文件名通配，用 {wildcards.sample}获取sample的实际文件名

• Snakemake用 expand()生成多个文件名，本质是Python的列表推导式

• Snakemake可以在规则外直接写Python代码，在规则内的 run里也可以写Python代码。

• Snakefile的第一个规则通常是 rule all，因为默snakemake默认执行第一条规则

最后附上代码如下, 希望你喜欢：

1. SAMPLES = ['A', 'B']

4. rule all:

5.  input:

6.  'report.html'

9. rule bwa_map:

10.  input:

11.  'data/genome.fa',

12.  'data/samples/{sample}.fastq'

13.  output:

14.  'mapped_reads/{sample}.bam'

15.  shell:

16.  'bwa mem {input} | samtools view -Sb - > {output}'

19. rule samtools_sort:

20.  input:

21.  'mapped_reads/{sample}.bam'

22.  output:

23.  'sorted_reads/{sample}.bam'

24.  shell:

25.  'samtools sort -T sorted_reads/{wildcards.sample} '

26.  '-O bam {input} > {output}'

29. rule samtools_index:

30.  input:

31.  'sorted_reads/{sample}.bam'

32.  output:

33.  'sorted_reads/{sample}.bam.bai'

34.  shell:

35.  'samtools index {input}'

38. rule bcftools_call:

39.  input:

40.  fa='data/genome.fa',

41.  bam=expand('sorted_reads/{sample}.bam', sample=SAMPLES),

42.  bai=expand('sorted_reads/{sample}.bam.bai', sample=SAMPLES)

43.  output:

44.  'calls/all.vcf'

45.  shell:

46.  'samtools mpileup -g -f {input.fa} {input.bam} | '

47.  'bcftools call -mv - > {output}'

50. rule report:

51.  input:

52.  'calls/all.vcf'

53.  output:

54.  'report.html'

55.  run:

56.  from snakemake.utils import report

57.  with open(input[0]) as vcf:

58.  n_calls = sum(1 for l in vcf if not l.startswith('#'))

60.  report('''

61.  An example variant calling workflow

62.  ===================================

64.  Reads were mapped to the Yeast

65.  reference genome and variants were called jointly with

66.  SAMtools/BCFtools.

68.  This resulted in {n_calls} variants (see Table T1_).

69.  ''', output[0], T1=input[0])

参考资料

• http://snakemake./en/latest/tutorial/basics.html

• http://sequana./en/master/

• https:///posts/2017-10-17/