A review of bioinformatics pipeline framework 的作者对已有的工具进行很好的分类
作者的看法:
implicit,也就是Make rule语法更适合用于整合不同执行工具
基于配置的流程更加稳定,也比较适合用于集群分配任务。
最后作者建议是:
如果实验室既不是纯粹的生物学试验(不需要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。
如下步骤会下载所需数据,并安装所需要的软件,并且启动工作环境。
wget https://bitbucket.org/snakemake/snakemake-tutorial/get/v3.11.0.tar.bz2
tar -xf v3.11.0.tar.bz2 --strip 1
cd snakemake-snakemake-tutorial-623791d7ec6d
conda env create --name snakemake-tutorial --file environment.yaml
source activate snakemake-tutorial
# 退出当前环境
source deactivate
当前环境下的所有文件
├── data
│ ├── genome.fa
│ ├── genome.fa.amb
│ ├── genome.fa.ann
│ ├── genome.fa.bwt
│ ├── genome.fa.fai
│ ├── genome.fa.pac
│ ├── genome.fa.sa
│ └── samples
│ ├── A.fastq
│ ├── B.fastq
│ └── C.fastq
├── environment.yaml
└── 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的规则就是下面的写法。
# 用你擅长的文本编辑器
vim 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}
'''
解释一下:这几行定义了一个规则(rule),在这个规则下,输入(input)有两个,而输出(output)只有一个,在 shell
中运行命令,只不过里面的文件都用 {}
形式替代。伪执行一下: snakemake -np mapped_reads/A.bam
检查一下是否会出错,真实运行情况如下(不带规则,默认执行第一个规则):
第二步:使用通配推广序列比对规则
如果仅仅是上面这样子处理一个文件,还无法体现 snakemake
的用途,毕竟还不如手动敲代码来的方便。 snakemake
的一个有点在于它能够使用文件名通配的方式对一类文件进行处理。将上面的 A
改成 {sample}
,就可以将符合 *.fastq
的文件处理成 *.bam
.
rule bwa_map:
input:
'data/genome.fa',
'data/samples/{sample}.fastq'
output:
'mapped_reads/{sample}.bam'
shell:
'''
bwa mem {input} | samtools view -Sb - > {output}
'''
那么,用 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。
第三步:比对后排序
比对后的文件还需要进一步的排序,才能用于后续分析,那么规则该如何写呢?
rule samtools_sort:
input:
'mapped_reads/{sample}.bam'
output:
'sorted_reads/{sample}.bam'
shell:
'samtools sort -T sorted_reads/{wildcards.sample}'
' -O bam {input} > {output}'
以之前的输出作为输出文件名,输出到另一个文件夹中。和之前的规则基本相同,只不过这里用到了 wildcards.sample
来获取通配名用作 -T
的临时文件的前缀 sample
实际名字。
运行 snakemake -np sorted_reads/B.bam
,你就会发现他就会非常智能的先比对再排序。这是因为 snakemake
会自动解决依赖关系,并且按照依赖的前后顺序进行执行。
第四步: 建立索引和对任务可视化
这里我们再写一个规则,对之前的排序后的BAM文件建立索引。
rule samtools_index:
input:
'sorted_reads/{sample}.bam'
output:
'sorted_reads/{sample}.bam.bai'
shell:
'samtools index {input}'
目前已经写了三个规则,那么这些规则的执行和依赖关系如何呢? snakemake
提供了 --dag
选项用于 dot
命令进行可视化
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
第五步:基因组变异识别
基因组变异识别需要整合之前所有的BAM文件,你可能会打算这样写
rule bcftools_call:
input:
fa='data/genome.fa',
bamA='sorted_reads/A.bam'
bamB='sorted_reads/B.bam'
baiA='sorted_reads/A.bam.bai'
baiB='sorted_reads/B.bam.bai'
output:
'calls/all.vcf'
shell:
'samtools mpileup -g -f {input.fa} {input.bamA} {input.bamB} | '
'bcftools call -mv - > {output}'
这样写的却没有问题,但是以后每多一个样本就需要多写一个输入,太麻烦了。这里就体现出Snakemake和Python所带来的特性了,我们可以用列表推导式的方法搞定。
['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)
那么最后的规则就是
SAMPLES=['A','B']
rule bcftools_call:
input:
fa='data/genome.fa',
bam=expand('sorted_reads/{sample}.bam', sample=SAMPLES),
bai=expand('sorted_reads/{sample}.bam.bai', sample=SAMPLES)
output:
'calls/all.vcf'
shell:
'samtools mpileup -g -f {input.fa} {input.bam} | '
'bcftools call -mv - > {output}'
小练习: 请用snakemake生成当前的DAG图。
第六步:编写报告
上面都是在规则里执行shell脚本,snakemake的一个优点就是可以在规则里面写Python脚本,只需要把 shell
改成 run
,此外还不需要用到引号。
rule report:
input:
'calls/all.vcf'
output:
'report.html'
run:
from snakemake.utils import report
with open(input[0]) as vcf:
n_calls = sum(1 for l in vcf if not l.startswith('#'))
report('''
An example variant calling workflow
===================================
Reads were mapped to the Yeast
reference genome and variants were called jointly with
SAMtools/BCFtools.
This resulted in {n_calls} variants (see Table T1_).
''', output[0], T1=input[0])
这里还用到了 snakemake
的一个函数,report,可以对markdown语法进行渲染生成网页。
第七步:增加目标规则
之前运行snakemake都是用的 snakemake 目标文件名
, 除了目标文件名外,snakemake还支持规则名作为目标。通常我们按照习惯定义一个 all
规则,来生成结果文件。
rule all:
input:
'report.html
Snakemake基于规则执行命令,规则一般由 input, output,shell
三部分组成。
Snakemake可以自动确定不同规则的输入输出的依赖关系,根据时间戳来判断文件是否需要重新生成
Snakemake 以{sample}.fa
形式进行文件名通配,用 {wildcards.sample}
获取sample的实际文件名
Snakemake用 expand()
生成多个文件名,本质是Python的列表推导式
Snakemake可以在规则外直接写Python代码,在规则内的 run
里也可以写Python代码。
Snakefile的第一个规则通常是 rule all
,因为默snakemake默认执行第一条规则
最后附上代码如下, 希望你喜欢:
SAMPLES = ['A', 'B']
rule all:
input:
'report.html'
rule bwa_map:
input:
'data/genome.fa',
'data/samples/{sample}.fastq'
output:
'mapped_reads/{sample}.bam'
shell:
'bwa mem {input} | samtools view -Sb - > {output}'
rule samtools_sort:
input:
'mapped_reads/{sample}.bam'
output:
'sorted_reads/{sample}.bam'
shell:
'samtools sort -T sorted_reads/{wildcards.sample} '
'-O bam {input} > {output}'
rule samtools_index:
input:
'sorted_reads/{sample}.bam'
output:
'sorted_reads/{sample}.bam.bai'
shell:
'samtools index {input}'
rule bcftools_call:
input:
fa='data/genome.fa',
bam=expand('sorted_reads/{sample}.bam', sample=SAMPLES),
bai=expand('sorted_reads/{sample}.bam.bai', sample=SAMPLES)
output:
'calls/all.vcf'
shell:
'samtools mpileup -g -f {input.fa} {input.bam} | '
'bcftools call -mv - > {output}'
rule report:
input:
'calls/all.vcf'
output:
'report.html'
run:
from snakemake.utils import report
with open(input[0]) as vcf:
n_calls = sum(1 for l in vcf if not l.startswith('#'))
report('''
An example variant calling workflow
===================================
Reads were mapped to the Yeast
reference genome and variants were called jointly with
SAMtools/BCFtools.
This resulted in {n_calls} variants (see Table T1_).
''', output[0], T1=input[0])
参考资料
http://snakemake.readthedocs.io/en/latest/tutorial/basics.html
http://sequana.readthedocs.io/en/master/
https://adoubi.life/posts/2017-10-17/
联系客服