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

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

1(中國科學院 計算機網絡信息中心,北京 100190)

2(中國科學院大學,北京 100049)

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

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

基因測序方法在一定程度上解決了單基因序列比對出現的問題,保證了系統發育樹的合理性.但是,與此同時,隨著序列數據的增加,計算時間呈指數式增長(如圖1所示),因此,計算效率成為了亟待解決的問題.

圖1 Genbank數據庫序列數目隨年份走勢圖

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

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

在眾多生物系統發育相關性水平指標中,16S rRNA基因序列具有如下特征:

1)普遍存在于一切細胞內;

2)機體生理功能穩定且重要;

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

4)編碼基因比較穩定;

5)序列相對保守;

6)相對分子量適中;

7)基因序列長度適中;

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

基于以上的各個特點,16S rRNA 基因序列具備最佳的鑒定特征,是本文改進向量空間模型算法的應用數據,可以為物種鑒定打下堅實的基礎.

綜上,本文以16S rRNA 基因序列為應用對象,使用改進向量空間模型算法為核心,以達到快速分類和保證分類質量的研究目的.

1 背景及相關工作

分子生物系統發展史的出現以及基因測序方面的進步,大大加深了人們對物種進化的理解.因此,物種分類和鑒定在分子水平上的進步已經為微生物的分類提供了一個具有實用價值的工具.

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

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

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

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

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

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

綜上所述,本文旨在對常用的VSM算法進行改進,將該改進VSM算法應用于基于16S rRNA序列的物種鑒定領域,達到運算效率和分類質量兩方面的提高效果.本文后續的內容邏輯為:在第2節介紹兩種VSM模型算法以及兩種算法的區別,一種是基于TFIDF檢索技術的VSM模型算法,另一種是借鑒經典文獻[1]后的改進VSM模型算法.此外,本文還給出了改進VSM模型算法中遺傳距離在巴拿赫空間下的等價替代公式,并給出了相關說明;同時,第2節給出本文為測試改進VSM算法運算效率和分類排名質量兩方面效果所使用的數據集來源,以及對應的運算時間和排名效果結果匯總及相應分析;第3節是對接下來研究工作的討論與未來展望.

2 VSM算法與改進VSM算法

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

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

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

圖2 序列相似度

圖2中,每一項的權重都由詞頻和逆文檔頻率綜合表示.

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

其中,序列dj中詞項wi的TF-IDF權重公式為:

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

2.2 改進VSM算法原理介紹

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

第一步.計算K串詞項出現的頻率.

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

第二步.計算隨機突變背景下的噪音頻率.

隨機突變在分子水平上或多或少以隨機的方式發生,而基因重組的選擇決定了進化的方向.以上因素導致了一些K串詞項產成了一定的隨機性.為了還原該K串詞項的原始頻率,本文需要對此噪音進行刻畫,根據最大熵原理的推導過程,我們得到了噪音頻率公式:

第三步.計算修正后K串詞項頻率.

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

第四步.計算每一個16S rRNA序列修正后的特征向量.

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

以此類推,對于物種B,我們仍有特征向量:

第五步.計算各序列間的遺傳距離.

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

由公式(9)知,夾角余弦數值的取值范圍為[–1,1].若將夾角余弦記做物種間的遺傳距離,則有:兩物種特征向量對應的遺傳距離越大,說明兩個物種之間的相關性越強;反之,遺傳距離越小,說明兩物種之間的相關性越弱.為了符合直觀,表達相關性強,對應遺傳距離小;相關性弱,則遺傳距離大的說法,本文對此距離公式進行標準化修正,公式為:

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

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

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

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

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

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

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

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

以及極化恒等式:實線性空間上的內積和范數有以下關系:

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

2.4 改進VSM模型算法[23,24]運算效率和排名結果分析

2.4.1 待測樣本與測試數據集

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

這里,本文從參考數據庫中隨機選取了8000條樣本序列.其中,將前6000條為參考序列樣本庫,剩余的2000條作為待測樣本進行測試.

為了簡化名稱,這里依次定義序列編號為G1,G2,…,G6000,G6001,…,G8000.其中,G1,G2,…,G6000為參考數據庫,G6001,…,G8000為待測樣本.

2.4.2 改進VSM算法運算效率與blast運算效率結果

結合第2.2節所述,可以發現,本文所選取的6000條參考樣本序列文本可以通過改進向量空間模型算法進行計算,得出對應的6000個特征向量,是本實驗的預處理階段.因此,以上6000個特征向量的運算時間完全不需要計入該算法的運算時間,這也是該算法提高運算效率的一大優勢.

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

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

說明.此階段為數據預處理階段,不占用算法計算時間;其中每一個特征向量的維數為:1×44即1×256.

第二步.i=6001,對待測樣本Gi計算出對應的修正頻率特征向量Ai.

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

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

第五步:i=i+1,直至計算至最后一個待測樣本G8000.

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

說明2.本文將i從6001依次逐個循環至8000進行操作運算,消耗的總時間除以2000,記為改進VSM模型算法的測試時間.

緊接著,本文使用BLAST本地構建6000條參考數據庫,使用blastn程序對2000條待測樣本進行逐一運算,輸出結果,其花費的時間除以2000,同樣記為blast算法測試時間.

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

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

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

說明3.本文使用改進VSM算法,使用的是c程序,改進VSM算法和blastn算法均在Ubuntu 12.04.4 LTS同一個操作環境下運行.

最后,綜合上述兩項內容的操作,本文給出了改進VSM模型算法和BLAST算法運行時間.

表1 改進VSM算法與BLAST算法運行效率(單位:ms)

2.4.3 改進VSM算法運算效率與BLAST算法排名結果

本文選擇輸出前10名用于比較兩種算法的鑒定效果,原因在于:本文使用的參考數據集序列數量為隨機抽樣后的6000條序列,序列數量相對較小;且在物種鑒定領域中,一般輸出BLAST相似度98%以上的排名結果,這里使用BLAST輸出序列相似度98%以上的序列數均小于10條,因此選擇前10名作為評價的參考標準.

按照第2.2節的操作進行,本文得出了對應的排名結果,這里以待測樣本G6001的排名結果為例,如表2所示.

表2 改進VSM算法與BLAST算法排名結果

本文K=8時,將改進VSM算法輸出排名與blast排名結果重復率進行統計,最終得出:所有的2000個待測16S rRNA樣本序列,通過使用改進VSM算法輸出的前十名排名結果,其檢出率已達到98.0%.

此外,若將輸出前10名序列信息,改為輸出前50,或前100名,我們發現檢出率和K相關,隨著K越大,算法的檢出率相對會越優;且當K=10,輸出前100名序列信息時,檢出率達到了97.6%,證明了該算法的收斂性.

2.4.4 改進VSM算法與blast算法對比綜合分析

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

根據表1,我們可以看出,隨著K的增大,運算時間也成約4倍遞增,然而當K=8時,BLAST運算時間約為改進VSM模型算法的50倍.

綜合以上運算效率和排名結果兩方面的分析,我們可以得出改進VSM算法維持了其計算效率的優越性,并改進了排名結果,提高了檢出率.

3 展望

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

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