基于16S rRNA序列物種鑒定的改進(jìn)向量空間模型算法①

祝 斌,亓合媛,馬俊才1,

1(中國(guó)科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

2(中國(guó)科學(xué)院大學(xué),北京 100049)

3(中國(guó)科學(xué)院 微生物研究所,北京 100101)

在過去的幾十年中,隨著生物學(xué)數(shù)據(jù)的大量累積,以及計(jì)算機(jī)技術(shù)、數(shù)學(xué)和生物學(xué)交叉學(xué)科的崛起,21世紀(jì)已經(jīng)進(jìn)入了云計(jì)算和大數(shù)據(jù)時(shí)代.云計(jì)算時(shí)代也為基因序列比對(duì)能夠在較短時(shí)間完成提供了堅(jiān)實(shí)的基礎(chǔ).物種鑒定是用來描述物種間近緣關(guān)系和進(jìn)化層次的非常有用的一種工具.最初,物種鑒定常常基于單個(gè)基因序列或是很少的幾個(gè)基因序列進(jìn)行對(duì)比,這種方法雖然簡(jiǎn)單易行,但是由于橫向基因轉(zhuǎn)移(Horizontal Gene Transfer,HGT)、并系同源基因(Paralog)以及物種進(jìn)化差異等因素的出現(xiàn),這種方法受到了質(zhì)疑.

基因測(cè)序方法在一定程度上解決了單基因序列比對(duì)出現(xiàn)的問題,保證了系統(tǒng)發(fā)育樹的合理性.但是,與此同時(shí),隨著序列數(shù)據(jù)的增加,計(jì)算時(shí)間呈指數(shù)式增長(zhǎng)(如圖1所示),因此,計(jì)算效率成為了亟待解決的問題.

圖1 Genbank數(shù)據(jù)庫(kù)序列數(shù)目隨年份走勢(shì)圖

近期,關(guān)于物種鑒定的方法逐步地出現(xiàn)了全基因組序列對(duì)比.在此期間,眾多學(xué)者提出:比較兩個(gè)完整的基因組意義并不大,原因在于:每個(gè)物種都有自己特定的基因含量和基因順序,此外基因組的數(shù)量是不同的;另外,微生物的基因數(shù)據(jù)也是需要人工處理操作的,不同的實(shí)驗(yàn)室處理數(shù)據(jù)的不同會(huì)造成結(jié)果在一定程度上的差異,進(jìn)而使得結(jié)果不具備更完整的說服力,非序列比對(duì)方法便應(yīng)運(yùn)而生.

非序列比對(duì)方法在計(jì)算效率上明顯優(yōu)于前者,操作也較簡(jiǎn)單,運(yùn)算效率較高.近幾年,關(guān)于非序列對(duì)比的方法也不斷更新.目前比較常見的方法有:K串組份向量方法,(0,l)序列法,DNA Walk,壓縮矩陣法,表示法,CGR方法,Nandy二維圖形[1]等.近年的新興算法在物種鑒定方面上的應(yīng)用上不夠廣泛,在所有的非序列比對(duì)算法中,使用最為廣泛且傳統(tǒng)的算法為基于TFIDF檢索技術(shù)的向量空間模型(Vector Space Model,VSM)算法,然而其物種鑒定的分類效果得不到保證[2],原因在于該算法沒有借鑒到微生物的背景,因此無法消除在基因突變和物種進(jìn)化的背景下,基因序列的噪音影響.因此,在確定非序列對(duì)比算法具備了提高運(yùn)算效率的優(yōu)點(diǎn),以及向量空間模型算法在眾多經(jīng)典和最新的文獻(xiàn)[3]中使用較為廣泛的特點(diǎn)之后,為此,本文以如何改進(jìn)向量空間模型算法,進(jìn)一步達(dá)到提高運(yùn)算效率和保證分類效果質(zhì)量?jī)煞矫鏋橹饕康?

在眾多生物系統(tǒng)發(fā)育相關(guān)性水平指標(biāo)中,16S rRNA基因序列具有如下特征:

1)普遍存在于一切細(xì)胞內(nèi);

2)機(jī)體生理功能穩(wěn)定且重要;

3)在微生物中含量高,且容易提取;

4)編碼基因比較穩(wěn)定;

5)序列相對(duì)保守;

6)相對(duì)分子量適中;

7)基因序列長(zhǎng)度適中;

8)既含有高度保守的序列區(qū)域,又含有高度變化的序列區(qū)域.

基于以上的各個(gè)特點(diǎn),16S rRNA 基因序列具備最佳的鑒定特征,是本文改進(jìn)向量空間模型算法的應(yīng)用數(shù)據(jù),可以為物種鑒定打下堅(jiān)實(shí)的基礎(chǔ).

綜上,本文以16S rRNA 基因序列為應(yīng)用對(duì)象,使用改進(jìn)向量空間模型算法為核心,以達(dá)到快速分類和保證分類質(zhì)量的研究目的.

1 背景及相關(guān)工作

分子生物系統(tǒng)發(fā)展史的出現(xiàn)以及基因測(cè)序方面的進(jìn)步,大大加深了人們對(duì)物種進(jìn)化的理解.因此,物種分類和鑒定在分子水平上的進(jìn)步已經(jīng)為微生物的分類提供了一個(gè)具有實(shí)用價(jià)值的工具.

目前分子系統(tǒng)發(fā)展史有兩大重要研究成果:一是線粒體和葉綠體之間具有內(nèi)生共體特性,二是目前為止,生物可劃分為古生菌,細(xì)菌,和真核生物三個(gè)生物領(lǐng)域.然而,隨著完整的微生物基因組數(shù)據(jù)的逐步添加,實(shí)驗(yàn)結(jié)果逐漸地對(duì)公眾預(yù)期提出了質(zhì)疑[4],在這一爭(zhēng)議過程中,仍有幾個(gè)實(shí)驗(yàn)試圖從完整的基因組中推斷出原核生物發(fā)展史.以上實(shí)驗(yàn)使用的方法包括利用基因含量[5],直系同源基因簇的存在/缺失值比例[6],父系樹[7],保存基因?qū)8]等方法.然而這些方法最終都依賴于序列比對(duì)這一傳統(tǒng)思路,到目前為止,還沒有一種能夠被廣泛接受且用于從完整基因組數(shù)據(jù)中推斷出系統(tǒng)發(fā)育樹的方法.

此后,逐漸出現(xiàn)了非序列比對(duì)的方法[9],計(jì)算效率和結(jié)果都得到了廣泛的認(rèn)可,因此成為了除BLAST算法以外物種分類與鑒定方面不可或缺的方法.而向量空間模型(VSM)算法在眾多前沿文獻(xiàn)中使用的頻率較高,由此可見,目前向量空間模型算法是非序列比對(duì)算法中構(gòu)建系統(tǒng)發(fā)育樹的主流算法.因此,對(duì)其算法的改進(jìn)具有重大意義.

根據(jù)相關(guān)文獻(xiàn)[10]的說明,截至目前,使用16S rRNA基因序列對(duì)物種進(jìn)行鑒定和分類的項(xiàng)目有:美國(guó)的Greengenes,RDP核糖體數(shù)據(jù)庫(kù),以及韓國(guó)的EzTaxon.以上項(xiàng)目的核心基礎(chǔ)仍是利用BLAST局部比對(duì)算法進(jìn)行快速分類,輸出初始排名結(jié)果,隨后使用雙序列全局比對(duì),給出在參考樣本數(shù)據(jù)庫(kù)中與待測(cè)序列最為接近的排名序列,以此作為參考,對(duì)樣本序列進(jìn)行鑒定和分類.

根據(jù)前面的分析,我們發(fā)現(xiàn),用于物種鑒定的主流算法仍是基于BLAST的序列比對(duì)算法,然而由于該算法出現(xiàn)計(jì)算量過于龐大,運(yùn)算效率低以及資源消耗較高等問題,使用VSM方法能夠有效地解決上述問題.

VSM算法的運(yùn)算效率相比于BLAST算法更優(yōu),此特點(diǎn)解決了BLAST算法的核心問題,但該算法的不足之處在于其分類效果遠(yuǎn)遠(yuǎn)沒有主流BLAST鑒定算法更為優(yōu)越.因此,對(duì)VSM算法的改進(jìn)就具有了現(xiàn)實(shí)意義,而改進(jìn)的VSM算法可以作為物種鑒定的另一種有效工具方便科研人員參考和使用.

此外,經(jīng)典文獻(xiàn)[11] 提到的K-String組份向量算法在病毒[12],原核生物[13–17],真菌[18],葉綠體序列[19]以及人體的腸道元基因組[20]有了成功的應(yīng)用.

綜上所述,本文旨在對(duì)常用的VSM算法進(jìn)行改進(jìn),將該改進(jìn)VSM算法應(yīng)用于基于16S rRNA序列的物種鑒定領(lǐng)域,達(dá)到運(yùn)算效率和分類質(zhì)量?jī)煞矫娴奶岣咝Ч?本文后續(xù)的內(nèi)容邏輯為:在第2節(jié)介紹兩種VSM模型算法以及兩種算法的區(qū)別,一種是基于TFIDF檢索技術(shù)的VSM模型算法,另一種是借鑒經(jīng)典文獻(xiàn)[1]后的改進(jìn)VSM模型算法.此外,本文還給出了改進(jìn)VSM模型算法中遺傳距離在巴拿赫空間下的等價(jià)替代公式,并給出了相關(guān)說明;同時(shí),第2節(jié)給出本文為測(cè)試改進(jìn)VSM算法運(yùn)算效率和分類排名質(zhì)量?jī)煞矫嫘Ч褂玫臄?shù)據(jù)集來源,以及對(duì)應(yīng)的運(yùn)算時(shí)間和排名效果結(jié)果匯總及相應(yīng)分析;第3節(jié)是對(duì)接下來研究工作的討論與未來展望.

2 VSM算法與改進(jìn)VSM算法

2.1 VSM算法原理[21,22]介紹

本文將以16S rRNA基因序列分析為研究背景,介紹VSM在該背景下的操作流程.

一個(gè)物種16S rRNA基因序列文本,其堿基只有AGCT四種,將堿基序列劃分為不同的K子串,那么此排列方式就有4K種可能,通過計(jì)算詞頻和逆文檔頻率,最終得到該16S rRNA序列文本對(duì)應(yīng)的權(quán)重向量,維數(shù)為 1×4K.

圖2 序列相似度

圖2中,每一項(xiàng)的權(quán)重都由詞頻和逆文檔頻率綜合表示.

假設(shè)有N條樣本序列,記為D={d1,d2,···,dN},詞頻fij表示K串詞項(xiàng)wi在序列dj中出現(xiàn)的次數(shù),ni為含文本w的數(shù)量,逆文本頻率計(jì)算公式為:

其中,序列dj中詞項(xiàng)wi的TF-IDF權(quán)重公式為:

最后,樣本序列相似度量值計(jì)算公式為:

2.2 改進(jìn)VSM算法原理介紹

該算法涉及6個(gè)步驟,將分別作出說明.

第一步.計(jì)算K串詞項(xiàng)出現(xiàn)的頻率.

長(zhǎng)度為L(zhǎng)的16S rRNA序列 α1α2···αL,選取長(zhǎng)度為K(K

第二步.計(jì)算隨機(jī)突變背景下的噪音頻率.

隨機(jī)突變?cè)诜肿铀缴匣蚨嗷蛏僖噪S機(jī)的方式發(fā)生,而基因重組的選擇決定了進(jìn)化的方向.以上因素導(dǎo)致了一些K串詞項(xiàng)產(chǎn)成了一定的隨機(jī)性.為了還原該K串詞項(xiàng)的原始頻率,本文需要對(duì)此噪音進(jìn)行刻畫,根據(jù)最大熵原理的推導(dǎo)過程,我們得到了噪音頻率公式:

第三步.計(jì)算修正后K串詞項(xiàng)頻率.

綜合前兩步,本文給出了修正后頻率計(jì)算公式:

第四步.計(jì)算每一個(gè)16S rRNA序列修正后的特征向量.

將每一個(gè)可能的序列子串 α1α2···αK的頻率作為一個(gè)物種的特征向量的元素.為了進(jìn)一步簡(jiǎn)化這一個(gè)定義,我們定義ai為所有排列好的K子串中第i種子串類型對(duì)應(yīng)特征向量中的第i個(gè)分量.這里i從1到4K循環(huán).因此,我們可以得出對(duì)于16S rRNA序列A的特征向量:

以此類推,對(duì)于物種B,我們?nèi)杂刑卣飨蛄?

第五步.計(jì)算各序列間的遺傳距離.

這里同樣以序列A,B為例,兩序列間的遺傳距離使用傳統(tǒng)的夾角余弦進(jìn)行表示,公式如下:

由公式(9)知,夾角余弦數(shù)值的取值范圍為[–1,1].若將夾角余弦記做物種間的遺傳距離,則有:兩物種特征向量對(duì)應(yīng)的遺傳距離越大,說明兩個(gè)物種之間的相關(guān)性越強(qiáng);反之,遺傳距離越小,說明兩物種之間的相關(guān)性越弱.為了符合直觀,表達(dá)相關(guān)性強(qiáng),對(duì)應(yīng)遺傳距離小;相關(guān)性弱,則遺傳距離大的說法,本文對(duì)此距離公式進(jìn)行標(biāo)準(zhǔn)化修正,公式為:

第六步.計(jì)算待測(cè)樣本與參考序列庫(kù)之間的遺傳距離,從小到大進(jìn)行排序,輸出前十名相關(guān)性最強(qiáng)的序列及其遺傳信息,以輔助科研人員參考和進(jìn)行物種分類和鑒定工作.

2.3 巴拿赫空間下等價(jià)替換的遺傳距離公式

遺傳距離的定義是計(jì)算分子生物學(xué)中一個(gè)重要環(huán)節(jié).該距離的定義需要滿足以下三個(gè)條件(記D(x,y)為兩個(gè)物種間的遺傳距離):

非負(fù)性:D(x,y)≥0,D(x,y)=0等價(jià)于x=y;

對(duì)稱性:D(y,x)=D(x,y);

三角形不等式:任意三個(gè)物種z,x,y,距離恒滿足:D(z,y)+D(y,x)≥D(x,y).

顯然,在本文的第2.2節(jié)中的公式(10)符合遺傳距離的定義.這里值得一提的是,夾角余弦公式使用的是內(nèi)積空間下的2-范數(shù).根據(jù)向量范數(shù)的等價(jià)性定理:

設(shè)||x||s,||x||t為Rn上向量的任意兩種范數(shù),則存在常數(shù)c1,c2＞0,使得對(duì)一切x∈Rn,有c1||x||s≤||x||t≤c2||x||s.

以及極化恒等式:實(shí)線性空間上的內(nèi)積和范數(shù)有以下關(guān)系:

綜合上述內(nèi)容,本文將公式(8)中內(nèi)積范數(shù)進(jìn)行重新定義,給出在1-范數(shù)和無窮范數(shù)下的計(jì)算公式,公式如下:

2.4 改進(jìn)VSM模型算法[23,24]運(yùn)算效率和排名結(jié)果分析

2.4.1 待測(cè)樣本與測(cè)試數(shù)據(jù)集

本文所使用的16S rRNA樣本序列數(shù)據(jù),來源于(863計(jì)劃,課題編號(hào):2014AA021501)中通過質(zhì)檢工具pipeline篩選整理出的高質(zhì)量16S rRNA基因序列參考數(shù)據(jù)庫(kù).

這里,本文從參考數(shù)據(jù)庫(kù)中隨機(jī)選取了8000條樣本序列.其中,將前6000條為參考序列樣本庫(kù),剩余的2000條作為待測(cè)樣本進(jìn)行測(cè)試.

為了簡(jiǎn)化名稱,這里依次定義序列編號(hào)為G1,G2,…,G6000,G6001,…,G8000.其中,G1,G2,…,G6000為參考數(shù)據(jù)庫(kù),G6001,…,G8000為待測(cè)樣本.

2.4.2 改進(jìn)VSM算法運(yùn)算效率與blast運(yùn)算效率結(jié)果

結(jié)合第2.2節(jié)所述,可以發(fā)現(xiàn),本文所選取的6000條參考樣本序列文本可以通過改進(jìn)向量空間模型算法進(jìn)行計(jì)算,得出對(duì)應(yīng)的6000個(gè)特征向量,是本實(shí)驗(yàn)的預(yù)處理階段.因此,以上6000個(gè)特征向量的運(yùn)算時(shí)間完全不需要計(jì)入該算法的運(yùn)算時(shí)間,這也是該算法提高運(yùn)算效率的一大優(yōu)勢(shì).

這里,本文首先按照第2.2節(jié)所述的操作步驟逐一進(jìn)行:(這里以K=4為例)

第一步.對(duì)前6000條16S rRNA參考樣本序列G1,G2,…,G6000,逐一帶入公式(6),計(jì)算出每一個(gè)序列文本對(duì)應(yīng)的修正頻率特征向量A1,A2,…,A6000.

說明.此階段為數(shù)據(jù)預(yù)處理階段,不占用算法計(jì)算時(shí)間;其中每一個(gè)特征向量的維數(shù)為:1×44即1×256.

第二步.i=6001,對(duì)待測(cè)樣本Gi計(jì)算出對(duì)應(yīng)的修正頻率特征向量Ai.

第三步.計(jì)算Gi與G1,G2,…,G6000序列之間的遺傳距離C(Gi,G1),C(Gi,G2),…,C(Gi,G6000),依次記為遺傳距離d1,d2,…,d6000.

第四步.對(duì)上述6000個(gè)遺傳距離d1,d2,…,d6000按照遞增的順序進(jìn)行排序,輸出相關(guān)性較高的前十條序列排名結(jié)果作為物種鑒定的參考初步排名結(jié)果.

第五步:i=i+1,直至計(jì)算至最后一個(gè)待測(cè)樣本G8000.

說明1.其中第2,3,4步為一個(gè)待測(cè)樣本與6000條參考樣本序列G1,G2,…,G6000的整個(gè)運(yùn)算過程,其花費(fèi)時(shí)間也為該改進(jìn)向量空間模型算法的單個(gè)樣本進(jìn)行排名輸出的運(yùn)算時(shí)間.

說明2.本文將i從6001依次逐個(gè)循環(huán)至8000進(jìn)行操作運(yùn)算,消耗的總時(shí)間除以2000,記為改進(jìn)VSM模型算法的測(cè)試時(shí)間.

緊接著,本文使用BLAST本地構(gòu)建6000條參考數(shù)據(jù)庫(kù),使用blastn程序?qū)?000條待測(cè)樣本進(jìn)行逐一運(yùn)算,輸出結(jié)果,其花費(fèi)的時(shí)間除以2000,同樣記為blast算法測(cè)試時(shí)間.

說明1.這里使用blastn命令:

blastn -query 6001.fa -db Sequence6000 -evalue 1e-5 -out blast6001.xls -outfmt 6 -num_alignments 10 -num_threads 1

說明2.以上參數(shù)中,-query6001.fa為待測(cè)樣本G6001的fa格式文件,-db Sequence6000表示6000條參考樣本的本地化數(shù)據(jù)庫(kù),-evalue 1e–5表示控制誤差,-outfmt 6表示輸出文件排版格式按照格式6進(jìn)行輸出,-num_alignments 10表示輸出排名前10的序列結(jié)果,-num_threads 1表示單線程.

說明3.本文使用改進(jìn)VSM算法,使用的是c程序,改進(jìn)VSM算法和blastn算法均在Ubuntu 12.04.4 LTS同一個(gè)操作環(huán)境下運(yùn)行.

最后,綜合上述兩項(xiàng)內(nèi)容的操作,本文給出了改進(jìn)VSM模型算法和BLAST算法運(yùn)行時(shí)間.

表1 改進(jìn)VSM算法與BLAST算法運(yùn)行效率(單位:ms)

2.4.3 改進(jìn)VSM算法運(yùn)算效率與BLAST算法排名結(jié)果

本文選擇輸出前10名用于比較兩種算法的鑒定效果,原因在于:本文使用的參考數(shù)據(jù)集序列數(shù)量為隨機(jī)抽樣后的6000條序列,序列數(shù)量相對(duì)較小;且在物種鑒定領(lǐng)域中,一般輸出BLAST相似度98%以上的排名結(jié)果,這里使用BLAST輸出序列相似度98%以上的序列數(shù)均小于10條,因此選擇前10名作為評(píng)價(jià)的參考標(biāo)準(zhǔn).

按照第2.2節(jié)的操作進(jìn)行,本文得出了對(duì)應(yīng)的排名結(jié)果,這里以待測(cè)樣本G6001的排名結(jié)果為例,如表2所示.

表2 改進(jìn)VSM算法與BLAST算法排名結(jié)果

本文K=8時(shí),將改進(jìn)VSM算法輸出排名與blast排名結(jié)果重復(fù)率進(jìn)行統(tǒng)計(jì),最終得出:所有的2000個(gè)待測(cè)16S rRNA樣本序列,通過使用改進(jìn)VSM算法輸出的前十名排名結(jié)果,其檢出率已達(dá)到98.0%.

此外,若將輸出前10名序列信息,改為輸出前50,或前100名,我們發(fā)現(xiàn)檢出率和K相關(guān),隨著K越大,算法的檢出率相對(duì)會(huì)越優(yōu);且當(dāng)K=10,輸出前100名序列信息時(shí),檢出率達(dá)到了97.6%,證明了該算法的收斂性.

2.4.4 改進(jìn)VSM算法與blast算法對(duì)比綜合分析

根據(jù)表2,我們可以看到:隨著K由4逐步遞增至K=8,其輸出的排名結(jié)果檢出率由30%上升至90%;此外,G6000的排名也逐漸靠前,以及排名第一的G5997和BLAST的第1名結(jié)果吻合.

根據(jù)表1,我們可以看出,隨著K的增大,運(yùn)算時(shí)間也成約4倍遞增,然而當(dāng)K=8時(shí),BLAST運(yùn)算時(shí)間約為改進(jìn)VSM模型算法的50倍.

綜合以上運(yùn)算效率和排名結(jié)果兩方面的分析,我們可以得出改進(jìn)VSM算法維持了其計(jì)算效率的優(yōu)越性,并改進(jìn)了排名結(jié)果,提高了檢出率.

3 展望

本文提出的改進(jìn)VSM算法,是將經(jīng)典文獻(xiàn)中的K串組份向量空間模型算法應(yīng)用于微生物16S rRNA序列的物種鑒定中,并對(duì)遺傳距離公式進(jìn)行改進(jìn),以期克服傳統(tǒng)VSM模型算法在物種鑒定方面上的不足,進(jìn)一步提高物種鑒定的檢出率,最終保證物種鑒定的質(zhì)量效果.

后續(xù)的研究工作還包括:改進(jìn)VSM模型算法多線程模板設(shè)置,進(jìn)一步提升該算法的運(yùn)算效率.初步設(shè)置思路為:將參考數(shù)據(jù)集劃分成多個(gè)模塊,然后將待測(cè)樣本分別與各個(gè)模塊進(jìn)行比對(duì),輸出各自的遺傳距離向量,接著將各個(gè)向量匯集成一個(gè)完整的向量,最終對(duì)該向量進(jìn)行排序輸出最終結(jié)果.