snakemake搭建生信分析流程

像经常处理生物信息大数据的同学，肯定有自己一套的分析脚本，也就是常说的pipeline了。在我刚开始的时候，也是按照大伙的分析流程，一步一个脚本，后来熟悉了一些，会把经常用到的脚本放在一起，再用变量名替换需要修改的内容，这样每次不用修改具体的文件或者软件参数，直接在脚本里面提供相应的变量名就可以了。如果举个例子的话，可以看看技能树里面的一个帖子史上最快的转录组定量流程，这里面就是把常用的脚本直接写在shell里面，而且很机智的使用了配置文件，脚本中也是读取配置文件中的配置信息，用于后续分析流程的配置。

snakemake介绍

上面例子中的方法当然可以应付大部分生物信息分析流程了，可是如果你需要分析群体数据时，而且这个群体样本数很多的时候，用简单的脚本处理，就比较麻烦了，而且有的时候，分析流程中的某一部还会莫名的出错而导致流程的中断，这个时候进行排查，也是一件挺麻烦的事情。

这里为大家提供snakemake工作流管理软件。Snakemake是用Python3写的一个工具，模仿makefile的格式（makefile就是我们安装软件 make make install时，使用的各种编译规则，不明白也没关系），处理流程中各步骤的依赖关系，这样就知道哪一步分析先运行，哪一步后运行，各个分析步骤就形成了一个有向无环图（DAG，就是GO分析里面的DOG喽），引用文档里面的一个图：

Snakemake里面最重要的有概念有3个，虽然很重要，可是真的很简单啊：

  • input files， 定义了输入文件

  • output files， 定义了输出文件

  • rules，定义了输入文件生成输出文件过程

下面介绍使用方法的时候，大家可要一定记住这3个概念。

snakemake使用

安装

首先是安装snakemake，snakemake已经整理成Python包，可以直接使用 pip进行安装，不过需要的Python3的环境，可以愉快地利用 conda进行安装：

  1. conda create --name Py3 python==3.5

  2. source activate Py3

  3. # 这里的pip是指pip3

  4. pip install snakemake

试试 snakemake -h看看安装成功没有？

使用说明

这里参照文档中的一个例子进行说明。首先看看一个Snakefile的样子：

  1. rule bwa_map:

  2.    input:

  3.        "data/genome.fa",

  4.        "data/samples/A.fastq"

  5.    output:

  6.        "mapped_reads/A.bam"

  7.    shell:

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

Skefile中只定义一个叫bwa_map的rule，rule里面分别写明了输入和输出文件，还有运行脚本的命令。在shell中用了{input}和{output}引用了输入和输出文件，对于下面打印出来的脚本，可以看到{input}里面有顺序的。

前面假装已经帮参考序列genome.fa建好了索引，这里直接使用 bwa就可以了。相关的数据文件已经保存到我的github上，大家直接下载使用就可以。

  1. cd snakemake-tutorial

  2. snakemake -np

这里不执行程序，而是把相关的脚本打印出来，同时还有任务统计信息：

  1. rule bwa_map:  

  2.    input: data/genome.fa, data/samples/A.fastq  

  3.    output: mapped_reads/A.bam  

  4. bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam  

  5. Job counts:  

  6.    count   jobs  

  7.    1   bwa_map  

  8.    1  

上面列出了rule中的input和output，最终执行的脚本以及任务的统计，这里只是简单的一个任务，只是简单说明了使用过程。

前面说到如果有很多文件的话，用Snakemake处理也很方便，具体可以看下面这个例子，相关文件保存于 multi_inputs中，文件 snakefile：

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

  2. rule all:

  3.    input:

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

  5. rule quantify_genes:

  6.    input:

  7.        genome = 'genome.fa',

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

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

  10.    output:

  11.        '{sample}.txt'

  12.    shell:

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

上面一共有 all和 quantify_genes两个rule，可以看到前面 all这个rule其实里面没什么具体的执行过程，只是提供一个输入文件，而这个输入文件是后面 quantify_genes的输出文件，这样通过输入和输出文件，把两个规则的前后关系理清了，Snakemake也知道应该先进行 quantify_genes再进行 all。根据上面的规则，我们可以尝试从最终结果逆推回去：

  1. 首先在snakefile最前面写了一个SAMPLES的列表，里面保存了2个样品名

  2. all中需要输入文件，这个文件列表由表达式 expand('{sample}.txt', sample=SAMPLES)生成， expand的用法可以参考下面，这里最终生成的结果是 [Sample1.txt, Sample2.txt]两个文件

  3. 如果 all中两个输入文件不存在的话，Snakemake还会从其他rule里面进行查找，看看有没有其他输出文件匹配这两个文件。刚好在 quantify_genes中的输出文件就是这两个文件，Snakemake目前把考虑 quantify_genes中的具体内容

  4. 在 quantify_genes中要想生成最终两个文件，还需要考虑自己需要哪些输入文件，刚好这里也定义了输入文件信息，这里定义输入文件的格式与上面的还不一样，采用了键值对的形式进行设置，而且在后面使用的时候，也是采用了类似Python属性调用的形式（或者使用索引也是可以的 input[0]，但可读性不好），里面的{sample}变量应该和前面 all中的{sample}来源不一样，Snakemake会根据 fastq目录中存在的文件进行匹配，匹配的最终也会得到 ['Sample1', 'Sample2']，最后的shell执行过程也只是很简单的 echo语句

上面的例子只是简单说明了不同rule之间是怎么进行联系，确定先后的运行顺序，其实就是很简单的一个脚本，偏偏写这么复杂，其实Snakefile应该有更简单的方法，读者学习后，可以想想要怎么写会更加方便。

例子中用到的 expand函数，可以方便地利用一些简单的列表和基本模板，得到文件列表，使用方法可以简单看下面这两个例子：

  1. expand("sorted_reads/{sample}.bam", sample=SAMPLES)

  2. # ["sorted_reads/A.bam", "sorted_reads/B.bam"]

  3. # 可以提供多个参数列表

  4. expand("sorted_reads/{sample}.{replicate}.bam", sample=SAMPLES, replicate=[0, 1])

  5. # ["sorted_reads/A.0.bam", "sorted_reads/A.1.bam", "sorted_reads/B.0.bam", "sorted_reads/B.1.bam"]

glob_wildcards的使用

在 fastq目录中一共有两组paired-end测序数据：

  1. fastq

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

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

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

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

类似这种文件名存在一定规律的文件，Snakemake提供了类似bash里的通配符一样的功能，把匹配上的信息提取出来。Snakemake这里提供了 glob_wildcards这个函数来实现。

glob_wildcards可以接受一个或多个通配符表达式，类似 {pattern}，最后返回一个由list组成的tuple，所以如果只有一个变量的话，别忘了添加逗号（具体可以参考下面的例子）。还有一点就是 glob_wildcards不支持bash中用的通配符（ ?, *)。

  1. # Example 1

  2. (SAMPLE_LIST,) = glob_wildcards("fastq/{sample}.fastq")

需要注意下SAMPLE_LIST后面的逗号。

  1. # Example 2

  2. (genome_builds, samples, ...) = glob_wildcards("{pattern1}/{pattern2}.whatever")

例子2说明，可以接收多个匹配模板。如果存在下面文件：

  1. grch37/A.bam

  2. grch37/B.bam

  3. grch38/A.bam

  4. grch38/B.bam

那 (targets, samples) = glob_wildcards("{target}/{sample}.bam")返回的结果将是：

  1. targets = ['grch37', 'grch38']

  2. samples = ['A', 'B']

在上面的规则中，我们直接定义了 SAMPLES = ['Sample1', 'Sample2']两个变量，现在尝试用通配符来自动匹配，文件 multi_inputs/snakefile2：

  1. from os.path import join

  2. # Globals ---------------------------------------------------------------------

  3. # Full path to a FASTA file.

  4. GENOME = 'genome.fa'

  5. # A Snakemake regular expression matching the forward mate FASTQ files.

  6. SAMPLES, = glob_wildcards(join('fastq', '{sample,Samp[^/]+}.R1.fastq.gz'))

  7. # Patterns for the 1st mate and the 2nd mate using the 'sample' wildcard.

  8. PATTERN_R1 = '{sample}.R1.fastq.gz'

  9. PATTERN_R2 = '{sample}.R2.fastq.gz'

  10. # Rules -----------------------------------------------------------------------

  11. rule all:

  12.    input:

  13.        'test2.txt'

  14. rule quantify_genes:

  15.    input:

  16.        genome = GENOME,

  17.        r1 = join('fastq', PATTERN_R1),

  18.        r2 = join('fastq', PATTERN_R2)

  19.    output:

  20.        '{sample}.txt'

  21.    shell:

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

  23. rule collate_outputs:

  24.    input:

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

  26.    output:

  27.        'test2.txt'

  28.    run:

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

  30.            for i in input:

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

  32.                for line in open(i):

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

上面的有几个需要注意的点：

  • 将匹配的结果保存在 SAMPLES中，可是后面的使用上有 {sample}，还有 SAMPLES，可以看到 SAMPLES是不能直接作为输入输出的变量，需要借助 expand

  • collate_outputs这块规则的后面运行程序是直接运行Python代码的，还可以运行R代码（参考链接2）

如果需要重新运行Snakemake的话，会先判断相关结果是否存在，而决定是否重新运行。也就是说，如果前面有部分流程中断，那重新运行，就会把一些中断流程中步骤重新运行。当然如果某个文件只生成了一半，也会当成运行成功的。这个时候可以用 --forcerules或 -R让Snakemake强制重新运行

  1. snakemake -R somerule

如果想重新运行上面的流程，可以使用：

  1. snakemake -R all --snakefile snakefile2

如果想重新某一步，可以指定相应的rule名称。

集群任务投递

如果服务器不是计算集群的话，可以直接使用 snakemake --snakefile Snakefile来进行投递，可以配合 nohup转到后台执行。 
对于使用 qsub调度系统的（e.g. the Sun Grid Engine），可以使用下面的命令进行提交任务：

  1. snakemake --cluster "qsub -l vf=4G,p=1"

最后

上面介绍的只是简单介绍了Snakemake的使用，其实Snakemake还有许多有意思的功能和运用，比如可以用图形界面查看运行状态，不同流程之间的引用和组合。如果有兴趣的话，可以多参考官方文档和其他用Snakemake写的流程，如果遇到问题也可以在StackOverflow上提问。

相关代码及文件已经保存到我的Github中，可以自由下载使用，如果有问题，也请添加相关的Issue。

参考

  1. Snakemake wiki  (http://snakemake.readthedocs.io/en/stable/index.html)

  2. Writing a RNA-Seq workflow with snakemakehttp://pedagogix-tagc.univ-mrs.fr/courses/ABD/practical/snakemake/snake_intro.html

  3. Build bioinformatics pipelines with Snakemake http://slowkow.com/notes/snakemake-tutorial/

赞 (0)

相关推荐