Snakemake：簡(jiǎn)單好用的生信流程搭建工具

Snakemake是什么？

 

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

這時(shí)候就需要 Snakemake 閃亮登場(chǎng)了，因為snakemake可能是處理這類(lèi)問(wèn)題最簡(jiǎn)單、有效的一個(gè)方法了。比如將流程像下圖這樣搭建起來(lái)，清晰明了：

PS：這張流程圖也是 snakemake 通過(guò)命令自動(dòng)生成的奧，具體命令后面會(huì )有介紹。

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

snakemake的用法介紹

 

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

01操作環(huán)境

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

02軟件安裝

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

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

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

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

 

03數據準備

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

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

 

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

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

PS：建議將工作流文件命名為 py 文件，這樣你在寫(xiě) Python 代碼時(shí)會(huì )有語(yǔ)法高亮奧。

 

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

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

 

06使用通配符升級規則

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

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

 

07添加 sort 規則

比對完成之后需要使用 samtools 進(jìn)行排序，那我們就可以繼續寫(xiě)下一個(gè)規則：

 

08建立索引

 

09添加 calling 規則

 

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

1. 首先在工作流文件開(kāi)頭定義一個(gè)變量：

2. 然后使用 snakemake 的內置函數 expand：

PS：這個(gè) expand 函數其實(shí)是 Python 的列表推導式轉化而來(lái)，有興趣的同學(xué)可以自己研究一下。

最后，calling 規則如下：

10調用自己寫(xiě)的腳本

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

 

11最重要的一條規則：all

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

 

12生成可視化流程圖

 

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

 

 

13完整工作流

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

參考資料：

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

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