|
像经常处理生物信息大数据的同学，肯定有自己一套的分析脚本，也就是常说的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}'
上面一共有 all和 quantify_genes两个rule，可以看到前面 all这个rule其实里面没什么具体的执行过程，只是提供一个输入文件，而这个输入文件是后面 quantify_genes的输出文件，这样通过输入和输出文件，把两个规则的前后关系理清了，Snakemake也知道应该先进行 quantify_genes再进行 all。根据上面的规则，我们可以尝试从最终结果逆推回去： 首先在snakefile最前面写了一个SAMPLES的列表，里面保存了2个样品名 all中需要输入文件，这个文件列表由表达式 expand('{sample}.txt', sample=SAMPLES)生成， expand的用法可以参考下面，这里最终生成的结果是 [Sample1.txt, Sample2.txt]两个文件
如果 all中两个输入文件不存在的话，Snakemake还会从其他rule里面进行查找，看看有没有其他输出文件匹配这两个文件。刚好在 quantify_genes中的输出文件就是这两个文件，Snakemake目前把考虑 quantify_genes中的具体内容 在 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)
上面的有几个需要注意的点： 如果需要重新运行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。 参考 Snakemake wiki (http://snakemake./en/stable/index.html) Writing a RNA-Seq workflow with snakemakehttp://pedagogix-tagc./courses/ABD/practical/snakemake/snake_intro.html Build bioinformatics pipelines with Snakemake http:///notes/snakemake-tutorial/
|
