接觸到一個(gè)生信分析項(xiàng)目碎罚，作者基于snakemake 串起了17個(gè)用python 腳本編寫的分析流程。自己在搭建的Linux 環(huán)境下纳像，測試了這些腳本荆烈。此文簡要介紹snakemake 的基礎(chǔ)知識(shí)和實(shí)際分析的細(xì)節(jié)。
snakemake 簡介
Snakemake工作流管理系統(tǒng)是一個(gè)用于創(chuàng)建可重復(fù)竟趾、可擴(kuò)展的數(shù)據(jù)分析的工具憔购。snakemake的工作流基于python語言來完成，具有可讀性強(qiáng)岔帽、能無縫地?cái)U(kuò)展到服務(wù)器玫鸟、集群、網(wǎng)格和云環(huán)境犀勒，而不需要修改工作流定義等特性屎飘。此外，Snakemake工作流可以包含所需軟件的描述贾费，這些軟件將自動(dòng)部署到任何執(zhí)行環(huán)境中钦购。
snakemake 的幾個(gè)特性
- 可讀性強(qiáng)：snakemake 通過制定每一條規(guī)則（rule）來定義每一個(gè)分析步驟，在規(guī)則中定義輸入（input）褂萧、輸出文件（output）押桃。不同規(guī)則的依賴性則自動(dòng)識(shí)別。
- 分析流程的模塊化：可以將數(shù)據(jù)分析拆分成不同模塊导犹，通過腳本或shell 命令完成凉袱。
- 可擴(kuò)展性：結(jié)合conda 環(huán)境和配置文件弯淘，可以輕松實(shí)現(xiàn)環(huán)境的導(dǎo)出，擴(kuò)展至其他環(huán)境下。
snakemake 語法基礎(chǔ)
安裝snakemake
建議通過conda 創(chuàng)建一個(gè)環(huán)境譬重，通過conda 完成安裝
https://snakemake.readthedocs.io/en/stable/getting_started/installation.html
創(chuàng)建 snakefile
rule bwa_map: # 定義規(guī)則名柳刮，也是分析的步驟名
input: # 指定輸入文件尝抖，多個(gè)輸入文件伴嗡，以逗號(hào)分隔
"data/genome.fa",
"data/samples/A.fastq"
output: # 指定輸出文件
"mapped_reads/A.bam"
shell: # 運(yùn)行shell 命令，有多行命令時(shí)，使用三引號(hào)纫骑。大括號(hào)實(shí)現(xiàn)自動(dòng)替換內(nèi)容
"bwa mem {input} | samtools view -Sb - > {output}"
若文件夾不存在時(shí)蝎亚，snakemake 會(huì)自動(dòng)創(chuàng)建文件夾。
運(yùn)行 snakefile
當(dāng)所在環(huán)境已安裝snakemake 后先馆，且該目錄下存在snakefile发框，執(zhí)行snakemake --cores，該命令會(huì)自動(dòng)匹配snakefile文件煤墙，并運(yùn)行其中的規(guī)則梅惯。snakemake 采取自頂向下的方式執(zhí)行分析步驟，也會(huì)識(shí)別其中的依賴關(guān)系仿野。
snakemake 執(zhí)行過程中有一些常見的參數(shù)铣减，能為分析過程中提供一些便利。
-s 指定Snakefile脚作，
-n 不真正執(zhí)行葫哗，
-p 輸出要執(zhí)行的shell命令
-r 輸出每條rule執(zhí)行的原因,默認(rèn)FALSE
-j 指定運(yùn)行的核數(shù)，若不指定球涛，則使用最大的核數(shù)
-f 重新運(yùn)行第一條rule或指定的rule
-F 重新運(yùn)行所有的rule劣针，不管是否已經(jīng)有輸出結(jié)果
snakemake 進(jìn)階實(shí)戰(zhàn)
聲明進(jìn)程數(shù)
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語法進(jìn)行寫，然后用在snakemake 文件中聲明亿扁，configfile: "FilePath/config.yaml"捺典。需要使用cofigfile 的信息時(shí)，通過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 文件內(nèi)容
samples:
A: data/samples/A.fastq
B: data/samples/B.fastq
基于snakemake 的snuupy 項(xiàng)目
https://github.com/ZhaiLab-SUSTech/snuupy/tree/master
該項(xiàng)目是一個(gè)對(duì)單細(xì)胞核測序數(shù)據(jù)的二代測序數(shù)據(jù)和三代測序數(shù)據(jù)的snakemake 的分析流程襟己，涉及17個(gè)python 腳本，完全通過python 完成各個(gè)分析步驟哄褒。
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 的部分內(nèi)容
configfile: "/public/home/liuzj/publicPipeline/snuupy/snakemake/config.yaml"
pipelineDir = config['pipelineDir']
# 對(duì)于包括多個(gè)規(guī)則的snakemake 文件，這一條rule all 是必須的煌张。這是一個(gè)特殊的規(guī)則呐赡，只有一個(gè)輸入文件，該文件是snakemake 分析流程的目標(biāo)文件骏融。
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 的命令行參數(shù)輸入來完成的，這主要通過一個(gè)python 模塊click 來完成档玻。在snuupy 文件中定義parseIllumina 函數(shù)需要的參數(shù)怀泊，在命令行中的解析這些參數(shù)，傳遞給執(zhí)行腳本误趴，完成分析步驟霹琼。
click 模塊主要包括兩個(gè)必要參數(shù)：
@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
"""
## 導(dǎo)入實(shí)際分析腳本
from scripts.parseIllumina import parseIllumina
parseIllumina(
secondBam, secondIndex, genomeIndex, windowSize, parsedIndex, useColumn
)
配置文件部分內(nèi)容
配置文件中也可以包括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/'
總結(jié)
在snuupy 項(xiàng)目中，一個(gè)分析步驟涉及多個(gè)輸入文件枣申，使用snakemake 流程具有獨(dú)特的優(yōu)勢(shì)售葡。而且分析過程中，多處涉及同一個(gè)分析工具忠藤，將其寫入配置文件挟伙，也是簡化輸入的方式。
因?yàn)樯婕暗胶芏鄠€(gè)分析步驟模孩，先將命令打印出來尖阔，再對(duì)個(gè)別步驟單獨(dú)測試，或者直接運(yùn)行特定snakemake rule榨咐，是比較好的測試思路介却。
參考
snakemake 官方文檔
使用Snakemake搭建分析流程
snakemake 使用筆記
https://github.com/ZhaiLab-SUSTech/snuupy/tree/master
Python Click 模塊
