接触到一个生信分析项目，作者基于snakemake 串起了17个用python 脚本编写的分析流程。自己在搭建的Linux 环境下，测试了这些脚本。此文简要介绍snakemake 的基础知识和实际分析的细节。
snakemake 简介
Snakemake工作流管理系统是一个用于创建可重复、可扩展的数据分析的工具。snakemake的工作流基于python语言来完成，具有可读性强、能无缝地扩展到服务器、集群、网格和云环境，而不需要修改工作流定义等特性。此外，Snakemake工作流可以包含所需软件的描述，这些软件将自动部署到任何执行环境中。
snakemake 的几个特性
- 可读性强：snakemake 通过制定每一条规则（rule）来定义每一个分析步骤，在规则中定义输入（input）、输出文件（output）。不同规则的依赖性则自动识别。
- 分析流程的模块化：可以将数据分析拆分成不同模块，通过脚本或shell 命令完成。
- 可扩展性：结合conda 环境和配置文件，可以轻松实现环境的导出，扩展至其他环境下。
snakemake 语法基础
安装snakemake
建议通过conda 创建一个环境，通过conda 完成安装
https://snakemake.readthedocs.io/en/stable/getting_started/installation.html
创建 snakefile
rule bwa_map: # 定义规则名，也是分析的步骤名
input: # 指定输入文件，多个输入文件，以逗号分隔
"data/genome.fa",
"data/samples/A.fastq"
output: # 指定输出文件
"mapped_reads/A.bam"
shell: # 运行shell 命令，有多行命令时，使用三引号。大括号实现自动替换内容
"bwa mem {input} | samtools view -Sb - > {output}"
若文件夹不存在时，snakemake 会自动创建文件夹。
运行 snakefile
当所在环境已安装snakemake 后，且该目录下存在snakefile，执行snakemake --cores，该命令会自动匹配snakefile文件，并运行其中的规则。snakemake 采取自顶向下的方式执行分析步骤，也会识别其中的依赖关系。
snakemake 执行过程中有一些常见的参数，能为分析过程中提供一些便利。
-s 指定Snakefile，
-n 不真正执行，
-p 输出要执行的shell命令
-r 输出每条rule执行的原因,默认FALSE
-j 指定运行的核数，若不指定，则使用最大的核数
-f 重新运行第一条rule或指定的rule
-F 重新运行所有的rule，不管是否已经有输出结果
snakemake 进阶实战
声明进程数
rule samtools_sort:
input:
"mapped_reads/{sample}.bam"
output:
"sorted_reads/{sample}.bam"
threads:
56
shell:
"samtools sort -T sorted_reads/{wildcards.sample} "
"-O bam {input} > {output}"
配置文件
配置文件可以用JSON或YAML语法进行写，然后用在snakemake 文件中声明，configfile: "FilePath/config.yaml"。需要使用cofigfile 的信息时，通过config[键]来获取。
包括config.yaml 的snakefile 文件
configfile: "config.yaml"
rule bcftools_call:
input:
fa="data/genome.fa",
bam=expand("sorted_reads/{sample}.bam", sample=config["samples"]), # config["samples"])
bai=expand("sorted_reads/{sample}.bam.bai", sample=config["smaples"])
output:
"calls/all.vcf"
shell:
'''
samtools mpileup -g -f {input.fa} {input.bam} |
bcftools call -mv - > {output}
'''
config.yaml 文件内容
samples:
A: data/samples/A.fastq
B: data/samples/B.fastq
基于snakemake 的snuupy 项目
https://github.com/ZhaiLab-SUSTech/snuupy/tree/master
该项目是一个对单细胞核测序数据的二代测序数据和三代测序数据的snakemake 的分析流程，涉及17个python 脚本，完全通过python 完成各个分析步骤。
Job counts:
count jobs
1 addGeneName
1 addPolyATag
1 addUnmappedBaseTag
1 all
1 barcodeAssignment
1 generateIlluminaWindow
1 generateMtx
1 generateNanoporeWindow
1 getMismatch
1 getSpliceInfo
1 minimapMappingPolished
1 minimapMappingRaw
1 parseIllumina
1 polishReads
1 polyAClusterDetected
1 runCellRanger
1 windowBlast
17
snuupy snakemake 的部分内容
configfile: "/public/home/liuzj/publicPipeline/snuupy/snakemake/config.yaml"
pipelineDir = config['pipelineDir']
# 对于包括多个规则的snakemake 文件，这一条rule all 是必须的。这是一个特殊的规则，只有一个输入文件，该文件是snakemake 分析流程的目标文件。
rule all:
input:
generateMtxFinished = f"{config['resultDir']}step16_generateMtx/generateMtxFinished.empty" # f'***'python 语法，格式化字符输出
rule parseIllumina:
input:
runCellRangerFinished = f"{config['resultDir']}step1_runCellRanger/runCellRangerFinished.empty"
output:
parseIlluminaResults = f"{config['resultDir']}step2_parseIllumina/parseIlluminaResults.index"
params:
genomeFai = config['genomeFai'],
useBarcodeGz = f"{config['resultDir']}step1_runCellRanger/test/outs/filtered_feature_bc_matrix/barcodes.tsv.gz",
useBarcode = f"{config['resultDir']}step1_runCellRanger/test/outs/filtered_feature_bc_matrix/barcodes.tsv",
runCellRangerBam = f"{config['resultDir']}step1_runCellRanger/test/outs/possorted_genome_bam.bam",
windowSize = 500,
gpu = "0"
threads:2
shell:
"""
cd {pipelineDir}
gzip -d -c {params.useBarcodeGz} > {params.useBarcode} &&
python ./snuupy/snuupy.py parseIllumina --bam {params.runCellRangerBam} --barcode {params.useBarcode} --genome {params.genomeFai} --window {params.windowSize} --parsed {output.parseIlluminaResults}
"""
以上文件中shell 命令是通过python 的命令行参数输入来完成的，这主要通过一个python 模块click 来完成。在snuupy 文件中定义parseIllumina 函数需要的参数，在命令行中的解析这些参数，传递给执行脚本，完成分析步骤。
click 模块主要包括两个必要参数：
@main.command("parseIllumina")
@click.option(***)
# snuupy.py
import click
@click.group()
def main():
pass
@main.command("parseIllumina")
@click.option("--bam", "secondBam", help="cellranger result")
@click.option(
"--barcode",
"secondIndex",
help="the filtered barcode list which is output by cellranger; format tsv NOT gz!!!",
)
@click.option("--genome", "genomeIndex", help="the fai format file")
@click.option(
"--genomeCounts",
"useColumn",
default=5,
type=int,
show_default=True,
help="chromosome counts",
)
@click.option("--window", "windowSize", type=int, default=500, help="window size")
@click.option("--parsed", "parsedIndex", help="the parsedIndex; end with .index")
def _parseIllumina(
secondBam, secondIndex, genomeIndex, windowSize, parsedIndex, useColumn
):
"""
\b
parse Illumina bam file and generate Illumina index
"""
## 导入实际分析脚本
from scripts.parseIllumina import parseIllumina
parseIllumina(
secondBam, secondIndex, genomeIndex, windowSize, parsedIndex, useColumn
)
配置文件部分内容
配置文件中也可以包括yaml 文件。
pipelineDir:
'/public/home/liuzj/publicPipeline/snuupy/'
resultDir:
'/public/home/liuzj/publicPipeline/snuupy/example/snakemakeResults/'
#inpuf files
rawNanoporeFa:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/nanopore.fa'
nanoporeWorkspace:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/workspace'
nanoporeSeqSummary:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/nanoporeSeqSummary.txt'
#genome files
geneAnnoBed:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/gene.bed'
geneAnnoRepreBed:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/araport11.representative.gene_model.bed'
geneNot12Bed:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/geneNotBed12.bed'
genomeFa:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/genome.fa'
genomeFai:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/genome.fa.fai'
usedIntron:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/ref/annoFile/select_intron_by_Cyto.id.txt'
#softwares path
minimap2Path:
'/public/home/liuzj/softwares/anaconda3/bin/minimap2'
picardPath:
'/public/home/liuzj/softwares/picard/picard.jar'
seqkitPath:
#cellranger configs
runCellRangerConfig:
'/public/home/liuzj/publicPipeline/snuupy/removeExonRegion/cellRangerParameter.yaml'
#reads path
illuminaFastqs:
'/public/home/liuzj/publicPipeline/snuupy/example/00_data/illuminaFastq/'
总结
在snuupy 项目中，一个分析步骤涉及多个输入文件，使用snakemake 流程具有独特的优势。而且分析过程中，多处涉及同一个分析工具，将其写入配置文件，也是简化输入的方式。
因为涉及到很多个分析步骤，先将命令打印出来，再对个别步骤单独测试，或者直接运行特定snakemake rule，是比较好的测试思路。
参考
snakemake 官方文档
使用Snakemake搭建分析流程
snakemake 使用笔记
https://github.com/ZhaiLab-SUSTech/snuupy/tree/master
Python Click 模块
