使用Snakemake搭建分析流程
目前已有的框架
A review of bioinformatics pipeline framework 的作者對已有的工具進行很好的分類
作者的看法：
- implicit慎璧，也就是Make rule語法更適合用于整合不同執(zhí)行工具
- 基于配置的流程更加穩(wěn)定，也比較適合用于集群分配任務(wù)吏奸。
最后作者建議是：
- 如果實驗室既不是純粹的生物學(xué)試驗（不需要workbench這種UI界面）晓勇，也不需要高性能基于類的流程設(shè)計毙籽， 不太好選役首， 主要原則是投入和產(chǎn)出比
- 如果實驗室進行的是重復(fù)性的研究，那么就需要對數(shù)據(jù)和軟件進行版本控制掺逼， 建議是 configuration-based pipelines
- 如果實驗室做的是探索性的概念證明類工作（exploratory proofs-of-concept)吃媒，那么需要的是 DSL-based pipeline。
- 如果實驗室用不到高性能計算機(HPC)吕喘，只能用云服務(wù)器赘那，就是server-based frameworks.
目前已有的流程可以在awesome-pipeline 進行查找。
就目前來看氯质，pipeline frameworks & library 這部分的框架中 nextflow 是點贊數(shù)最多的生物學(xué)相關(guān)框架募舟。只可惜nextflow在運行時需要創(chuàng)建fifo，而在NTFS文件系統(tǒng)上無法創(chuàng)建闻察，所以我選擇 snakemake , 一個基于Python寫的DSL流程框架拱礁。
環(huán)境準備
為了能夠順利完成這部分的教程，請準備一個Linux環(huán)境辕漂，如果使用Windows呢灶，則按照biostarhandbook(一)分析環(huán)境和數(shù)據(jù)可重復(fù) 部署一個虛擬機，并安裝miniconda3钉嘹。
如下步驟會下載所需數(shù)據(jù)鸯乃，并安裝所需要的軟件，并且啟動工作環(huán)境跋涣。
wget https://bitbucket.org/snakemake/snakemake-tutorial/get/v3.11.0.tar.bz2
tar -xf v3.11.0.tar.bz2 --strip 1
cd snakemake-snakemake-tutorial-623791d7ec6d
conda env create --name snakemake-tutorial --file environment.yaml
source activate snakemake-tutorial
# 退出當前環(huán)境
source deactivate
當前環(huán)境下的所有文件
├── data
│ ├── genome.fa
│ ├── genome.fa.amb
│ ├── genome.fa.ann
│ ├── genome.fa.bwt
│ ├── genome.fa.fai
│ ├── genome.fa.pac
│ ├── genome.fa.sa
│ └── samples
│ ├── A.fastq
│ ├── B.fastq
│ └── C.fastq
├── environment.yaml
└── README.md
基礎(chǔ)：一個案例流程
如果你編譯過軟件缨睡，那你應(yīng)該見過和用過make, 但是你估計也沒有仔細想過make是干嘛用的鸟悴。Make是最常用的軟件構(gòu)建工具，誕生于1977年宏蛉，主要用于C語言的項目遣臼，是為了處理編譯時存在各種依賴關(guān)系，尤其是部分文件更新后拾并，Make能夠重新生成需要更新的文件以及其對應(yīng)的文件揍堰。
Snakemake和Make功能一致，只不過用Python實現(xiàn)嗅义，增加了許多Python的特性屏歹，并且和Python一樣非常容易閱讀。下面將使用Snakemake寫一個變異檢測流程之碗。
第一步：序列比對
Snakemake非常簡單蝙眶，就是寫各種rule來完成不同的任務(wù)。我們第一條rule就是將序列比對到參考基因組上褪那。如果在命令行下就是bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam幽纷。 但是按照Snakemake的規(guī)則就是下面的寫法。
# 用你擅長的文本編輯器
vim Snakefile
# 編輯如下內(nèi)容
rule bwa_map:
input:
"data/genome.fa",
"data/samples/A.fastq"
output:
"mapped_reads/A.bam"
shell:
"""
bwa mem {input} | samtools view -Sb - > {output}
"""
解釋一下：這幾行定義了一個規(guī)則(rule)博敬，在這個規(guī)則下友浸，輸入(input)有兩個，而輸出(output)只有一個偏窝，在shell中運行命令收恢，只不過里面的文件都用{}形式替代。偽執(zhí)行一下:snakemake -np mapped_reads/A.bam檢查一下是否會出錯祭往，真實運行情況如下（不帶規(guī)則伦意，默認執(zhí)行第一個規(guī)則）:
第二步：推廣序列比對規(guī)則
如果僅僅是上面這樣子處理一個文件，還無法體現(xiàn)snakemake的用途硼补，畢竟還不如手動敲代碼來的方便驮肉。snakemake的一個有點在于它能夠使用文件名通配的方式對一類文件進行處理。將上面的A改成{sample},就可以將符合*.fastq的文件處理成*.bam.
rule bwa_map:
input:
"data/genome.fa",
"data/samples/{sample}.fastq"
output:
"mapped_reads/{sample}.bam"
shell:
"""
bwa mem {input} | samtools view -Sb - > {output}
"""
那么已骇，用snakemake -np mapped_reads/{A,B,C}.bam离钝，就會發(fā)現(xiàn)，他非常機智的就比對了B.fastq和C.fastq疾捍，而不會再比對一遍A.fastq, 也不需要你寫一堆的判斷語句去手動處理奈辰。
當然，如果你用touch data/samples/A.fastq改變A.fastq的時間戳乱豆，他就會認位A.fastq文件發(fā)生了改變奖恰，那么重復(fù)之前的命令就會比對A.fastq。
第三步：比對后排序
比對后的文件還需要進一步的排序，才能用于后續(xù)分析瑟啃，那么規(guī)則該如何寫呢论泛？
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}"
以之前的輸出作為輸出文件名，輸出到另一個文件夾中蛹屿。和之前的規(guī)則基本相同屁奏，只不過這里用到了wildcards.sample來獲取通配名用作-T的臨時文件的前綴sample實際名字。
運行snakemake -np sorted_reads/B.bam错负，你就會發(fā)現(xiàn)他就會非常智能的先比對再排序坟瓢。這是因為snakemake會自動解決依賴關(guān)系，并且按照依賴的前后順序進行執(zhí)行犹撒。
第四步： 建立索引和對任務(wù)可視化
這里我們再寫一個規(guī)則折联，對之前的排序后的BAM文件建立索引。
rule samtools_index:
input:
"sorted_reads/{sample}.bam"
output:
"sorted_reads/{sample}.bam.bai"
shell:
"samtools index {input}"
目前已經(jīng)寫了三個規(guī)則识颊，那么這些規(guī)則的執(zhí)行和依賴關(guān)系如何呢诚镰？ snakemake提供了--dag選項用于dot命令進行可視化
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
第五步：基因組變異識別
基因組變異識別需要整合之前所有的BAM文件，你可能會打算這樣寫
rule bcftools_call:
input:
fa="data/genome.fa",
bamA="sorted_reads/A.bam"
bamB="sorted_reads/B.bam"
baiA="sorted_reads/A.bam.bai"
baiB="sorted_reads/B.bam.bai"
output:
"calls/all.vcf"
shell:
"samtools mpileup -g -f {input.fa} {input.bamA} {input.bamB} | "
"bcftools call -mv - > {output}"
這樣寫的卻沒有問題祥款，但是以后每多一個樣本就需要多寫一個輸入清笨，太麻煩了。這里就體現(xiàn)出Snakemake和Python所帶來的特性了刃跛，我們可以用列表推導(dǎo)式的方法搞定抠艾。
["sorted_reads/{}.bam".format(sample) for sample in ["A","B"]]
進一步，可以在規(guī)則外定義SAMPLES=["A","B"]奠伪，則規(guī)則內(nèi)的輸入可以寫成bam=["sorted_reads/{}.bam".format(sample) for sample in SAMPLES]. 由于列表推導(dǎo)式比較常用跌帐，但是寫起來有點麻煩首懈，snakemake定義了expand進行簡化, 上面可以繼續(xù)改寫成expand("sorted_reads/{sample}.bam", sample=SAMPLES)
那么最后的規(guī)則就是
SAMPLES=["A","B"]
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}"
小練習(xí)： 請用snakemake生成當前的DAG圖绊率。
第六步：編寫報告
上面都是在規(guī)則里執(zhí)行shell腳本，snakemake的一個優(yōu)點就是可以在規(guī)則里面寫Python腳本究履，只需要把shell改成run滤否，此外還不需要用到引號。
rule report:
input:
"calls/all.vcf"
output:
"report.html"
run:
from snakemake.utils import report
with open(input[0]) as vcf:
n_calls = sum(1 for l in vcf if not l.startswith("#"))
report("""
An example variant calling workflow
===================================
Reads were mapped to the Yeast
reference genome and variants were called jointly with
SAMtools/BCFtools.
This resulted in {n_calls} variants (see Table T1_).
""", output[0], T1=input[0])
這里還用到了snakemake的一個函數(shù)最仑，report藐俺，可以對markdown語法進行渲染生成網(wǎng)頁。
第七步：增加目標規(guī)則
之前運行snakemake都是用的snakemake 目標文件名, 除了目標文件名外泥彤，snakemake還支持規(guī)則名作為目標欲芹。通常我們按照習(xí)慣定義一個all規(guī)則，來生成結(jié)果文件吟吝。
rule all:
input:
"report.html
基礎(chǔ)部分小結(jié)：
總結(jié)下學(xué)習(xí)過程菱父，知識點如下：
- Snakemake基于規(guī)則執(zhí)行命令，規(guī)則一般由
input, output,shell三部分組成。 - Snakemake可以自動確定不同規(guī)則的輸入輸出的依賴關(guān)系浙宜，根據(jù)時間戳來判斷文件是否需要重新生成
- Snakemake
以{sample}.fa形式進行文件名通配官辽，用{wildcards.sample}獲取sample的實際文件名 - Snakemake用
expand()生成多個文件名，本質(zhì)是Python的列表推導(dǎo)式 - Snakemake可以在規(guī)則外直接寫Python代碼粟瞬，在規(guī)則內(nèi)的
run里也可以寫Python代碼同仆。 - Snakefile的第一個規(guī)則通常是
rule all，因為默snakemake默認執(zhí)行第一條規(guī)則
進階：對流程進一步修飾
在基礎(chǔ)部分中裙品，我們完成了流程的框架俗批，下一步則是對這個框架進行不斷完善，比如說編寫配置文件市怎，聲明不同rule的消耗資源扶镀，記錄運行日志等。
第一步： 聲明所需進程數(shù)
對于一些工具焰轻，比如說bwa臭觉，多進程或者多線程運行能夠大大加速計算。snakemake使用threads來定義當前規(guī)則所用的進程數(shù)辱志，我們可以對之前的bwa_map增加該指令蝠筑。
rule bwa_map:
input:
"data/genome.fa",
"data/samples/{sample}.fastq"
output:
"mapped_reads/{sample}.bam"
threads:8
shell:
"bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
聲明threads后，Snakemake任務(wù)調(diào)度器就會在程序運行的時候是否并行多個任務(wù)揩懒。這主要和參數(shù)中的--cores相關(guān)什乙。比如說
snakemake --cores 10
由于總體上就分配了10個核心，于是一次就只能運行一個需要消耗8個核心的bwa_map已球。但是當其中一個bwa_map運行完畢臣镣，這個時候snakemaek就會同時運行一個消耗8個核心的bwa_map和沒有設(shè)置核心數(shù)的samtools_sort,來保證效率最大化。因此對于需要多線程或多進程運行的程序而言智亮，將所需的進程單獨編碼忆某，而不是硬編碼到shell命令中，能夠更有效的使用資源阔蛉。
第二步：配置文件
之前的SAMPLES寫在了snakefile弃舒，也就是意味這對于不同的項目，需要對snakefile進行修改状原，更好的方式是用一個配置文件聋呢。配置文件可以用JSON或YAML語法進行寫，然后用configfile: "config.yaml"讀取成字典颠区，變量名為config削锰。
config.yaml內(nèi)容為:
samples:
A: data/samples/A.fastq
B: data/samples/B.fastq
YAML使用縮進表示層級關(guān)系，其中縮進必須用空格毕莱，但是空格數(shù)目不重要器贩，重要的是所今后左側(cè)對齊测暗。上面的YAML被Pytho讀取之后，以字典保存磨澡，形式為
{'samples': {'A': 'data/samples/A.fastq', 'B': 'data/samples/B.fastq'}}
而snakefile也可以改寫成
configfile: "config.yaml"
...
rule bcftools_call:
input:
fa="data/genome.fa",
bam=expand("sorted_reads/{sample}.bam", sample=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}"
雖然sample是一個字典碗啄，但是展開的時候，只會使用他們的key值部分稳摄。
關(guān)于YAML格式的教程稚字，見阮一峰的博客：http://www.ruanyifeng.com/blog/2016/07/yaml.html
第三步：輸入函數(shù)
既然已經(jīng)把文件路徑都存入到配置文件中，那么可以進一步的改寫之前的bwa_map里的輸入部分厦酬。也就是從字典里面提取到存放的路徑胆描。最開始我就是打算這樣寫
rule bwa_map:
input:
"data/genome.fa",
config['samples']["{sample}"]
output:
"mapped_reads/{sample}.bam"
threads:8
shell:
"bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
畢竟"{sample}"從理論上應(yīng)該得到sample的名字。但是snakemake -np顯示出現(xiàn)錯誤
KeyError in line 11 of /home6/zgxu/snakemake-snakemake-tutorial-623791d7ec6d/Snakefile:
'{sample}'
這可能是{sample}的形式只能在匹配的時候使用仗阅，而在獲取值的時候應(yīng)該用基礎(chǔ)第三步的wildcards.sample形式昌讲。于是繼續(xù)改成config["samples"][wildcards.sample]。然而還是出現(xiàn)了錯誤减噪。
name 'wildcards' is not defined
為了理解錯誤的原因短绸，并找到解決方法，我們需要理解Snakemake工作流程執(zhí)行的一些原理筹裕，它執(zhí)行分為三個階段
- 在初始化階段醋闭，工作流程會被解析，所有規(guī)則都會被實例化
- 在DAG階段朝卒，也就是生成有向無環(huán)圖证逻，確定依賴關(guān)系的時候，所有的通配名部分都會被真正的文件名代替抗斤。
- 在調(diào)度階段囚企，DAG的任務(wù)按照順序執(zhí)行
也就是說在初始化階段，我們是無法獲知通配符所指代的具體文件名瑞眼，必須要等到第二階段龙宏，才會有wildcards變量出現(xiàn)。也就是說之前的出錯的原因都是因為第一個階段沒通過负拟。這個時候就需要輸入函數(shù)推遲文件名的確定烦衣，可以用Python的匿名函數(shù)歹河，也可以是普通的函數(shù)
rule bwa_map:
input:
"data/genome.fa",
lambda wildcards: config["samples"][wildcards.sample]
output:
"mapped_reads/{sample}.bam"
threads: 8
shell:
"bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
第四步：規(guī)則參數(shù)
有些時候掩浙，shell命令不僅僅是由input和output中的文件組成，還需要一些靜態(tài)的參數(shù)設(shè)置秸歧。如果把這些參數(shù)放在input里厨姚，則會因為找不到文件而出錯，所以需要專門的params用來設(shè)置這些參數(shù)键菱。
rule bwa_map:
input:
"data/genome.fa",
lambda wildcards: config["samples"][wildcards.sample]
output:
"mapped_reads/{sample}.bam"
threads: 8
params:
rg="@RG\tID:{sample}\tSM:{sample}"
shell:
"bwa mem -R '{params.rg}' '-t {threads} {input} | samtools view -Sb - > {output}"
寫在rule中的params的參數(shù)谬墙，可以在shell命令中或者是run里面的代碼進行調(diào)用今布。
第五步： 日志文件
當工作流程特別的大，每一步的輸出日志都建議保存下來拭抬，而不是輸出到屏幕部默，這樣子出錯的時候才能找到出錯的所在。snakemake非常貼心的定義了log,用于記錄日志造虎。好處就在于出錯的時候傅蹂，在log里面定義的文件是不會被snakemake刪掉，而output里面的文件則是會被刪除算凿。繼續(xù)修改之前的bwa_map.
rule bwa_map:
input:
"data/genome.fa",
lambda wildcards: config["samples"][wildcards.sample]
output:
"mapped_reads/{sample}.bam"
params:
rg="@RG\tID:{sample}\tSM:{sample}"
log:
"logs/bwa_mem/{sample}.log"
threads: 8
shell:
"(bwa mem -R '{params.rg}' -t {threads} {input} | "
"samtools view -Sb - > {output}) 2> {log}"
這里將標準錯誤重定向到了log中份蝴。
第六步：臨時文件和受保護的文件
由于高通量測序的數(shù)據(jù)量通常很大，因此很多無用的中間文件會占據(jù)大量的磁盤空間氓轰。而特異在執(zhí)行結(jié)束后寫一個shell命令清除不但寫起來麻煩婚夫，而且也不好管理。Snakemake使用temp()來將一些文件標記成臨時文件署鸡，在執(zhí)行結(jié)束后自動刪除案糙。
rule bwa_map:
input:
"data/genome.fa",
lambda wildcards: config["samples"][wildcards.sample]
output:
temp("mapped_reads/{sample}.bam")
params:
rg="@RG\tID:{sample}\tSM:{sample}"
log:
"logs/bwa_mem/{sample}.log"
threads: 8
shell:
"(bwa mem -R '{params.rg}' -t {threads} {input} | "
"samtools view -Sb - > {output}) 2> {log}"
修改之后的代碼，當samtools_sort運行結(jié)束后就會把"mapped_reads"下的BAM刪掉靴庆。同時由于比對和排序都比較耗時侍筛，得到的結(jié)果要是不小心被誤刪就會浪費大量計算時間，最后的方法就是用protected()保護起來
rule samtools_sort:
input:
"mapped_reads/{sample}.bam"
output:
protected("sorted_reads/{sample}.bam")
shell:
"samtools sort -T sorted_reads/{wildcards.sample} "
"-O bam {input} > {output}"
最后撒穷，snakemake就會在文件系統(tǒng)中對該輸出文件寫保護匣椰，也就是最后的權(quán)限為-r--r--r--, 在刪除的時候會問你rm: remove write-protected regular file ‘A.bam’?.
進階部分小結(jié)
- 使用
threads:定義不同規(guī)則所需線程數(shù)，有利于snakemake全局分配任務(wù)端礼，最優(yōu)化任務(wù)并行 - 使用
configfile:讀取配置文件禽笑，將配置和流程分離 - snakemake在DAG階段才會知道通配的具體文件名，因此在input和output出現(xiàn)的
wildcards就需要推遲到第二步蛤奥。 - 在
log里定義的日志文件佳镜，不會因任務(wù)失敗而被刪除 - 在
params定義的參數(shù)，可以在shell和run中直接調(diào)用 -
temp()中的文件運行結(jié)束后會被刪除凡桥，而protected()中的文件會有寫保護蟀伸，避免意外刪除。
高級：實現(xiàn)流程的自動部署
上面的分析流程都是基于當前環(huán)境下已經(jīng)安裝好要調(diào)用的軟件缅刽，如果你希望在新的環(huán)境中也能快速部署你的分析流程啊掏，那么你需要用到snakmake更高級的特性，也就是為每個rule定義專門的運行環(huán)境衰猛。
全局環(huán)境
我建議你在新建一個snakemake項目時迟蜜，都先用conda create -n 項目名 python=版本號創(chuàng)建一個全局環(huán)境，用于安裝一些常用的軟件啡省，例如bwa娜睛、samtools髓霞、seqkit等。然后用如下命令將環(huán)境導(dǎo)出成yaml文件
conda env export -n 項目名 -f environment.yaml
那么當你到了一個新的環(huán)境畦戒，你就可以用下面這個命令重建出你的運行環(huán)境
conda env create -f environment.yaml
局部環(huán)境
當然僅僅依賴于全局環(huán)境或許還不夠方库，對于不同的規(guī)則(rule)可能還有Python2和Python3的區(qū)別，所以你還得為每個規(guī)則創(chuàng)建環(huán)境障斋。
snakemake有一個參數(shù)--use-conda,會解析rule中的conda規(guī)則薪捍，根據(jù)其提供的yaml文件安裝特定版本的工具，以基礎(chǔ)第一步的序列比對為例配喳，
rule bwa_map:
input:
"data/genome.fa",
"data/samples/A.fastq"
output:
"mapped_reads/A.bam"
conda:
"envs/map.yaml"
shell:
"""
mkdir -p mapped_reads
bwa mem {input} | samtools view -Sb - > {output}
"""
隨后在snakemake執(zhí)行的目錄下創(chuàng)建envs文件夾酪穿，增加map.yaml, 內(nèi)容如下
name: map
channels:
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
- https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
- defaults
dependencies:
- bwa=0.7.17
- samtools=1.9
show_channel_urls: true
注意: YAML文件的name行不是必要的，但是建議加上晴裹。
那么當你用snakmake --use-conda執(zhí)行時被济，他就會在.snakemake/conda下創(chuàng)建專門的conda環(huán)境用于處理當前規(guī)則。對于當前項目涧团，該conda環(huán)境創(chuàng)建之后就會一直用于該規(guī)則只磷，除非yaml文件發(fā)生改變。
如果你希望在實際運行項目之前先創(chuàng)建好環(huán)境泌绣，那么可以使用--create-envs-only參數(shù)钮追。
由于默認情況下，每個項目運行時只會在當前的.snakemake/conda查找環(huán)境或者安裝環(huán)境阿迈，所以在其他目錄執(zhí)行項目時元媚，snakemake又會重新創(chuàng)建conda環(huán)境，如果你擔心太占地方或者環(huán)境太大苗沧，安裝的時候太廢時間刊棕，你可以用--conda-prefix指定專門的文件夾。
代碼總結(jié)
最后的代碼如下
configfile: "config.yaml"
rule all:
input:
"report.html"
rule bwa_map:
input:
"data/genome.fa",
lambda wildcards: config["samples"][wildcards.sample]
output:
temp("mapped_reads/{sample}.bam")
params:
rg="@RG\tID:{sample}\tSM:{sample}"
log:
"logs/bwa_mem/{sample}.log"
threads: 8
shell:
"(bwa mem -R '{params.rg}' -t {threads} {input} | "
"samtools view -Sb - > {output}) 2> {log}"
rule samtools_sort:
input:
"mapped_reads/{sample}.bam"
output:
protected("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=config["samples"]),
bai=expand("sorted_reads/{sample}.bam.bai", sample=config["samples"])
output:
"calls/all.vcf"
shell:
"samtools mpileup -g -f {input.fa} {input.bam} | "
"bcftools call -mv - > {output}"
rule report:
input:
"calls/all.vcf"
output:
"report.html"
run:
from snakemake.utils import report
with open(input[0]) as vcf:
n_calls = sum(1 for l in vcf if not l.startswith("#"))
report("""
An example variant calling workflow
===================================
Reads were mapped to the Yeast
reference genome and variants were called jointly with
SAMtools/BCFtools.
This resulted in {n_calls} variants (see Table T1_).
""", output[0], T1=input[0])
執(zhí)行snakemake
寫完Snakefile之后就需要用snakemake執(zhí)行待逞。snakemake的選項非常多甥角，這里列出一些比較常用的運行方式。
運行前檢查潛在錯誤：
snakemake -n
snakemake -np
snakemake -nr
# --dryrun/-n: 不真正執(zhí)行
# --printshellcmds/-p: 輸出要執(zhí)行的shell命令
# --reason/-r: 輸出每條rule執(zhí)行的原因
直接運行:
snakemake
snakemake -s Snakefile -j 4
# -s/--snakefile 指定Snakefile识樱，否則是當前目錄下的Snakefile
# --cores/--jobs/-j N: 指定并行數(shù)嗤无，如果不指定N，則使用當前最大可用的核心數(shù)
強制重新運行：
snakemake -f
# --forece/-f: 強制執(zhí)行選定的目標怜庸，或是第一個規(guī)則当犯，無論是否已經(jīng)完成
snakemake -F
# --forceall/-F: 也是強制執(zhí)行，同時該規(guī)則所依賴的規(guī)則都要重新執(zhí)行
snakemake -R some_rule
# --forecerun/-R TARGET: 重新執(zhí)行給定的規(guī)則或生成文件休雌。當你修改規(guī)則的時候灶壶，使用該命令
可視化：
snakemake --dag | dot -Tsvg > dag.svg
snakemake --dag | dit -Tpdf > dag.pdf
# --dag: 生成依賴的有向圖
snakemake --gui 0.0.0.0:2468
# --gui: 通過網(wǎng)頁查看運行狀態(tài)
集群執(zhí)行：
snakemake --cluster "qsub -V -cwd -q 投遞隊列" -j 10
# --cluster /-c CMD: 集群運行指令
## qusb -V -cwd -q， 表示輸出當前環(huán)境變量(-V),在當前目錄下運行(-cwd), 投遞到指定的隊列(-q), 如果不指定則使用任何可用隊列
# --local-cores N: 在每個集群中最多并行N核
# --cluster-config/-u FILE: 集群配置文件
