趙 芳,高 靜,劉振羽,王永軍,鄔學(xué)敏

(內(nèi)蒙古農(nóng)業(yè)大學(xué)計算機(jī)與信息工程學(xué)院內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)大數(shù)據(jù)研究與應(yīng)用重點實驗室,中國內(nèi)蒙古呼和浩特010018)

轉(zhuǎn)錄組差異性表達(dá)分析是RNA測序數(shù)據(jù)分析的核心內(nèi)容,也是轉(zhuǎn)錄組研究的基本目標(biāo)之一。其通過分析不同條件下的轉(zhuǎn)錄組數(shù)據(jù),進(jìn)而識別差異基因或異構(gòu)體,對于揭示基因功能和調(diào)控規(guī)律具有重要意義[1]。但該過程是一個多步迭代的復(fù)雜流程,涉及多軟件、多步驟,軟件之間還具有復(fù)雜的數(shù)據(jù)依賴關(guān)系。傳統(tǒng)方式下,研究人員手動執(zhí)行各環(huán)節(jié),前一任務(wù)的輸出即為后一任務(wù)的輸入,只有前一任務(wù)執(zhí)行完畢后才可串行執(zhí)行后續(xù)依賴任務(wù)[2]。而且,各軟件的輸入輸出格式不盡相同,均需研究人員進(jìn)行手動的格式轉(zhuǎn)換。該過程中,復(fù)雜的命令行操作以及大量的參數(shù)配置很大程度上增加了轉(zhuǎn)錄組差異性表達(dá)分析的難度,導(dǎo)致效率低、出錯率高、可重復(fù)性差。此外,環(huán)境的不易遷移性也進(jìn)一步導(dǎo)致數(shù)據(jù)分析問題成為生物信息學(xué)的研究瓶頸。

近年來,針對轉(zhuǎn)錄組差異性表達(dá)分析,人們開發(fā)了一些高性能的自動化分析框架,如system-PipeR[3]、SAMSA[4]、miRPursuit[5]。這些框架實現(xiàn)了轉(zhuǎn)錄組差異性分析的自動化,一定程度上克服了傳統(tǒng)方式需大量手動操作的弊端。但在框架部署與初始配置時,僅提供命令行或特定的編程語言接口,仍需用戶進(jìn)行較多的手動操作,不利于缺乏編程基礎(chǔ)的生物信息人員使用。此外,移植性較差,需用戶在不同的環(huán)境下進(jìn)行多次重復(fù)性部署,且用戶無法按需集成新的工具,無法滿足研究人員的實際分析需求。

因此,文中針對以上低效重復(fù)的手動操作、環(huán)境的不可遷移以及無法按需集成等問題,提出了一種基于Docker容器技術(shù)的轉(zhuǎn)錄組差異性表達(dá)分析方法,形成一套可重復(fù)、易移植、高效自動化的輕量級分析框架。

1 概念定義與原理

1.1 Docker容器技術(shù)

Docker是基于Go語言實現(xiàn)的云開源項目,其主要目標(biāo)是通過對應(yīng)用組件的封裝、分發(fā)、部署、運行等生命周期的管理,使用戶的App(Web應(yīng)用或數(shù)據(jù)庫應(yīng)用等)及其運行環(huán)境能夠做到“一次封裝,到處運行”。Docker具有部署速度快、開發(fā)測試敏捷、提高系統(tǒng)利用率、降低資源成本、跨環(huán)境可移植性好等優(yōu)點。

因此,本研究首次將云開源項目Docker虛擬化技術(shù)應(yīng)用于生物信息學(xué)領(lǐng)域,用于解決復(fù)雜分析的不可重復(fù)性、不易遷移性以及無法按需集成等核心難題,使經(jīng)驗較少或非生物信息專業(yè)的研究人員執(zhí)行復(fù)雜的流程分析時就像運行一個簡單的生物信息學(xué)工具一樣。

1.2 最佳實踐流程

轉(zhuǎn)錄組差異性表達(dá)分析流程較為復(fù)雜,涉及多步驟的處理,如序列比對、轉(zhuǎn)錄本組裝、轉(zhuǎn)錄本合并等。各處理過程需借助不同的工具軟件來實現(xiàn),但目前各步驟的處理軟件層出不窮,軟件之間復(fù)雜的數(shù)據(jù)依賴關(guān)系使得研究人員往往存在著軟件選擇難的問題,面臨著各軟件相互銜接難的困境[6]。

2016年,Ghosh和Cean針對轉(zhuǎn)錄組差異性表達(dá)分析提出了一套最佳實踐流程[7],該流程在提出的幾年里被生物信息學(xué)領(lǐng)域?qū)＜宜J(rèn)可并延續(xù)至今仍被廣泛使用。因此,文中將該最佳實踐流程Bowtie2[8]-TopHat2[9]-Cufflinks[10]-cummeRbund[11]與數(shù)據(jù)預(yù)處理環(huán)節(jié)相整合,在已構(gòu)建好參考基因庫的前提下,將轉(zhuǎn)錄組差異性表達(dá)分析主要分為8個處理步驟,形成圖1所示完整流程。

圖1 轉(zhuǎn)錄組差異性表達(dá)分析處理流程Fig.1 The transcriptome differential expression analysis process

1)數(shù)據(jù)預(yù)處理:首先,用SRA-Toolkit(https://hpc.nih.gov/apps/sratoolkit.html)將sra格式的原始測序數(shù)據(jù)轉(zhuǎn)換為fastq格式;其次,通過FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc)對測序數(shù)據(jù)進(jìn)行質(zhì)量檢測;最后,通過Trimmomatic質(zhì)量檢測[12]將質(zhì)量低或帶有接頭(a-dapters)的測序數(shù)據(jù)進(jìn)行修剪和過濾,進(jìn)而得到高質(zhì)量測序數(shù)據(jù)。

2)序列比對:TopHat2(依賴Bowtie2,比對過程中會調(diào)用Bowtie2)將測序數(shù)據(jù)與參考基因庫進(jìn)行比對,輸出結(jié)果文件為accepted_hits.bam。

3)轉(zhuǎn)錄本組裝:Cufflinks將TopHat2/Bowtie2比對結(jié)果accepted_hits.bam及參考基因庫作為輸入文件構(gòu)建轉(zhuǎn)錄本,并計算出各轉(zhuǎn)錄本FPKM值。

4)轉(zhuǎn)錄本合并:Cuffmerge將多個轉(zhuǎn)錄本合并成一套轉(zhuǎn)錄本集合,將Cufflinks生成的gtf文件融合為一個更加全面的注釋文件merged.gtf。

5)表達(dá)水平量化:Cuffquant對單個bam文件表達(dá)水平進(jìn)行定量,用于后續(xù)Cuffdiff分析。

6)表達(dá)水平標(biāo)準(zhǔn)化:將Cuffquant的輸出結(jié)果作為Cuffnorm的輸入文件,對轉(zhuǎn)錄組標(biāo)準(zhǔn)化后的表達(dá)水平進(jìn)行計算,可得到一系列可比較的基因、轉(zhuǎn)錄組、CDS組和TSS組的表達(dá)值。

7)差異表達(dá)分析:通過Cuffdiff計算不同條件下轉(zhuǎn)錄本表達(dá)水平的顯著性差異,進(jìn)而尋找差異基因。

8)結(jié)果可視化:采用cummeRbund軟件對Cuffdiff輸出結(jié)果以熱圖形式進(jìn)行可視化呈現(xiàn),以更直觀的方式描述目標(biāo)差異基因。

轉(zhuǎn)錄組差異性表達(dá)分析流程中所需軟件如表1所示。

2 基于Docker的轉(zhuǎn)錄組差異性表達(dá)分析方法

2.1 層次構(gòu)建

基于Docker容器技術(shù)的轉(zhuǎn)錄組差異性表達(dá)分析方法形成的整體框架層次主要分為基礎(chǔ)層、核心層、調(diào)度層以及應(yīng)用層(圖2)。基礎(chǔ)層除硬件設(shè)施與操作系統(tǒng)外,還需為上層提供運行環(huán)境的Docker引擎;核心層為轉(zhuǎn)錄組差異性表達(dá)分析容器,屬于核心功能組件,可按需集成不同的復(fù)雜分析流程;調(diào)度層即通過多腳本嵌套觸發(fā)對服務(wù)、資源的協(xié)調(diào)和控制;應(yīng)用層即為用戶接口層。該框架各層次之間協(xié)同配合,進(jìn)而實現(xiàn)高效自動化的轉(zhuǎn)錄組差異性表達(dá)分析。

圖2 整體框架層次圖Fig.2 Framework hierarchy diagram

2.2 核心層

鑒于Docker容器可重復(fù)、配置部署快速簡單、易于遷移、易于擴(kuò)展的輕量級特性,文中將1.2節(jié)轉(zhuǎn)錄組最佳實踐流程涉及的軟件及其依賴環(huán)境封裝預(yù)配置到Docker容器中,形成一個輕量級自動化的轉(zhuǎn)錄組差異性表達(dá)分析容器,并將其命名為DEA-container(differential expression analysis container)。容器內(nèi)部主要通過多腳本聯(lián)動的方式觸發(fā)啟動DEA-container,控制和協(xié)調(diào)各組件的執(zhí)行(核心腳本信息如表2所示),從而實現(xiàn)轉(zhuǎn)錄組差異性表達(dá)的自動化分析。

表1 最佳實踐流程軟件列表Table 1 List of best practice softwares

表2 核心腳本列表Table 2 The list of core script

此外,DEA-container還為用戶提供了友好的Web圖形化界面。我們將全過程封裝整合,僅對外提供單一的數(shù)據(jù)輸入和輸出端點。用戶僅需將目標(biāo)分析序列提交。全過程無需人工干預(yù),經(jīng)后臺服務(wù)器分析完畢后會自動將分析結(jié)果文件通過Web界面反饋至用戶,供用戶下載。

2.3 多腳本聯(lián)動過程

根據(jù)整體的框架構(gòu)建,各層次之間的協(xié)調(diào)配合需要通過多腳本聯(lián)動方式對轉(zhuǎn)錄組差異性表達(dá)分析的自動化進(jìn)行實現(xiàn),具體實現(xiàn)分為以下處理流程,其執(zhí)行過程如圖3所示。

1)數(shù)據(jù)獲取:將用戶基于Web圖形界面提交的任務(wù)及分析參數(shù)形成任務(wù)表單,采用get_input.sh腳本讀取用戶數(shù)據(jù)及相關(guān)分析參數(shù)并將數(shù)據(jù)存儲于指定位置。

2)初始化檢測:通過init_env.sh腳本對宿主主機(jī)Docker服務(wù)進(jìn)行狀態(tài)監(jiān)測,觸發(fā)啟動DEA-container并對環(huán)境進(jìn)行初始化,主要完成DEA-container配置參數(shù)、最佳實踐流程組件以及任務(wù)數(shù)據(jù)狀態(tài)的檢測。

3)設(shè)置分析參數(shù):對用戶提交的任務(wù)表單進(jìn)行分析參數(shù)的讀取,即分析數(shù)據(jù)類型,包括分析物種以及物種的名稱。通過set_parameters.sh腳本對目標(biāo)參考基因組、任務(wù)數(shù)據(jù)的存放位置、輸入輸出參數(shù)進(jìn)行指定。

4)差異性表達(dá)分析:參數(shù)設(shè)置完畢后,通過pipeline.sh腳本觸發(fā)最佳實踐分析流程,分別進(jìn)行數(shù)據(jù)預(yù)處理、序列比對、轉(zhuǎn)錄本組裝、轉(zhuǎn)錄本合并、差異表達(dá)分析等核心處理步驟,最終實現(xiàn)轉(zhuǎn)錄組差異性表達(dá)的自動化分析,以更直觀、專業(yè)的熱圖方式描述差異基因。

5)結(jié)果反饋:將結(jié)果文件及其熱圖通過Web圖形化界面呈現(xiàn)給用戶。同時,為保證研究人員可進(jìn)行后續(xù)研究與分析,通過pull_output.sh腳本對核心結(jié)果文件進(jìn)行提取并以下載列表的形式反饋至用戶,方便下載和進(jìn)行二次處理分析。

以上各環(huán)節(jié)均存在檢錯、糾錯機(jī)制,當(dāng)某一功能環(huán)節(jié)檢測報錯后會選取最近點任務(wù)進(jìn)行斷點恢復(fù),這種實時的狀態(tài)檢測與糾錯可有效避免任務(wù)的中斷和失敗。在處理過程中,各腳本均會對上一步驟執(zhí)行狀態(tài)進(jìn)行檢測,并根據(jù)不同的返回值進(jìn)行相應(yīng)的操作。若返回值為1,表明上一步驟執(zhí)行異常,存在錯誤,程序自動進(jìn)入糾錯環(huán)節(jié),將上一步驟產(chǎn)生的結(jié)果文件刪除后重新執(zhí)行該步驟,運行完畢會再次檢測,檢測通過后才可執(zhí)行下一處理步驟。若返回值為0,則表示上一步驟執(zhí)行成功,可以順序執(zhí)行下一處理步驟。

圖3 多腳本聯(lián)動過程Fig.3 The process of scripts linkage

本研究向用戶提供兩種服務(wù)界面:1)用戶圖形界面:該方式僅需用戶將sra格式的原始序列通過Web用戶圖形界面進(jìn)行提交,所需的計算資源由該平臺的后臺服務(wù)器提供,全過程無需用戶干預(yù),處理完畢后分析結(jié)果會通過Web界面反饋至用戶并供用戶下載;2)命令行界面:本研究發(fā)布了DEA-container以及相關(guān)的一鍵化部署腳本,用戶將其下載至本地并通過一鍵化部署腳本可在任意物理主機(jī)基于命令行方式進(jìn)行特定的數(shù)據(jù)分析。二者相比,前者更適合于一些經(jīng)驗較少的生物研究人員或非生物信息學(xué)研究者進(jìn)行通用常規(guī)的分析,后者更適合于具備豐富的數(shù)據(jù)分析經(jīng)驗且需進(jìn)行特定參數(shù)分析的研究人員。

3 實驗結(jié)果與性能分析

為多方面驗證DEA-container性能,本研究采用蒙古羊亞系小尾寒羊的高通量測序數(shù)據(jù),探討FSHR基因在不同繁殖力的小尾寒羊中的差異性表達(dá)。我們分別選取SRR5786000、SRR5786-001、SRR5786002、SRR5786003 雙端數(shù)據(jù)(共計59.4 GB)作為樣本數(shù)據(jù)集進(jìn)行分析。樣本數(shù)據(jù)集的詳細(xì)信息如表3所示。

表3 樣本數(shù)據(jù)集Table 3 List of test samples

3.1 準(zhǔn)確性驗證

準(zhǔn)確性是評估一切方法框架最基本的前提。為驗證DEA-container框架分析結(jié)果的準(zhǔn)確性,首先采用兩種方式對同一原始數(shù)據(jù)進(jìn)行處理分析,即傳統(tǒng)人工方式和DEA-container框架全自動方式;隨后將兩種方式的結(jié)果數(shù)據(jù)進(jìn)行收集與對比,通過驗證結(jié)果文件的一致性與完整性來確保數(shù)據(jù)的準(zhǔn)確性。文中主要選取的核心輸出文件有TopHat2 結(jié)果文件:accepted_hits.bam、Cuffmerge結(jié)果文件:merged.gtf、Cuffnorm結(jié)果文件:genes.fpkm_table、Cuffdiff結(jié)果文件:isoforms.fpkm_tracking。驗證方法如圖4所示,結(jié)果表明不同處理方式所產(chǎn)生的結(jié)果文件內(nèi)容完全一致,證明DEA-container框架具有準(zhǔn)確性。

圖4 準(zhǔn)確性驗證方法Fig.4 The method for verifying accuracy

3.2 高效性驗證

為測試DEA-container的處理效率,本研究一方面在實驗室的惠普工作站(16 GB RAM,4 CPU Intel Xeon)實驗平臺對小尾寒羊中目標(biāo)基因的差異性表達(dá)分析任務(wù)通過Web端進(jìn)行6次提交和執(zhí)行;另一方面,采用傳統(tǒng)人工方式對相同數(shù)據(jù)集進(jìn)行人工的分析處理。我們分別記錄了不同處理方式下每個完整分析周期的總運行時間并求取平均值,結(jié)果如表4、圖5所示。

圖5 高效性對比驗證Fig.5 The comparison of the efficiency

表4 處理時間統(tǒng)計表Table 4 Statistics of the processing time

在傳統(tǒng)人工方式處理過程中,TopHat2耗時約 26 h,Cufflinks耗時約 9.9 h,Cuffnorm、Cuffquant、Cuffdiff總耗時約7.4 h,一個完整周期約耗時43 h 18 min。其處理效率差,耗時長,需研究人員手動部署和執(zhí)行各個處理步驟,工作繁瑣且難度大易出錯,在大數(shù)據(jù)時代下該方式不可取[13]。文中采用DEA-container框架全自動化方式進(jìn)行的6次重復(fù)測試均實現(xiàn)了轉(zhuǎn)錄組差異性表達(dá)分析全流程的自動化執(zhí)行,對6次執(zhí)行效率進(jìn)行計算后取平均值,得到一次完整的分析周期約為12 h 6 min。從以上分析可知,DEA-container的處理效率與傳統(tǒng)人工方式相比提升了2倍多,分析時間縮短約72%,很大程度提升了分析效率,在有效降低成本的同時保證了數(shù)據(jù)分析的高效性。

3.3 部署時間對比驗證

DEA-container的突出優(yōu)勢為一次封裝、隨處運行、可跨平臺、不依賴任何編程語言。為了驗證這些優(yōu)勢,本研究采用兩種方式,即DEA-container一鍵化腳本部署方式和傳統(tǒng)人工方式,將轉(zhuǎn)錄組差異性表達(dá)分析最佳實踐流程部署至11臺主機(jī)上。但由于傳統(tǒng)人工方式基于Linux的命令行,使大量的參數(shù)配置以及復(fù)雜的命令行操作均受到有無Linux基礎(chǔ)的限制,導(dǎo)致部署時間受到較大的影響,因此我們分別選取Linux技能熟練程度不同的人員進(jìn)行操作,對兩種方式的部署時間進(jìn)行記錄,最后求取平均值,結(jié)果如圖6所示。

由圖6可知,傳統(tǒng)人工方式的部署時間受Linux技能熟練程度影響較大,其完整部署時間平均約為26.09 h,大量時間耗費在手動的環(huán)境配置與軟件的部署。與傳統(tǒng)人工方式不同,DEA-contianer一鍵化腳本部署與研究人員Linux技能熟練程度無關(guān),該方式將手動操作的全過程封裝且自動化,完整的腳本執(zhí)行時間平均約26.7 min,高效而自動化,完全解決了軟件的環(huán)境配置難、安裝難等問題,成功地避免了研究人員重復(fù)繁瑣的手動操作,很大程度上減少了用戶工作量、降低了用戶工作難度,有效提升了分析效率。同時,DEA-container一鍵化腳本部署方式降低了對生物研究人員計算機(jī)專業(yè)技能的要求,更為簡單易用。

圖6 部署時間對比驗證Fig.6 The time comparison for deployment

3.4 實驗小結(jié)

兩種處理方式的比對結(jié)果如表5所示,與傳統(tǒng)人工方式相比,DEA-container方式具有如下優(yōu)勢:1)分析效率顯著提升,是傳統(tǒng)人工方式的2倍多,分析時間縮短約72%;2)具備較好的移植性,可跨平臺隨處運行;3)環(huán)境部署快速便捷,僅需執(zhí)行一鍵化部署腳本即可;4)易擴(kuò)展,可按需集成多種分析流程與工具軟件;5)將復(fù)雜流程全面封裝,極大簡化了用戶操作。

4 結(jié)語

本文針對轉(zhuǎn)錄組差異性表達(dá)分析復(fù)雜流程形成了一個簡單易用的框架,用戶只需要通過Web界面對目標(biāo)分析序列進(jìn)行提交,無需做其他額外操作,而且后臺處理分析完畢后會直接給出分析結(jié)果供用戶下載,極大地簡化了操作流程,提高了分析效率,使得研究人員可以更好地進(jìn)行數(shù)據(jù)分析工作,從海量數(shù)據(jù)中得到更可靠并具有生物學(xué)意義的結(jié)果。

表5 全局對比Table 5 The global contrast

目前,我們已將轉(zhuǎn)錄組差異性表達(dá)分析容器及相關(guān)內(nèi)容公布發(fā)放至GitHub(https://github.com/fullblossom/the-Differential-expression-analysiscontainer),可供用戶下載使用。在后續(xù)的工作中我們將按需求形成多種分析容器,并將其公布至平臺供用戶下載和使用,以期為生物學(xué)研究者提供簡單易用且高效的分析工具。