介绍

凭借经济高效的全基因组测序和广泛公开的数据，新一代生物信息学方法进入了超高通量领域。例如，最近发表的覆盖度超过30x的 10 多个小麦基因组的组装或用泛基因组取代单一参考基因组。然而，这些空前数量的数据的下游处理和分析往往很复杂。生物信息学方法通常由多个步骤组成，利用大量公开可用的第三方软件进行生物分析，这些软件在 Unix 兼容操作系统上执行。然而，每个步骤的软件通常使用不同的编程语言，可用在几个版本中，或者随着时间的推移而变得过时。管道中顺序使用的不同软件包的交叉兼容性问题已通过多种方式解决。软件存储库采用了维护每个包的不同版本的策略，并允许通过独立的虚拟环境使用不同的安装。 Anaconda 是最流行的版本控制和软件存储库之一，可帮助用户安装和维护不同版本的 Python、PERL、R 和 Unix 包。

由于生物信息学的发展，数据大小和随之而来的计算需求增加了多个数量级。单个输入数据文件的范围通常高达数百 GB，处理它们可能需要在单个内核上花费数月的计算时间。大多数生物信息学软件都可以通过 pthreads（POSIX 线程）多线程进行扩展，它将计算工作分配到多个内核上。然而，现代数据集的处理通常会取代台式计算机的大部分计算资源。高性能计算 (HPC) 集群提供数百个内核、高内存和主要数据存储选项。虽然使用 HPC 集群是新一代生物信息学任务的最佳选择，但 HPC 集群的可用性较差。计算任务（称为作业）需要在 HPC 集群的工作节点上执行，这些工作节点采用独立的文件系统和环境。然后，这些作业需要在 HPC 集群的作业调度程序（例如 SLURM 工作负载管理器）中排队，以便在下一个可用的工作节点上运行。虽然作业调度通常由调度程序有效处理，但为作业分配的计算资源不能自动确定，必须由用户预定义。这给已经错综复杂的生物信息学管道工作流程增加了另一层复杂性。

为了简化链接生物软件和在 HPC 集群上的数百个输入文件上运行正确版本的过程，已经开发了多种管道自动化工具。 Snakemake 是 GNU Make 的演变，用于在 UNIX 系统上编译软件。虽然代码通常在自上而下的系统中运行，但 Snakemake 会根据已完成的依赖项执行步骤。一旦它们的依赖关系都得到满足，步骤就会启动 [7]。 Snakemake 还与 SLURM 等 HPC 集群作业调度程序交互，以自动提交作业。设置和运行 Snakemake 将在以下部分使用示例管道进行讨论。该示例将涉及使用 Bowtie2 [8] 将整个基因组测序文件与参考比对，使用 Samtools [9] 将 BAM 文件转换为 FASTQ 文件，并使用 Samtools 创建每个碱基的覆盖率表。

材料

输入文件

此示例的输入文件是 fastq 格式的测序读取：sample-A.R1.fq、sample-A.R2.fq、sample-B.R1.fq 和 sample-B。 R2.fq。这些应该放在名为 input_data 的数据目录中。
此外，应在工作目录中创建一个名为 snakefile 的新空白文件，稍后将对其进行修改

硬件

该协议旨在用于 HPC 集群。虽然具体要求因任务而异，但处理 WGS 数据通常需要 GB 到 TB 范围内的磁盘空间和至少具有八个内核的多个节点才能有效地扩展流程。

软件

所有软件都在 HPC 集群上的 UNIX 环境中运行。包管理器 Anaconda 用于安装 Snakemake。可以在以下位置获取 Anaconda：https://www.anaconda.com/products/individual 一旦 Anaconda 被初始化，它可以用来安装 Snakemake 和对各个步骤进行故障排除。

方法

安装软件

Anaconda

wget https://repo.anaconda.com/archive/Anaconda3-2020.07-Linux-x86_64.sh 
chmod +x Anaconda3-2020.07-Linux-x86_64.sh
binbash Anaconda3-2020.07-Linux-x86_64.sh

创建一个新的 Anaconda 环境，激活它，并安装 Snakemake。

conda create -n snakemake 
conda activate snakemake 
conda install -c bioconda snakemake

为了保持 Snakemake 和 HPC 集群之间的稳定交互，应该使用 Anaconda 的更新功能定期更新 Snakemake

使用 Snakemake 运行流程

Snakemake 使用蛇文件来指定在哪个数据文件上运行哪个命令。 snakefile 由在每个管道步骤执行的规则组成。每个规则指定输入文件、输出文件和使用输入文件生成输出文件的命令（在本例中为 bash shell）。

基本snakefile示例

rule bowtie2_alignment: 
    input: 
    R1=’input_data/sample-A_R1.fq’ 
    R2=’input_data/sample-A_R2.fq’ 
    output: 
    sam=’bowtie2_output/sample-A.sam’ 
    stats=’bowtie2_output/sample-A.stats’ 
    shell: 
    ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}’’’

规则 bowtie2_alignment 搜索两个输入文件 sample-A.R1.fq 和 sample-A.R2.fq。如果找到这两个文件，则规则可以执行并运行 Bowtie2 比对。除了执行命令的规则之外，每个 Snakemake 运行都需要一个“全部规则（rule all）”，它表示所有输入文件都已处理并且 Snakemake 可以关闭。 “Rule all”只接受输入文件，一旦找到“rule all”的所有输入文件，就会自动终止 Snakemake 运行。 “rule all”的输入文件可以通过python 样式列表传递。

添加rule all

rule all: 
    input: 
        [’bowtie2_output/sample-A.sam’, ’bowtie2_output/sample-A.stats’] 
    rule bowtie2_alignment: 
        input: 
            R1=’input_data/sample-A_R1.fq’ 
            R2=’input_data/sample-A_R2.fq’ 
            output:
            sam=’bowtie2_output/sample-A.sam’ 
            stats=’bowtie2_output/sample-A.stats’ 
            shell: 
            ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}’’’

rule all 和规则 bowtie2_alignment 保存在一个名为 snakefile 的文件中，现在可以执行 Snakemake：

$ Snakemake -j 4

Snakemake 现在将执行所有指定的作业，在四个本地内核上运行。

Snakemake 中的通配符

在基本示例中，Snakemake 将使用 sample-A.R1.fq 和 sample-A.R2.fq 作为输入文件。 Snakemake 无需手动指定所有输入文件，而是理解通配符输入。要在 Snakemake 中使用通配符，使用 glob_wildcards 函数创建通配符列表，然后在单独的规则中使用这些通配符。 glob_wildcards 函数适用于大多数 python 正则表达式模式，但具有清晰的文件结构使用户能够创建简单有效的通配符。

将 glob_wildcards 函数添加到 snakefile

IDS, = glob_wildcards(’input_data/{id}_R1.fq’)

glob_wildcards 现在将通过筛选与 glob_wildcards 输入匹配的所有文件来生成对象 IDS。这将在给定示例中找到两个文件：

input_data/sample-A.r1.fq 
input_data/sample-B.r1.fq.

在找到文件之后，glob_wildcards 提取替换变量部分 sample-A 和 sample-B，并自动将它们添加到 IDS 的 python 风格列表中。通配符列表 IDS 现在可以在 snakefile 中使用了，它的语法略有改变。

通配符 snakefile 示例

IDS, = glob_wildcards(’input_data/{id}_R1.fq’) 
rule all: 
    input: 
    expand[’aligned/{id}.sam’, id={IDS}] 
rule bowtie2_alignment: 
    input: 
    R1=’input_data/{id}_R1.fq’ 
    R2=’input_data/[2]_R2.fq’ 
    output: 
    sam=’bowtie2_output/[2].sam’ 
    stats=’bowtie2_output/[2].stats’ 
    shell: 
    ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}’’’

这里通配符有两种使用方式。首先，通配符是单独指定的，例如在规则 bowtie2_alignment 中。这将分别对每个匹配的通配符运行规则。或者，由于 Snakemake 接受 python 样式列表作为输入和输出变量，我们可以使用 expand 函数聚合 IDS 列表的所有变量，如 rule_all 所示。可以在优秀的 Snakemake 小插图中找到有关 Snakemake 通配符的综合指南：https://Snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards

日志和基准

单个 Snakemake 作业的日志和基准有助于排除故障和优化作业执行。两者都需要指定为 Snakemake 选项，并且必须通过在命令中指定 2 > {log} 将每个命令的 std.out 传递给日志对象。
日志和基准示例 snakefile

IDS, = glob_wildcards(’input_data/{id}_R1.fq’) 
rule all: 
    input: 
    expand[’aligned/{id}.sam’, id={IDS}] 
rule bowtie2_alignment: 
    input: 
    R1=’input_data/{id}_R1.fq’
    R2=’input_data/{id}_R2.fq’ 
    output: 
    sam=’bowtie2_output/{id}.sam’ 
    stats=’ bowtie2_output/{id}.stats’ 
    conda: ’path/to/conda-config.yaml’ 
    shell: 
    ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {log}’’’

Snakemake 中的版本控制

Snakemake 将不同软件的无限命令链接在一起的能力增加了软件或依赖项冲突的风险。为了避免冲突，Snakemake 可以使用 Anaconda 为每个规则创建一个临时的 conda 环境。
运行与 Anaconda 相关的 Snakemake 需要三件事。
1.创建一个.yaml文件，指定该规则的conda环境应该安装哪些软件包：

vim bowtie2_conda.yaml

2. 将需要的环境添加到bowtie2_conda.yaml。

channels: - bioconda dependencies: - bowtie2=2.4.2

这个 Anaconda 配置文件现在使用 conda 选项传递给 snakefile。
Anaconda snakefile 示例

IDS, = glob_wildcards(’input_data/{id}_R1.fq’) rule all: input: expand[’aligned/[2].sam’, id={IDS}] rule bowtie2_alignment: input: R1=’input_data/{id}_R1.fq’ R2=’input_data/{id}_R2.fq’ output: sam=’bowtie2_output/{id}.sam’ stats=’ bowtie2_output/{id}.stats’ conda: ’path/to/conda-config.yaml’ shell: ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}

3. 在 Snakemake 执行中使用 --use-conda。

Snakemake --j 4 --use-conda

执行Snakemake后，创建指定的Anaconda环境并执行工作。

HPC 环境中的 Snakemake

在前面的示例中，Snakemake 在本地系统上运行。虽然 snakefile 和数据结构在 HPC 系统上是相同的，但需要调整作业执行以适应 HPC 作业调度程序（本例中为 SLURM）。使用 Snakemake 在 SLURM 上安排作业的最简单方法是使用 Snakemake 配置文件。配置文件是 .yaml 格式的文件，其中包含传递给 Snakemake 执行的标志。由于 config.yaml 文件是 YAML 格式，因此需要调整使用的选项的语法，详见注释 2。在此示例中，配置文件将被称为 slurm，使用它的 Snakemake 调用将如下所示：

$Snakemake --profile slurm

Snakemake 现在将在 slurm/ 文件夹中查找 config.yaml 文件。可以将 slurm/ 文件夹的绝对路径传递给 --profile 标志，或者在预定义位置创建文件夹，默认位置为 HOME 目录中创建了 slurm/ 文件夹和 config.yaml 文件：

$ mkdir $HOME/.config/Snakemake/slurm 
$ touch $HOME/.config/Snakemake/slurm/config.yaml

此示例的 config.yaml 文件

jobs: 50 
cluster: ’sbatch -A ’HPC-account-name’ 
--cpus-per-task {resources.cpus} 
-p {resources.p} 
-n 1
 -t {resources.time} 
 --mem {resources.mem} 
 --cluster {resources.cluster} 
 -o ’slurm_logs/%j.{rule}_{wildcards}.out’ 
 -e ’slurm_logs/%j.{rule}_{wildcards}.err’ ’ 
    use-conda: true 
    keep-going: true 
    rerun-incomplete: true 
    default-resources: [cpus=20, mem=’58G’, time=’24:00:00’, p=’workq’]

通常，配置文件包含 Snakemake 标志，例如包含 srun 标志的集群选项，以通过 SLURM 控制作业提交。使用标志的细分：

jobs: 50

Snakemake 将向集群并行提交最多 50 个作业，并在达到 50 个作业后暂停提交。作业完成后，将提交更多作业以再次达到限制。

use-conda: true

启用 Snakemake conda 支持。对应于 --useconda 标志。

keep-going: true

默认情况下，一旦出现错误，Snakemake 将停止提交作业。如果 keep-going 设置为 true，独立于失败作业以外的作业仍将被提交。

rerun-incomplete: true

如果 Snakemake 运行中断，Snakemake 将删除不完整的文件，任何中断的作业将在下一次 Snakemake 执行时重新运行。

cluster:

集群标志控制 HPC 提交（本例中为 SLURM）。通常，SLURM 使用 sbatch 函数提交作业，sbatch 选项可以在 SLURM 文档中找到。

-A {HPC-user-account} 
--cpus-per-task {integer, number of cores requested for a job} 
-p {queue, e.g., ’workq’} 
-t {requested wall time in ’hh:mm:ss’ format} 
--mem {requested job memory} 
-o {path/to/job-stdout-file} 
-e {path/to/job-error-file} ’

集群参数可以是固定值，也可以传递给一个对象，该对象是示例中的资源对象。资源对象本身包含在 config.yaml 中指定的值。例如，t {resources.time} 会指示 Snakemake 在对象资源中查找时间值并将其传递给 SLURM 的 -t wall time 标志。通过添加资源选项，为 snakefile 中的每个规则单独创建资源对象。

HPC 资源 snakefile 示例

IDS, = glob_wildcards(’input_data/{id}_R1.fq’) 
rule all: 
input: 
    expand[’aligned/{id}.sam’, id={IDS}] 
rule bowtie2_alignment: 
    input: 
    R1=’input_data/{id}_R1.fq’ 
    R2=’input_data/{id}_R2.fq’ 
    output: 
    sam=’bowtie2_output/{id}.sam’ 
    stats=’ bowtie2_output/{id}.stats’ 
    conda: ’path/to/conda-config.yaml’ 
    resources: time = ’02:00:00’ 
    shell: 
    ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}’’’

资源选项指定 2 小时的时间，然后将其传递给 SLURM。为避免指定每次创建规则时传递给 SLURM 的所有资源标志，可以为 SLURM 标志指定默认值。这是通过将 defaultresources 选项添加到 config.yaml 来完成的，如上面的示例 config.yaml 文件所示。

构建 Snakemake流程

Snakemake 的核心功能是按顺序执行多个任务，从而生成自动化管道。这是通过 Snakefile 中的输入和输出文件链管理的。一旦规则的输入文件存在，规则就会被执行，这会产生任务链。例如，如果规则 B 使用规则 A 的输出文件作为输入文件，规则 B 将在规则 A 完成后执行。

IDS, = glob_wildcards(’input_data/{id}_R1.fq’) 
rule all: 
    input: 
    expand[’depth/{id}.per-base.bed’, id={IDS}] 
rule bowtie2_alignment: 
    input: 
    R1=’input_data/{id}_R1.fq’ 
    R2=’input_data/{id}_R2.fq’ 
    output: 
    sam=’bowtie2_output/{id}.sam’ 
    stats=’ bowtie2_output/{id}.stats’ 
    conda: ’path/to/bowtie2-conda-config.yaml’ 
    resources: time = ’02:00:00’ 
    shell: 
    ’’’bowtie2 -X 1000 -x reference-genome -1 {input.R1} -2 {input.R2} -S {output.sam} --end-to-end --sensitive -threads 20 2> {output.stats}’’’

rule sam2bam: 
    input: 
    sam=’aligned/{id}.sam’ 
    output:
    bam=’aligned/{id}.bam’, 
    benchmark: 
    ’benchmarks/{id}.sam2bam.benchmark.txt’ 
    log: 
    ’logs/{id}.sam2bam.log’, 
    conda: 
    ’path/to/samtools-conda-config.yaml’ 
    resources: 
    cpus=20, 
    time=’00:30:00’ 
    shell: 
    ’’’samtools view -bS -@ 20 {input.sam} > {output.bam}’’’ 
rule sort_bam: 
    input: 
    bam=’aligned/{id}.bam’, 
    output: 
    sorted_bam=’aligned/{id}.sorted.bam’, 
    conda: ’path/to/samtools-conda-config.yaml’ 
    log: ’logs/{id}.sort_bam.log’, 
    benchmark: ’benchmarks/{id}.sort_bam.benchmark’ 
    resources: 
    cpus=20, 
    time=’00:30:00’ 
    shell: 
    ’’’samtools sort -@ 20 {input.bam} -o {output.sorted_bam} 2> {log}’’’ 
rule samtools_depth: 
    input: 
    bam=’aligned/{id}.sorted.bam’, 
    output: 
    base=’depth/{id}.per-base.bed’, 
    conda:
    ’path/to/samtools-conda-config.yaml’ 
    log: 
    ’logs/{id}.samtools-depth.log’,
    benchmark: 
    ’benchmarks/{id}.samtools-depth.benchmark’ 
    resources: 
    cpus=1, 
    time=’01:00:00’ 
    shell: 
    ’’’samtools depth -a {input.bam} > {output.base} 2> {log}’’’

在示例中，样本序列使用 bowtie2 与参考序列比对，.sam 文件被转换为.bam 文件，.bam 文件被排序，并计算它们的深度。一旦计算出通过代表通配符指定的所有文件的深度，rule all 将终止 Snakemake 运行。

与上面概述的作业顺序相反，Snakemake 基于依赖关系管理作业执行。首先，规则所有检查是否存在深度文件。如果它们不存在，Snakemake 会检查它们是否在另一个规则的输出中，即示例中的规则 samtools_depth。同样，Snakemake 检查规则 samtools_depth 的输入文件是否存在，这些文件存在或被指定为另一个规则的输出。重复此操作，直到找到现有的输入文件，然后执行一系列规则以创建深度文件。

附加功能

Snakemake 提供了丰富的附加功能，包括工作流的图形表示（见注释 3 和 4）。本指南说明了在 HPC 环境中构建和执行 Snakemake 工作流的简单示例，但无法涵盖 Snakemake 的全部功能。要了解有关 Snakemake 特定附加功能的更多信息，请参阅 Snakemake 小插图（请参阅note 5）。

Note

1.虽然 Anaconda 有助于维护软件包和 Snakemake 本身的版本，但它不能自动执行此操作。维护 Snakemake 的最新版本对于与 HPC 集群上使用的作业调度程序进行稳定的交互是必要的。要使用 anaconda 更新 Snakemake，可以使用以下命令：

conda update -n envname snakemake

或者有时全新的 Snakemake 安装可以解决兼容性问题。

conda remove -n snakemake snakemake 
conda activate snakemake 
conda install snakemake

2.集群配置文件 config.yaml 可以包含 Snakemake 和 SLURM 选项。但是，这些以 YAML 格式传递，这略微改变了语法。例如，Snakemake 选项 --keep_going 变为 keep-going: true。 SLURM 运行参数在 cluster: config.yaml 的选项中指定，并且是 srun 参数。所有 srun 参数均可在以下位置找到：https://slurm.schedmd.com/srun
3.Snakemake 还可以创建设计工作流程的图形表示。这可以通过两种方式实现。通过使用 --dag 选项执行 Snakemake 并使用管道将其发送到图形设备，可以生成显示工作流如何处理每个样本的有向无directed acyclic graph (DAG)：

snakemake --dag | dot -Tpdf > dag.pdf

4..或者，可以使用 --rulegraph 选项生成仅显示规则的图表：

snakemake --rulegraph | dot -Tpdf > rulegraph.pdf

本指南概述了如何在 HPC 环境中自动化管道。虽然这些信息可能足以满足大多数项目目标，但 Snakemake 提供了许多附加功能，例如内置的 html 报告。要探索此附加功能，请参阅出色的 Snakemake 插图：https://Snakemake.readthedocs.io/en/stable/