基于Snakemake構建非模式生物轉錄組分析框架

李岸林 劉長寧

摘要：指出了在生命科學領域，非模式生物具有越來越重要的地位，對非模式生物的研究大幅拓寬了生命科學研究的廣度，促進發現更多新的科學問題，彌補模式生物在研究中的不足之處。傳統的轉錄組分析存在流程操作復雜、工作效率低、入門門檻高等問題。基于Snakemake的轉錄組分析框架，能夠運用數據標簽批量處理多個數據，無需使用者多次設置參數，大大提高了工作效率。利用Snakemake初步搭建轉錄組分析框架，在野生型和FUL基因型番茄樣品的RNA--Seq轉錄本數據進行測試。本研究工作為非計算機背景的研究人員進行轉錄組分析提供了方便，對相關轉錄組數據分析研究具有一定參考意義。

關鍵詞：分析框架;非模式生物;Snakemake轉錄組分析

中圖分類號：Q75文獻標識碼：A 文章編號：1674-9944（2019）14-0207-04

1引言

當前模式生物已經廣泛地應用于各個研究領域之中，為揭示中心法則以及基因和轉錄本等方面的一些基本的生物學規律做出了重要貢獻。經典模式生物雖然已被廣泛應用于各研究領域之中，但其在生物多樣性和功能性等方面存在著很多問題。對非模式生物的研究可以拓展生命科學研究的廣度，促進發現更多新的科學問題，彌補模式生物在研究中的不足之處。在遺傳和進化方面，非模式生物有著重要的研究價值。

在傳統的轉錄組數據分析流程中，一般分為多個步驟進行分析，在每個步驟中會包含一個或者多個分析軟件。研究人員需要通過編程腳本來銜接各個步驟，并且在分析的過程中會有大量重復性的數據分析操作，不僅對于研究人員要求有一定計算機技能要求，還會耗費大量的時間和精力在重復性的工作上。當前測序數據量飛速增長，在大數據時代下，需要更加強大的轉錄組數據分析工具，以Snakemake為基礎的分析框架可以大幅提高轉錄組數據分析效率，且簡單易用，為非計算機背景研究者提供了方便。

Snakemake是基于Python3編寫的一個流程管理工具，Snakemake將每個步驟簡化為一個rule，利用shell命令或者Python代碼進行文件輸入和輸出文件，還可以利用Python語言對軟件運行的過程進行管理。利用Snakemake搭建轉錄組分析可以提高轉錄組分數據分析效率。本文初步搭建基于Snakemake的轉錄組分析框架，并利用測試組數據對框架進行測試。

2框架搭建流程

在框架搭建過程中，需要安裝許多轉錄組分析相關的軟件，本框架利用Conda自動下載安裝需要的軟件，避免了手動下載安裝軟件耗時多、軟件版本錯誤等方面的問題。.Conda還可以配置相關依賴項，完成在本地計算機系統中運行環境的創建、保存、加載和切換等一系列的操作。在框架參數和軟件參數、文件路徑等設置方面，本框架通過配置文件config.yaml進行統一設置，無需重復對相同的參數進行設置。同時用戶也可以通過配置文件對框架進行設置管理，白定義轉錄組分析模式。在軟件選擇方面，對轉錄組分析相關軟件進行研究和對比，選擇性能較好、應用廣泛、穩定性強的主流軟件應用于框架中。在框架搭建中將所有軟件轉化為Snakemake中的rules，然后利用Python語言對rules進行優化和串聯并形成Snakefile文件。根據測序數據類型，將單端數據和雙端數據進行分別處理，在流程上分為有參轉錄組分析和無參轉錄組分析兩種類型。共搭建4個分析流程，包括Refer_SE.sk、Refer PE.sk、noRefer_SE.sk和noRefer_PE.sk，分別用于單端有參轉錄組分析、雙端有參轉錄組分析、單端無參轉錄組分析和雙端無參轉錄組分析。具體的分析流程如圖1所示。構建框搭建過程中所使用的Snakefile文件及其他配置可以在GitHub上獲取。

3數據測試與應用

本測的測試數據來自番茄（Lycopersicon esculen-turn MillJ）的轉錄組測序數據，為有參單端數據測試。番茄在成熟期間，其果實肉質、色素沉著、糖含量、香氣和營養品質等方面經歷復雜的生理和生化變化。深人探究其果實成熟的機制，有助于我們優化番茄品種的質量。在果安成熟的機制中，MADS基因是果實成熟中心調節的基因，FUL基因是MADS的基因中的一個關鍵的調節催熟的基因，將探究FUL基因的存在是否會影響不同階段的番茄基因差異表達。

3.1測試數據

本次測試數據的參考基因組來自于Sol GenomicNetwork網站中的International Tomato Genome Se-quencing Project中的最新的番茄全基因組，其版本號為SL3.0番茄基因組（GCA_000188115.3）。RNA-seq數據下載于National Center for Biotechnology In-formation（NCBI）數據庫，測序設備為Illumina HiSeq2000，數據量大約為20G。樣本包括野生型和FUL基因型兩種番茄樣品，同時前期和后期番茄果實樣品各3個，共12個樣品。樣本的具體信息如表1所示。

3.2數據質量控制結果

在數據質量控制部分，本框架利用FastQC軟件進行質量查看，Trimmomatic進行序列的修正和處理。通過Snakemake轉錄組分析框架，成功獲得12個測試數據的轉錄組質量分析結果，如表2所示。表中顯示數據的整體質量較好，質量分數中值（Mean Quality Scores）均大于20;不能分辨的堿基N幾乎為O%;讀段長度均為98到100bp，長度較為一致。

3.3基因的表達量計算

基因表達量分析最直接的方法是計算比對到每個基因上的讀段的數量，在轉錄組測序中通常稱其為量化（count）。通過使用基因組比對的方法，本框架的比對分析過程是利用HISAT2軟件進行的，比對分析獲得每條reads比對到基因組上的位置，以及每個基因比對的reads數量，這些信息通常以GTF格式文件進行儲存。HISAT2軟件在比對分析過程中，結合GTF文件中提供的轉錄本剪接信息，在一定程度上提高了讀段比對的準確性。通過比對，總共大約有147百萬讀段比對到了基因組上，對每個基因比對的reads數量進行統計，如表3所示。

3.4基因的聚類分析

聚類分析熱圖可以根據轉錄本的表達量高低，對轉錄本和樣本進行聚類分組，框架利用DEseq2進行聚類分析，R包來完成結果的可視化如圖2所示。在聚類分析熱圖中，通過顏色的深淺變化代表每個基因或轉錄本在不同樣品中的表達情況，紅色越深表達量越高，藍色反之。通過顏色的深淺就可以直觀的看到不同樣品或分組之間表達譜的差異。此外，將聚類（Clustering）的結果表示在熱圖的上方，以進化樹的結構展示。若基因或樣品的表達趨勢越接近，則在聚類結果中也更為接近，通常認為能聚到一起的基因簇或樣品簇有著更為接近的生物學特征，以此可以推斷未知基因的功能。

3.5基因的差異性表達

差異表達分析的目的是分析在不同的環境或者不同組織中的基因表達差異。通過對基因的表達量來發現是否存在基因的差異性表達。衡量基因差異性表達首先需要將基因的表達量進行標準化，接著利用統計模型進行差異性表達檢驗。差異性表達分析主要指標是Fold-change，Fold-change指的是兩個樣本之間表達量差異的倍數。p值（p-value）是反映了結果的可信度，當p值越小時，基因在不同條件下的差異性表達越可信。而adjusted P或q值（q-value）或是False dis-covery rate（FDR）也叫假陽性率（False Positive Rate）修正下的p值，通常使用q值或FDR來衡量表達差異的可信度。

通過框架DEseq2軟件計算得到Fold-change以及p值或q值后，R包繪制得到火山圖（Volcano Plot），如圖3所示。在生成的火山圖中，橫軸logFC代表Foldchange的log值，縱軸代表adjusted p值的loglO的值，通過對橫縱坐標閾值的限定，可以非常直觀篩選出在兩樣本間發生差異表達的基因。本框架配置文件中設置p值默認值為0.05。而Fold Change值只有大于2時才認為是顯著差異表達，所以橫坐標值logFC應大于1。通過這樣的篩選條件，可以篩選出顯著差異表達的相關基因作為研究目標。

4總結與展望

基于Snakemake搭建的轉錄組框架降低了流程對于數據格式轉化、腳本編寫的要求，可以簡便、高效、方便地處理轉錄本數據。本研究初步搭建了非模式生物分析的框架，利用本框架可以高效完成一整套轉錄組數據分析的流程，獲得相應的數據和可視化結果。基于此框架可以使轉錄組分析更容易進行操作，為非計算機背景的研究人員進行轉錄組分析提供了方便。基于Snakemake進行轉錄本數據分析框架的搭建，能夠為后續相關分析框架的構建提供參考。

在Snakemake框架的搭建過程中，提高相關參數的普適性是一個難題，特別是對于沒有參考基因序列的非模式生物，無法確定固定的參數，需要用戶不斷地進行嘗試。在框架的功能方面，現在只有轉錄組分析的基本功能，框架的功能的豐富性有所欠缺，在后期需要進一步豐富框架的功能。另一方面，需要增強框架應用的易操作性，盡量優化用戶操作步驟。在框架的軟件選擇方面，在后期會添加備選軟件供用戶選擇。在后續工作中，會不斷完善框架內容和參數調配，并且針對不同的非模式生物匹配相應的配套參數。對于后續更新的一些框架內容，會定期上傳到Github上以供下載。