像经常处理生物信息大数据的同学,肯定有自己一套的分析脚本,也就是常说的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进行安装:

conda create --name Py3 python==3.5 
source activate Py3
# 这里的pip是指pip3
pip install snakemake

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

使用说明

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

rule bwa_map:
    input:
        "data/genome.fa",
        "data/samples/A.fastq"
    output:
        "mapped_reads/A.bam"
    shell:
        "bwa mem {input} | samtools view -Sb - > {output}"

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

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

cd snakemake-tutorial
snakemake -np

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

rule bwa_map:
input: data/genome.fa, data/samples/A.fastq
output: mapped_reads/A.bam

bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam
Job counts:
count jobs
1 bwa_map
1

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


前面说到如果有很多文件的话,用Snakemake处理也很方便,具体可以看下面这个例子,相关文件保存于multi_inputs中,文件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}'

上面一共有allquantify_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函数,可以方便地利用一些简单的列表和基本模板,得到文件列表,使用方法可以简单看下面这两个例子:

expand("sorted_reads/{sample}.bam", sample=SAMPLES)
# ["sorted_reads/A.bam", "sorted_reads/B.bam"]

# 可以提供多个参数列表
expand("sorted_reads/{sample}.{replicate}.bam", sample=SAMPLES, replicate=[0, 1])
# ["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测序数据:

fastq
├── Sample1.R1.fastq.gz
├── Sample1.R2.fastq.gz
├── Sample2.R1.fastq.gz
└── Sample2.R2.fastq.gz

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

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

# Example 1
(SAMPLE_LIST,) = glob_wildcards("fastq/{sample}.fastq")

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

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

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

grch37/A.bam
grch37/B.bam
grch38/A.bam
grch38/B.bam

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

targets = ['grch37', 'grch38']
samples = ['A', 'B']

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

from os.path import join

# Globals ---------------------------------------------------------------------

# Full path to a FASTA file.
GENOME = 'genome.fa'

# A Snakemake regular expression matching the forward mate FASTQ files.
SAMPLES, = glob_wildcards(join('fastq', '{sample,Samp[^/]+}.R1.fastq.gz'))

# Patterns for the 1st mate and the 2nd mate using the 'sample' wildcard.
PATTERN_R1 = '{sample}.R1.fastq.gz'
PATTERN_R2 = '{sample}.R2.fastq.gz'

# Rules -----------------------------------------------------------------------

rule all:
    input:
        'test2.txt'

rule quantify_genes:
    input:
        genome = GENOME,
        r1 = join('fastq', PATTERN_R1),
        r2 = join('fastq', PATTERN_R2)
    output:
        '{sample}.txt'
    shell:
        'echo {input.genome} {input.r1} {input.r2} > {output}'

rule collate_outputs:
    input:
        expand('{sample}.txt', sample=SAMPLES)
    output:
        'test2.txt'
    run:
        with open(output[0], 'w') as out:
            for i in input:
                sample = i.split('.')[0]
                for line in open(i):
                    out.write(sample + ' ' + line)

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

  • 将匹配的结果保存在SAMPLES中,可是后面的使用上有{sample},还有SAMPLES,可以看到SAMPLES是不能直接作为输入输出的变量,需要借助expand
  • collate_outputs这块规则的后面运行程序是直接运行Python代码的,还可以运行R代码(参考链接2)

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

snakemake -R somerule

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

snakemake -R all --snakefile snakefile2

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

集群任务投递

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

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

对于其他类型,可以参考这里

最后

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

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

参考:

  1. Snakemake wiki
  2. Writing a RNA-Seq workflow with snakemake
  3. Build bioinformatics pipelines with Snakemake