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群文件夹，欢迎下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
欢迎小伙伴留言评论！
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀MeRIP-seq 数据分析之 homer 富集 motif
◀MeRIP-seq 数据分析之 homer 注释 peaks
◀MeRIP-seq 数据分析之 callpeak 及 peak 可视化
◀…
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
