snakemake 初探

前方的路必定充满荆棘

1引言

什么是 snakemake?

Snakemake 工作流管理系统是一种创建 可重复 和 可扩展数据分析 的工具。工作流程是可读性的，基于 Python 语言描述的。它们可以无缝缩放到服务器，群集，网格和云环境，而无需修改工作流定义。最后，Snakemake 工作流程可能所需的软件，它将自动部署到任何执行环境。

对于比较长的数据分析流程或者经常用的流程，使用 snakemake 将会大大减少分析时间和减少代码管理。

我们跟随 snakemake 官网教程学习一些基础知识和用法。

2安装

安装 mamba ，它是 conda 的一个替代品，也是安装和管理软件的工具：

$ conda install -n base -c conda-forge mamba

snakemake 基于 python3 环境，创建一个 snakemake 环境并安装 snakemake：

$ mamba create -c conda-forge -c bioconda -n snakemake snakemake

激活环境，查看帮助：

$ conda activate snakemake
$ snakemake --help
usage: snakemake [-h] [--snakefile FILE] [--gui [PORT]] [--cores [N]]
                 [--local-cores N] [--resources [NAME=INT [NAME=INT ...]]]
                 [--config [KEY=VALUE [KEY=VALUE ...]]] [--configfile FILE]
                 [--list] [--list-target-rules] [--directory DIR] [--dryrun]
                 [--printshellcmds] [--debug-dag] [--dag]
                 [--force-use-threads] [--rulegraph] [--d3dag] [--summary]
                 [--detailed-summary] [--archive FILE] [--touch]
                 [--keep-going] [--force] [--forceall]
                 [--forcerun [TARGET [TARGET ...]]]
                 [--prioritize TARGET [TARGET ...]]
                 [--until TARGET [TARGET ...]]
                 [--omit-from TARGET [TARGET ...]] [--allow-ambiguity]
                 [--cluster CMD | --cluster-sync CMD | --drmaa [ARGS]]
                 [--drmaa-log-dir DIR] [--cluster-config FILE]
                 [--immediate-submit] [--jobscript SCRIPT] [--jobname NAME]
                 [--reason] [--stats FILE] [--nocolor] [--quiet] [--nolock]
                 [--unlock] [--cleanup-metadata FILE [FILE ...]]
                 [--rerun-incomplete] [--ignore-incomplete]
                 [--list-version-changes] [--list-code-changes]
                 [--list-input-changes] [--list-params-changes]
                 [--latency-wait SECONDS] [--wait-for-files [FILE [FILE ...]]]
                 [--benchmark-repeats N] [--notemp] [--keep-remote]
                 [--keep-target-files] [--keep-shadow]
                 [--allowed-rules ALLOWED_RULES [ALLOWED_RULES ...]]
                 [--max-jobs-per-second MAX_JOBS_PER_SECOND]
                 [--restart-times RESTART_TIMES] [--timestamp]
                 [--greediness GREEDINESS] [--no-hooks] [--print-compilation]
                 [--overwrite-shellcmd OVERWRITE_SHELLCMD] [--verbose]
                 [--debug] [--profile FILE] [--mode {0,1,2}]
                 [--bash-completion] [--use-conda] [--conda-prefix DIR]
                 [--wrapper-prefix WRAPPER_PREFIX]
                 [--default-remote-provider {S3,GS,SFTP,S3Mocked}]
                 [--default-remote-prefix DEFAULT_REMOTE_PREFIX] [--version]
                 [target [target ...]]

3基本参数说明

常用命令：

• --dryrun, -n： 不真正执行。
• --printshellcmds, -p： 输出 shell 命令。
• --reason, -r： 打印执行规则的原因。
• --snakefile FILE, -s： 指定 snakemake 文件。

检查测试：

$ snakemake -n
$ snakemake -np
$ snakemake -nr

运行 pipeline：

$ snakemake # 默认在当前目录下的snakemake文件
$ snakemake -s ./my_test_snakemake.py

生成流程图：

# 生成svg图片
$ snakemake --dag  | dot -Tsvg > dag.svg
# 生成pdf文件
$ snakemake --dag  | dit -Tpdf > dag.pdf
# 通过网页查看状态
$ snakemake --gui 0.0.0.0:2468

snakemake 最主要的是需要我们我们去 书写 snakemake 文件，里面定义了我们的 规则（rule）及分析流程的代码。

4示例

比对单个文件

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

创建一个以 .py 结尾的文件，里面输入以上内容，每个模块以 rule + 模块名称 开始，后面最基本的组成为 输入（input）、输出（output） 和 命令（shell）。此外区分 tab 补齐。

上面定义了一个 bwa 比对的模块，输入文件为 data/genome.fa 和 data/samples/A.fastq ，输出文件为 mapped_reads/A.bam ，输出文件夹如果不存在则会自动创建，运行的 shell 命令则是 bwa mem {input} | samtools view -Sb - > {output} 比对输出为 sam 文件。

值得注意的是，我们需要在当前环境已经安装好 bwa 和 samtools 软件，可以使用 conda 安装。

比对多个文件

在数据分析中我们往往会有十几个甚至几百个数据，Snakemake 允许使用命名通配符（named wildcards） 概括规则。只需用通配符 {Sample} 替换第二个输入文件中的 A 和输出文件即可：

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

通过 {Sample} 即可达到匹配多个 fastq 文件的目的。

sam 排序

使用 samtools sort 命令对 bam 文件排序，对于的模块代码：

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}"

同样的定义了 samtools_sort 的规则模块，输入为 mapped_reads/{sample}.bam ，输出为 sorted_reads/{sample}.bam ，将会自动创建 sorted_reads 文件夹，-T 参数为输出文件名的前缀。

5建索引

给 bam 文件建立索引：

rule samtools_index:
    input:
        "sorted_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam.bai"
    shell:
        "samtools index {input}"

注意这里的 shell 命令里可以只含有 {input}，但是上面的 input 和 output 必须有。

6输出流程图

$ snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg

7Calling genomic variants

Snakemake 提供了 expand 函数来聚合多个文件，类似于列表推导式：

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

# 对应结果
["sorted_reads/A.bam", "sorted_reads/B.bam"]

当模式包含多个通配符时，该功能特别有用。例如：

SAMPLES = ["A", "B"]
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"]

使用 bcftools ：

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}"

这里输入文件有多个，fa、bam 和 bai 文件，我们可以使用 {input.fa}、{input.bam}的格式进行获取，也可以使用 {input[0]}、{input[1]} 来获取。

此外，输入文件使用了 expand 扩展函数，说明输入文件例如 bam 文件将是 多个 bam 文件，对应的 shell 命令应该是：

$ samtools mpileup -g -f \
           data/genome.fa \
           sorted_reads/A.bam sorted_reads/B.bam |
           bcftools call -mv - > {output}

到这里对应的流程图：

8使用自定义脚本

snakemake 还可以支持自定义脚本，如 python 和 R ，也可以直接写 python 代码。

示例：

rule plot_quals:
    input:
        "calls/all.vcf"
    output:
        "plots/quals.svg"
    script:
        "scripts/plot-quals.py"

这里定义了一个绘图的模块，使用 calls/all.vcf 文件用 plot-quals.py 脚本绘图，输出 plots/quals.svg 图片。

python 脚本代码如下：

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from pysam import VariantFile

quals = [record.qual for record in VariantFile(snakemake.input[0])]
plt.hist(quals)

plt.savefig(snakemake.output[0])

我们使用 snakemake.input和 snakemake.output 作为脚本的输入输出。

---

使用 R 脚本：

rule NAME:
    input:
        "path/to/inputfile",
        "path/to/other/inputfile"
    output:
        "path/to/outputfile",
        "path/to/another/outputfile"
    script:
        "scripts/script.R"

R 脚本代码：

do_something <- function(data_path, out_path, threads, myparam) {
    # R code
}

do_something(snakemake@input[[1]], snakemake@output[[1]], snakemake@threads, snakemake@config[["myparam"]])

这里需要使用双[[]] 来获取输入或输出。

9添加目标规则

到目前为止，我们总是通过在命令行指定目标文件来执行工作流。除了文件名，如果请求的规则没有通配符，Snakemake 还接受规则名称作为目标。因此，可以编写收集所需结果或所有结果的特定子集的目标规则。此外，如果在命令行中没有给出目标，Snakemake 将定义 Snakefile 的第一条规则作为目标。因此，最佳做法是在工作流的顶部设置一条规则，该规则将所有通常需要的目标文件作为输入文件：

rule all:
    input:
        "plots/quals.svg"

10所有完整流程代码

SAMPLES = ["A", "B"]


rule all:
    input:
        "plots/quals.svg"


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 plot_quals:
    input:
        "calls/all.vcf"
    output:
        "plots/quals.svg"
    script:
        "scripts/plot-quals.py"

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,代码已上传至QQ群文件夹，欢迎下载。

群二维码：

老俊俊微信：

知识星球：

所以今天你学习了吗？

欢迎小伙伴留言评论！

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！

 往期回顾 

◀...