基于python的snakemake框架搭建生物信息学分析流程笔记总结

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