A review of bioinformatics pipeline framework 文献总结了不同流程的优缺点,可以看下。
如果实验室既不是纯粹的生物学试验(不需要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流程的使用做了一些笔记,如下:
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
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以{sample}.fa形式进行文件名通配,用{wildcards.sample}获取sample的实际文件名
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tpdf > dag.pdf
输出结构图
"samtools mpileup -g -f {input.fa} {input.bamA} {input.bamB} | " | bcftools call -mv - > {output} “
["sorted_reads/{}.bam".format(sample) for sample in ["A","B"]]
expand("sorted_reads/{sample}.bam", sample=SAMPLES)
bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
from snakemake.utils import report
with open(input[0]) as vcf:
report””” 内容 """, output[0], T1=input[0])
还用到了snakemake的一个函数,report,可以对markdown语法进行渲染生成网页。
snakemake使用threads来定义当前规则所用的进程数,
samples:
A: data/samples/A.fastq
B: data/samples/B.fastq
configfile: "config.yaml"
虽然sample是一个字典,但是展开的时候,只会使用他们的key值部分。
Input:
lambda wildcards: config["samples"][wildcards.sample]
也就是说在初始化阶段,我们是无法获知通配符所指代的具体文件名,必须要等到第二阶段,才会有wildcards变量出现。也就是说之前的出错的原因都是因为第一个阶段没通过。这个时候就需要输入函数推迟文件名的确定,可以用Python的匿名函数,也可以是普通的函数
有些时候,shell命令不仅仅是由input和output中的文件组成,还需要一些静态的参数设置。如果把这些参数放在input里,则会因为找不到文件而出错,所以需要专门的params用来设置这些参数。
由于高通量测序的数据量通常很大,因此很多无用的中间文件会占据大量的磁盘空间。而特异在执行结束后写一个shell命令清除不但写起来麻烦,而且也不好管理。Snakemake使用temp()来将一些文件标记成临时文件,在执行结束后自动删除。
同时由于比对和排序都比较耗时,得到的结果要是不小心被误删就会浪费大量计算时间,最后的方法就是用protected()保护起来
上面的分析流程都是基于当前环境下已经安装好要调用的软件,如果你希望在新的环境中也能快速部署你的分析流程,那么你需要用到snakmake更高级的特性,也就是为每个rule定义专门的运行环境。
我建议你在新建一个snakemake项目时,都先用conda create -n 项目名 python=版本号创建一个全局环境,用于安装一些常用的软件,例如bwa、samtools、seqkit等。然后用如下命令将环境导出成yaml文件
conda env export -n 项目名 -f environment.yaml
那么当你到了一个新的环境,你就可以用下面这个命令重建出你的运行环境
conda env create -f environment.yaml
snakemake有一个参数--use-conda,会解析rule中的conda规则,根据其提供的yaml文件安装特定版本的工具,以基础第一步的序列比对为例,
output:
"mapped_reads/A.bam"
conda:
"envs/map.yaml"
shell:
随后在snakemake执行的目录下创建envs文件夹,增加map.yaml, 内容如下
name: map
channels:
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
- defaults
dependencies:
- bwa=0.7.17
- samtools=1.9
show_channel_urls: true
注意: YAML文件的name行不是必要的,但是建议加上。
那么当你用snakmake --use-conda执行时,他就会在.snakemake/conda下创建专门的conda环境用于处理当前规则。对于当前项目,该conda环境创建之后就会一直用于该规则,除非yaml文件发生改变。
如果你希望在实际运行项目之前先创建好环境,那么可以使用--create-envs-only参数。
由于默认情况下,每个项目运行时只会在当前的.snakemake/conda查找环境或者安装环境,所以在其他目录执行项目时,snakemake又会重新创建conda环境,如果你担心太占地方或者环境太大,安装的时候太废时间,你可以用--conda-prefix指定专门的文件夹。
snakemake -R some_rule
# --forecerun/-R TARGET: 重新执行给定的规则或生成文件。当你修改规则的时候,使用该命令
snakemake --dag | dot -Tsvg > dag.svg
snakemake --dag | dit -Tpdf > dag.pdf