snakemake使用教程¶

snakemake简介¶

Snakemake 工作流程管理系统是一种用于创建可重复和可扩展数据分析的工具，工作流程是通过一种基于 Python3 的人类可读语言来描述的。撰写的工作流可以无缝扩展到服务器、集群、网格和云环境，而无需修改工作流定义。最后，Snakemake 工作流程可以包含对所需软件的描述，这些软件将自动部署到任何执行环境中。

特点¶

并行运算，支持各种任务调度管理器(SLURM/SGE/PBS/LSF等)；

断点运行，多种断点检测机制;

流程易于控制，且可读性强，简单易懂，rule之间的输入和输出可以无缝结合；

计算内存控制，进程/线程/内存；

运行时间控制；

提供详细运行记录，详细记录rule下的系统/CPU运行时间和内存等；

可拓展性强，环境易于构建&管理，支持conda/docker/singularity等多种容器；

支持模块化构建，模型可以重复使用，同时也方便管理，类似于一个函数/包，需要的地方直接调用；

支持嵌入python代码进行逻辑处理，更加灵活；

......

基本元素¶

snakemake流程必须具备三个有效元素才能构建，需要悉知。

input files: 定义了输入文件
output files: 定义了输出文件
rules: 定义了输入文件生成输出文件的处理计算过程

snakemake安装¶

snakemake目前是一个非常成熟的python包，官方社区积极维护和更新(基本上每月会更新2+版本)，安装方式：

####以下安装方案选择自己最方便的即可
####pip安装github包
pip install git+https://github.com/snakemake/snakemake
###pip安装pypi包
pip install -i http://mirrors.aliyun.com/pypi/simple/ --trusted-host mirrors.aliyun.com snakemake
###conda安装，这里同时创建了一个独立的conda环境，你可以去掉环境创建，推荐独立出来，避免包版本“污染冲突”
mamba create -c conda-forge -c bioconda -n snakemake snakemake

对于snakemake版本的选择，建议不要安装老版本(eg:v5、v6等版本)，否则某些功能无法使用甚至出问题，例如：snakemake中如果需要使用conda环境，snakemake和mamba必须是最新版本或者高版本，亲测：snakemak-6.9.1+mamba-0.15.3可行，用V5版本就不行；目前snakemake已经v8.16.0，我推荐不要使用最新的，可能不稳定，使用\v7的版本就够了。

此外对于老版本的snakemake，流程报错代码定位不准确，进行debug很头疼，自己曾经被v6版本的snakemake折磨过，建议还是安装v7+的snakemake

理解rule¶

为了清楚的理解snakemake的rule的使用，个人认为以下三个关键原始必须要悉知。

1.rule定义¶

snakemake中的rule定义一个从输入文件到输出文件的操作&处理，因此我认为理解为任务/task更易于理解，rule可以分为两个大类，我定义为以下两种(官方没有命名)：

1. 全局rule：rule的名称是固定的，必须是all，通常一个流程中通常只定义一个all rule，且必须指定一个input关键字，全局rule的作用：汇总其他rule的输出结果，用于判断流程是否分析完成，其次会根据汇总文件的时间戳等信息检测哪些rule需要重新run；此外，我们也可以利用该机制，当我们想重跑某个分析时，或者修改了某个rule时，可以将对应rule的文件结果删除或者重命名，重新投递snakemake命令行，流程将会自动从对应的rule重跑以及续跑关联的rule任务。
2. 任务rule：一个 任务rule** 可以看作是一个操作步骤，比如序列质控、基因组比对等调用特定软件的步骤就可以定义一个 *任务rule* ，任务rule的命名是任意的，某些情况下，可以省略，表示一个匿名规则，但是推荐与当前rule下做的处理操作相符，例如：对bam文件排序，可以取名为：SortBam；与全局rule不用的是：任务rule必须指定input,output,run/shell/notebook/script(四个执行方式任选其一）等三个关键字内容，其余关键字可选指定；此外，不同任务rule之间的依赖关系是通过输入和输出文件名称进行匹配，隐式地进行处理。

全局rule和任务rule，以及任务rule之间的关系可以理解为下图：

flowchart LR
      x[input data]-- input -->A([rule A])
      A-- output注册登记 -->y[rule all]
      A-- output结果 -->B([rule B])
      B-- output注册登记 -->y
      B-- output结果 -->C([rule ...])
      C-- output注册登记 -->y

2. 通配符说明¶

1. 通配符是什么？

通俗解释：通常一批数据很少只有一个样本，不同批次的项目分析样本名称和分组等信息均是不一样的，因此每次分析的输入或输出指定为一种可变的，或者可替换的字符是很有必要的；这里的样本名和组名/比较组等信息/物种/参考基因组等都可以设置为通配符，用于动态扩展&匹配每次分析的动态变化信息。

其实通配符还有一种通俗理解，你可以将它理解为编程语言中的正则表达式的捕获组(python ,R，Perl，bash都支持) 2. 如何使用

通配符是一个wildcards是一个通配符对象，Snakemake gets wildcards from the parsing the input/output file names. 通过解析各个rule中的input和output元素。例如我们在input中定义了一个sample的通配符变量，之后就可以通过wildcards进行引用glob_wildcards可以定义全局通配符，其实就是针对某个特定的路径进行正则匹配。

glob_wildcards可以接受一个或多个通配符表达式，类似{pattern}，最后返回一个由list组成的tuple，所以如果只有一个变量的话，别忘了添加逗号（具体可以参考下面的例子）。还有一点就是glob_wildcards不支持bash中用的通配符。

通配符可以在input/output/params等关键字下用花括号进行定义，不能在run/shell等字段下定义，否则会报错；这里写一个简单示例：{sample}代表通配符

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

## Example 1 glob_wildcards自动获取一个通配符
#(SAMPLES,) = glob_wildcards("data/samples/{sample}.fastq")

## Exapple 2 glob_wildcards获取两个通配符
#(SAMPLES, LAYOUT, ...) = glob_wildcards("data/samples/{sample}.{layout}.fastq")

# Example 3 直接list定义
SAMPLES=["A","B"]

rule all:
    input:
        expand("mapped_reads/{sample}.bam", sample=SAMPLES),

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

3. 关键字说明¶

snakemake的关键字也非为两种类型：全局关键字和任务关键字(也可称为局部关键字，只能在一个rule内访问），所谓的全局关键字就是全局变量，可以在流程中直接使用&访问获取需要的信息，其中最关键全局关键字有如下:

• config：config关键字其实就是一个全局字典，该字典可以用JSON 和 YAML 格式的文件配置获取，允许使用配置文件来定义一些用户自定义变量的值，如样本信息、数据路径、软件参数，软件路径(不使用conda等环境时)和一些流程控制参数，参考基因组信息等等，可以使工作流程更加灵活；示例如下：

- wildcards:通配符对象；

任务关键字的通俗解释：个人认为snakemake的每一个rule其实就是一个字典(就是python中的dict)，字典中的key值就是如下的十几个关键字(表格第一列)，每一个key值对应的value可以是list或者dict；最后，所有的rule又形成一个大的字典，大字典的key值就是rule的名称。rule中常涉及到的关键词主要包括如下15+个，新版本可能还会增加，后续再更新讲解。

其他个人不常用rule关键字包括：template_engine/group等等，详细使用见官方文档

示例&参数注解¶

这里以fastp质控原始数据&统计为例：

import os
import json
import getpass
import pandas as pd
from os.path import join as opj
from snakemake.utils import min_version

#### snakemake不同版本的参数/环境都会存在微调，版本之间可能存在不兼容情况，这里可以设置开发测试时的snakemake版本
# min_version("7.23.1")

#### 导入配置文件信息，基因组，软件，输入输出，工作目录等都可以在里面定义
# configfile: os.path.join(workflow.basedir, "../config/config.yaml")

# # project dir
# PATH = config.get('root_dir')
# # output
# OUTPUT = []

#### 导入其他已经写好的模块
# # ---------------------- prepare data ---------------------- #
# include: os.path.join(PATH, "rules/common/prepare.smk")
# # ---------------------- include rules --------------------- #
# # include pipeline rules
# include: os.path.join(PATH, f"rules/{PIPENAME}/{PIPENAME}.smk")

#### 输出全局参数，便于记录提示，推荐开放该功能
# onstart:
#     if "verbose" in config and config["verbose"]:
#         print(" Workflow Parameters ".center(80, '-'))
#         print("Workflow:          ", PIPELINE)
#         print("Sample list:       ", DATA.index.tolist())
#         print("Number of samples: ", len(SAMPLES))
#         print("-" * 80, "\n")

#         print(" Environment ".center(80, '-'))
#         print("Username:          ", getpass.getuser())
#         print("Project Directory: ", os.path.abspath(PATH))
#         print("Work Directory:    ", config['workdir'])
#         print("Output Directory:  ", config['outdir'])
#         print("-" * 80, "\n")

#### 生成样本名称通配符
(SAMPLES,) = glob_wildcards("data/samples/{sample}.fastq.gz") 
extra_params="-q 20 -u 30 -n 15 -l 50"  ####外置参数定义方案1——全局定义

#### 定义一些外置的全局变量，后续可以在input/output/param,run，shell等字段下直接使用，直接使用：{xx}即可使用,xx是变量名
FASTP="/xxx/xxx/xx/fastp"
QCSummary_py="xxxxx/xxxx/xxx/xx/QCSummary.py" ###外置的脚本


####以下时rule的定义

#定义rule all，rule名称不能更改，用于定义整个流程的目标结果文件，其他rule的结果均可以在这个rule下注册，后续流程根据注册的结果进行检查是否完整、时间点等和重跑
rule all:                        
    input:                           #rule all必须指定input
        expand(fq="QC/{sample}.fq.gz",sample=SAMPLES),
        expand(html="QC/{sample}.html", sample=SAMPLES),
        expand(json="QC/{sample}.json", sample=SAMPLES),
        "QC/QC.summary.csv"

#定义一个fastp质控处理任务
rule ReadsQC:
  input:                           #定义rule输入文件，以下是以键值对形式定义，其中{sample}是通配符，指代每一个样本名称
    fq="data/samples/{sample}.fastq.gz"
  output:                          #定义rule输出结果，以下是以键值对形式定义，其中{sample}是通配符，指代每一个样本名称
    fq="QC/{sample}.fq.gz",
    html="QC/{sample}.html",
    json="QC/{sample}.json"
  threads: 8                       #指定当前rule的最大线程数
  resources：                      #指定当前任务资源的可用量   
    tmpdir="/data/tmp",             ##指定当前任务snakemake的临时目录，这里可以进行定义单独rule的零时目录，也可以在后续的命令行上进行统一定义
    mem_mb=2000,                    ##指定当前任务的内存上限，单位为MB，这里设置内存需求为2G
        # nvidia_gpu=1                  ###指定GPU数
  log:                             #定义rule日志文件，用于debug
    "QC/{sample}.ReadsQC.log"
  params:                          #定义额外参数，只要是input和output无法定义的变量参数，均可在params下定义
    extra_params="-q 20 -u 30 -n 15 -l 50"               ####外置参数定义方案2——局部定义
  shell:                           #定义运行的脚本命令，用shell方式提交
    #这里只有分析命令
    "{FASTP} -w {threads} -h {output.html} -j {output.json} --in1 {input.fq} --out1 {output.fq} {extra_params} > {log} 2>&1"
    #为了演示，这里模拟多行命令
    # '''
    #   {FASTP} -w {threads} -h {output.html} -j {output.json} --in1 {input.fq} --out1 {output.fq} {extra_params} > {log} 2>&1
    #   wc -l {output.fq}
    # '''
################################################################
  # run:                           #定义运行的脚本命令，用run方式提交
  #   shell("{FASTP} -w {threads} -h {output.html} -j {output.json} --in1 {input.fq} --out1 {output.fq} {extra_params} > {log} 2>&1")
  #   shell("wc -l {output.fq}")

#####汇总所有样本fastp的json文件，统计结果生成QC统计表
rule QCSummary:
  input:
    expand("QC/{sample}.json")
  output:
    "QC/QC.summary.csv"
  ##其他关键字没必要配置，因为当前rule只是做一个json汇总，不怎么耗费资源
   params:
     all_json_file=",".join(["QC/"+x+".json"for x in SAMPLES])
  log:
    "QC/QCSummary.log"
  shell:
    "{QCSummary_py} {params.all_json_file} {output} > {log} 2>&1"

#####统计每一个样本fastp的json文件，统计结果生成单独QC统计表
rule QCSummary1:
  input:
    rules.QCSummary.output.json   ####该方案等价于：QC/{sample}.json
  output:                               
    "QC/QC.{summary}.summary.csv"   #####这里时每一个样本均统计结果
  log:
    "QC/QCSummary.{sample}.log"     #####注意这里和QCSummary rule的区别
  shell:
    "{QCSummary_py} {params.all_json_file} {output} > {log} 2>&1"

特别提醒：对于rule的定义和撰写时，注意缩进和末尾的,,特别提醒：rule all中input只有一行时，记得末尾加,,多行时，最后一行可以不写,,其他行必须写,

流程项目结构¶

官方给的规范化流程目录结构：

├── workflow
│   ├── rules              #独立存放路径， 所谓的独立模块就是一个科重复使用的snakmake流程，其他模块可以直接用include参数调用 
|   │   ├── module1.smk    ###独立的snakemake模块，官方推荐明明是smk后缀，个人倾向于设置为py后缀，为了编辑器可以语法高亮
|   │   └── module2.smk
│   ├── envs               #环境配置文件存放路径，例如conda的环境配置文件
|   │   ├── tool1.yaml
|   │   └── tool2.yaml
│   ├── scripts            #外置脚本的存放路径，推荐是独立（即直接命令行给参数就可以运行的脚本）的脚本
|   │   ├── script1.py
|   │   └── script2.R
│   ├── notebooks          #notebook文件的存放路径
|   │   ├── notebook1.py.ipynb
|   │   └── notebook2.r.ipynb
│   ├── utils              #常用函数定义存放路径，这里需要与scripts区别一下，scripts时独立可运行的脚本，utils定义的是可重复使用的函数模块
|   │   ├── utils.py
|   │   ├── utils.R
│   ├── wrapper            #存放个性化定义或者官方给定的wrapper模块
|   │   ├── samtools
|   │   └── bwa
|   └── Snakefile          #主流程文件，snakemake命令行运行的流程文件
├── config                 #流程配置文件存储路径
│   ├── config.yaml        #主要的配置文件
│   └── some-sheet.tsv     #其他文件，可以在config.yaml中进行指定
├── results                #结果存储目录，可以不要
└── resources              #流程资源配置路径，可以存放资源配置文件，例如

规范化的文档还包括一些其他东西，例如类似wrapper的描述信息，其实每一个rule都应该配置一个描述信息，让接盘的同事或者其他人可以理解你的流程，wrap描述或者rule的描述一般命名为：meta.yaml，示例如下：

name: samtools depth
description: Compute the read depth at each position or region using samtools.
authors:
  - xxx
notes: |
  * The `extra` param allows for additional program arguments (not `-@/--threads` or `-o`).
  * For more information see, http://www.htslib.org/doc/samtools-depth.html

如何成为高级玩家？¶

1. 学习创建wrapper模块

The Snakemake Wrappers repository | Snakemake wrappers https://snakemake-wrappers.readthedocs.io/en/stable/ 2. 学习使用snakemake内置的API 3. 学习其他大佬的流程

这里给个示例：https://github.com/dxsbiocc/Workflow 4. 多实践

运行流程¶

snakemake命令行参数解释¶

snakemake的命令行参数非常多，也非常复杂，可能也只有官方开发人员能力灵活运用每一个参数了，常用参数及其解释如下：

--dry-run, --dryrun, -n       不执行，仅仅检查该流程是否存在依赖错误
--cores,-c                    指定整个流程同时最多使用的CPU核心数或者任务数，如果指定-c,但不加数字，将默认使用所有CPU核心数资源
--jobs,-j                     在云或者集群上同时运行的任务数
-p                            输出每个rule需要执行的shell命令
-r                            输出每条rule执行的原因
-s,--snakefile                指定snakmake配置文件，如果不指定该参数，默认使用的当前目录下的Snakefile文件，如果该文件不存在，将会报错
--config                      设置或覆盖工作流配置对象中的值

--forece,-f                   强制重新执行指定rule
--forceall,-F                 强制执行指定rule和其依赖rule
--forcerun,-R                 当修改规则时，强制执行指定的规则，-R all 重跑所有流程，-R QCSummary ，重跑QCSummary任务
--unlock                      流程报错或者被强行中断后，很多任务没完成，snakemake检测到相关结果目录不完整，会将相关相关目录上锁，需要使用--unlock解锁目录才能重新投递
--rerun-incomplete            流程报错后，指定该参数可以自动重跑和续跑             
--default-resources           修改全局资源信息，例如修改： --default-resources "tmpdir='/data/tmp'"
--profile                     snakemake 包含非常多的参数，如果都写在命令行将会非常不便，尤其是报错了，修改命令行很麻烦可以使用--profile 参数来设置一个配置文件，将参数写入到文件中方便使用和查看
--touch                       更新文件的时间戳（不会重新跑）
--list                        展示smk脚本中所能获得的所有rule
--dag                         生成流程逻辑框架图
--keep-going                  在某个rule运行失败后，继续运行其他的独立任务
--stats                       记录任务的执行状态，输出到指定文件，默认是输出到屏幕    

这里注意—cores和—job的区别：--cores经常在本地运行snakemake时用到，而--jobs则经常在需要提交任务到集群的情况中用到；这里这里额外提一下本地运行CPU核心数资源的设置：如果将--cores指定为10，但你的某一个 rule需要 8 个线程，此时该rule的并行任务将会是1，多余的两个CPU核资源，内部调度程序将尝试使用其他rule来占用多余CPU核，达到充分利用资源。

详细解释见：snakemake -h或者官方在线文档。

snakemake的运行步骤¶

Snakemake 工作流一般分为三个阶段执行：如下：

1. 在初始化阶段，定义工作流的文件被解析并实例化所有规则；
2. 在 DAG 阶段，通过填充通配符并将输入文件与输出文件匹配来构建所有作业的有向无环依赖图；
3. 在调度阶段，执行作业的 DAG，根据可用资源启动作业。

单节点投递¶

单节点投递也叫做本地投递，详细的运行模式可以如下：

• 运行前检查潜在错误，有助于发现错误
  Text Only

  snakemake -n
  snakemake -np
  snakemake -nr
  # --dryrun/-n: 不真正执行
  # --printshellcmds/-p: 输出要执行的shell命令
  # --reason/-r: 输出每条rule执行的原因
  snakemake --dag  | dit -Tpdf > dag.pdf   ##########生成流程逻辑框架图(有向无环图，就是类似于GO的有向无环图)，便于查看流程处理是否有问题
  

• 流程无误后，正常投递
  Text Only

  snakemake -s QC.smk -j 30 --configfile ./config.yaml -p --default-resources "tmpdir='/data/tmp'"
  ##--configfile ./config.yaml 指定自定义的配置文件，如果不指定，可以在代码中进行指定读取
  
  ###################### 其他修改参数示例############
  #为特定的规则重新指定线程数和资源
  snakemake --cores 4 --set-threads QCSummary=2
  

• 流程报错异常

流程在报错或者强制中断后，未执行完整的结果目录会被"锁上"，此时需要对目录进行解锁，参考命令如下

snakemake -s QC.smk -j 30 --configfile ./config.yaml -p --default-resources "tmpdir='/data/tmp'" --unlock

###调试修改报错后，直接投递原来的命令或者添加"--rerun-incomplete" 即可续跑
snakemake -s QC.smk -j 30 --configfile ./config.yaml -p --default-resources "tmpdir='/data/tmp'" 

snakemake -f
# --forece/-f: 强制执行选定的目标，或是第一个规则，无论是否已经完成
snakemake -F
# --forceall/-F: 也是强制执行，同时该规则所依赖的规则都要重新执行
snakemake -R some_rule
# --forecerun/-R TARGET: 重新执行给定的规则或生成文件。当你修改规则的时候，使用该命令

集群投递¶

• qsub投递
  Text Only

  snakemake --cluster "qsub -V -cwd -q 投递队列" -j 10
  # --cluster /-c CMD: 集群运行指令
  ## qusb -V -cwd -q， 表示输出当前环境变量(-V),在当前目录下运行(-cwd), 投递到指定的队列(-q), 如果不指定则使用任何可用队列
  # --local-cores N: 在每个集群中最多并行N核
  # --cluster-config/-u FILE: 集群配置文件
  

• 其他集群任务调度器投递

Cluster Execution — Snakemake 7.19.1 documentation https://snakemake.readthedocs.io/en/v7.19.1/executing/cluster.html

后续个人将会基于个人实际经验完善snakemake集群核云的投递配置。

断点检测机制¶

很多时候如果第一个rule的输入只是单个文件，输出只是固定输出（例如某个文本）情况下，每次运行以下snakemake命令时会发现流程会从头开始运行，

此时重跑的原因是：Params have changed since last execution，即发现第一个rule的输入参数发生了变化，但是我们第一个rule的输入明明没有改变，造成该现象的原因是因为在新版本的snakemake(v7.8+)中，他的rule重跑触发机制发生了变化，snakemake会有很多的重跑触发机制，默认所有的重跑机制都会运行，保证结果的一致性。

debug经验¶

• 个人debug经验
• 检查语法是否错误，特别是末尾的逗号是否加上了
• 检查log关键字是否都输入到了一个log文件
• 检查python语法错误
• 检查包的兼容性，尤其是snakemake的版本问题
• 官方FAQ

Frequently Asked Questions | Snakemake 8.16.0 documentation https://snakemake.readthedocs.io/en/stable/project_info/faq.html#id7

参考¶

1. https://snakemake.readthedocs.io/en/stable/
2. http://hpc.ncpgr.cn/app/093-snakemake/
3. https://mp.weixin.qq.com/s/9dTtpufm8KZDz7-OwRAPkg