Snakemake是什么？

 

大多數(shù)生物信息初學者較早接觸的生物信息分析流程，大多數(shù)是用腳本語言（Perl、shell、R、Python等）將生物信息流程中的各個功能寫成一個個腳本，再用大腳本套小腳本的方式將整個流程打包封裝。這種方式對編程能力有很高的要求，否則很容易出現(xiàn)流程bug、無法監(jiān)控任務運行狀態(tài)、無法有效查看任務日志等問題。等實驗室終于出了個大牛將這些問題都解決了，可能大牛又要畢業(yè)了，新來的生信小白又只能看著層層疊疊的嵌套腳本“望洋興嘆”……

這時候就需要 Snakemake 閃亮登場了，因為snakemake可能是處理這類問題最簡單、有效的一個方法了。比如將流程像下圖這樣搭建起來，清晰明了：

PS：這張流程圖也是 snakemake 通過命令自動生成的奧，具體命令后面會有介紹。

Snakemake是一個工作流程管理系統(tǒng)。它是基于Python的、用于創(chuàng)建可重現(xiàn)和可擴展的數(shù)據(jù)分析的工具（當然現(xiàn)在也可以直接將它當做Python的一個模塊）。Snakemake所創(chuàng)建的流程還可以無縫擴展到服務器、集群和云環(huán)境等不同環(huán)境，當然前提是你需要提前將所需的軟件和依賴配置好，一起打包封裝在conda環(huán)境中。

snakemake的用法介紹

 

Snakemake 是基于 Python 的一款工具，所以它也繼承了 Python 語言簡單易讀、邏輯清晰、便于維護的特點，同時它還支持 Python 語法，非常適合新手用戶。snakemake 的基本組成單位叫“規(guī)則”，即 rule；每個 rule 里面又有多個元素（input、output、run等）。它的執(zhí)行邏輯就是將各個 rule 利用 input/output 連接起來，形成一個完整的工作流，當檢測到 input，就執(zhí)行相應 rule；檢測到 output，就跳過相應rule，根據(jù)這一規(guī)則，snakemake 還可以實現(xiàn)斷點續(xù)投。

01操作環(huán)境

本教程全程在Linux環(huán)境下完成，為了順利完成教程內(nèi)容，請準備一個Linux環(huán)境。如果是windows用戶，可以安裝一個虛擬機；如果是win10用戶，可以使用win10自帶的子系統(tǒng)（WSL）安裝Linux系統(tǒng)。

02軟件安裝

這里建議使用conda安裝軟件，可以輕松解決各種依賴問題。同時建議將本教程中用到的軟件安裝在一個新的conda環(huán)境中，以避免不同版本軟件間的依賴沖突。由于篇幅有限，conda的安裝方法不再贅述。

該命令創(chuàng)建了一個名為snakemake_env的conda環(huán)境，并在該環(huán)境中安裝了 snakemake 。

進入 snakemake_env 環(huán)境，查看上面要安裝的軟件，可以發(fā)現(xiàn)都已將安裝完成。

如果你的conda安裝很慢或者獲取不了最新版本的軟件，可以使用下面的命令安裝：

 

03數(shù)據(jù)準備

我們直接下載官網(wǎng)提供的測試數(shù)據(jù)：

解壓后得到一個 data 文件夾，目錄如下：

 

04創(chuàng)建工作流文件：Snakefile

當然工作流文件 Snakefile 文件名并不是固定的，位置也不是固定的。當你的工作流文件就在當前目錄下，且名稱正好是 Snakefile 時，就可以像上面的示例一樣，不用指定具體的工作流文件，snakemake 會自動調(diào)取當前路徑的名為 Snakefile 的文件去執(zhí)行。如果你要執(zhí)行其他路徑的 Snakefile 或者其他的文件名的工作流文件，可以使用 -s 參數(shù)，例如：

PS：建議將工作流文件命名為 py 文件，這樣你在寫 Python 代碼時會有語法高亮奧。

 

05創(chuàng)建第一條規(guī)則 bwa

我們繼續(xù) dry run ?一下，可以加上 -p 參數(shù)讓終端打印出 shell 運行的命令：

 

06使用通配符升級規(guī)則

上面的 bwa 規(guī)則僅僅比對了一個樣本，可是實際項目中有幾十上百的樣本時，我們就不能這樣直接寫樣本名來運行了，這樣很不 pythonic。snakemake 允許使用通配符（wildcard）來批量運行命令，示例如下：

其中，我們使用通配符 {sample} 來匹配所有樣品，花括號里的 sample 也可以使用其他字符。然后再次 dry run：

 

07添加 sort 規(guī)則

比對完成之后需要使用 samtools 進行排序，那我們就可以繼續(xù)寫下一個規(guī)則：

 

08建立索引

 

09添加 calling 規(guī)則

 

前面三步都是各個樣品批量、并行運行同樣的步驟，所以可以全部使用通配符 {sample} 完成匹配；但在變異檢測這一步需要將所有樣本的bam文件傳遞給一個命令而不再并行，這種方法就不在適用了。針對這種情況，snakemake 有它獨特的解決辦法：

1. 首先在工作流文件開頭定義一個變量：

2. 然后使用 snakemake 的內(nèi)置函數(shù) expand：

PS：這個 expand 函數(shù)其實是 Python 的列表推導式轉(zhuǎn)化而來，有興趣的同學可以自己研究一下。

最后，calling 規(guī)則如下：

10調(diào)用自己寫的腳本

當然，有時候在某個步驟想實現(xiàn)的功能并不能用簡單的幾句命令來實現(xiàn)，而且你自己也已經(jīng)有了一個現(xiàn)成的腳本了，那么你就可以直接調(diào)用該腳本實現(xiàn)相應功能：

 

11最重要的一條規(guī)則：all

前面我們已經(jīng)講過，snakemake 通過各個規(guī)則的輸入、輸出來確定執(zhí)行順序和依賴關系，snakemake 每次默認執(zhí)行第一個規(guī)則之后，會通過該規(guī)則的依賴關系（輸入）來決定接下來執(zhí)行哪個規(guī)則。所以，為了將所有規(guī)則都調(diào)用起來，我們需要在工作流文件最開始寫下第一個規(guī)則：all 規(guī)則，來調(diào)取各支線最后一條規(guī)則的輸出，這樣層層遞推，就可以將所有關聯(lián)規(guī)則都調(diào)用起來。示例如下：

 

12生成可視化流程圖

 

Snakemake 可以將整個工作流以流程圖的形式導出（結(jié)合 dot 命令），文件格式可以使png、pdf等。命令如下：

 

 

13完整工作流

對本文有疑問，請點擊上方按鈕聯(lián)系我們

關注百邁客云，隨時隨地學習

參考資料：

1. snakemake官方文檔 https://snakemake.readthedocs.io/en/stable/

2.snakemake github https://snakemake.github.io