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

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

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

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

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

1 概念定義與原理

1.1 Docker容器技術(shù)

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

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

1.2 最佳實(shí)踐流程

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

2016年,Ghosh和Cean針對(duì)轉(zhuǎn)錄組差異性表達(dá)分析提出了一套最佳實(shí)踐流程[7],該流程在提出的幾年里被生物信息學(xué)領(lǐng)域?qū)＜宜J(rèn)可并延續(xù)至今仍被廣泛使用。因此,文中將該最佳實(shí)踐流程Bowtie2[8]-TopHat2[9]-Cufflinks[10]-cummeRbund[11]與數(shù)據(jù)預(yù)處理環(huán)節(jié)相整合,在已構(gòu)建好參考基因庫(kù)的前提下,將轉(zhuǎn)錄組差異性表達(dá)分析主要分為8個(gè)處理步驟,形成圖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格式的原始測(cè)序數(shù)據(jù)轉(zhuǎn)換為fastq格式;其次,通過(guò)FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc)對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量檢測(cè);最后,通過(guò)Trimmomatic質(zhì)量檢測(cè)[12]將質(zhì)量低或帶有接頭(a-dapters)的測(cè)序數(shù)據(jù)進(jìn)行修剪和過(guò)濾,進(jìn)而得到高質(zhì)量測(cè)序數(shù)據(jù)。

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

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

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

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

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

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

8)結(jié)果可視化:采用cummeRbund軟件對(duì)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)?；A(chǔ)層除硬件設(shè)施與操作系統(tǒng)外,還需為上層提供運(yùn)行環(huán)境的Docker引擎;核心層為轉(zhuǎn)錄組差異性表達(dá)分析容器,屬于核心功能組件,可按需集成不同的復(fù)雜分析流程;調(diào)度層即通過(guò)多腳本嵌套觸發(fā)對(duì)服務(wù)、資源的協(xié)調(diào)和控制;應(yīng)用層即為用戶接口層。該框架各層次之間協(xié)同配合,進(jìn)而實(shí)現(xiàn)高效自動(dòng)化的轉(zhuǎn)錄組差異性表達(dá)分析。

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

2.2 核心層

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

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

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

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

2.3 多腳本聯(lián)動(dòng)過(guò)程

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與性能分析

為多方面驗(yàn)證DEA-container性能,本研究采用蒙古羊亞系小尾寒羊的高通量測(cè)序數(shù)據(jù),探討FSHR基因在不同繁殖力的小尾寒羊中的差異性表達(dá)。我們分別選取SRR5786000、SRR5786-001、SRR5786002、SRR5786003 雙端數(shù)據(jù)(共計(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)確性驗(yàn)證

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

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

3.2 高效性驗(yàn)證

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

圖5 高效性對(duì)比驗(yàn)證Fig.5 The comparison of the efficiency

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

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

3.3 部署時(shí)間對(duì)比驗(yàn)證

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

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

圖6 部署時(shí)間對(duì)比驗(yàn)證Fig.6 The time comparison for deployment

3.4 實(shí)驗(yàn)小結(jié)

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

4 結(jié)語(yǔ)

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

表5 全局對(duì)比Table 5 The global contrast

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