基于主成分特征投影法建立利用16S rRNA基因序列進(jìn)行物種分類方法的研究

任 瑩， 畢春霞， 秦江楠， 郭雙雙， 王 斌， 宋旭霞， 閆志勇*

(1.青島大學(xué)醫(yī)學(xué)院 微生物教研室，山東 青島 266071；2.青島市市立醫(yī)院，山東 青島 266071)

?

基于主成分特征投影法建立利用16S rRNA基因序列進(jìn)行物種分類方法的研究

任 瑩1， 畢春霞2， 秦江楠1， 郭雙雙1， 王 斌1， 宋旭霞1， 閆志勇1*

(1.青島大學(xué)醫(yī)學(xué)院 微生物教研室，山東 青島 266071；2.青島市市立醫(yī)院，山東 青島 266071)

提出一種有別于系統(tǒng)發(fā)育樹的根據(jù)16S rRNA基因序列進(jìn)行物種分類的新方法。首先將基因的堿基字母形式轉(zhuǎn)換成數(shù)字形式，構(gòu)建多維向量。然后根據(jù)主成分分析方法將該向量向數(shù)據(jù)分布最大方向投影，將原數(shù)據(jù)用幾個(gè)“主成分”線性表出，而不丟失原數(shù)據(jù)的信息，采用主成分的顯示功能作出三維主成分特征投影視圖，達(dá)到分類的目的。在雙歧桿菌和腸球菌的分類識(shí)別中得到較好的應(yīng)用。

物種分類；堿基序列；向量；主成分分析

隨著測(cè)序技術(shù)的快速發(fā)展，現(xiàn)代細(xì)菌分類學(xué)多是建立在對(duì)16S rRNA基因序列的系統(tǒng)發(fā)育學(xué)分析基礎(chǔ)上，并結(jié)合多種微生物信息的多相分類學(xué)[1]。1990年Worse[2]利用16S rRNA序列的同源性研究原核生物的分類和系統(tǒng)發(fā)育，取得了巨大成功，2004年Li和Ren[3]首先采用16S rRNA序列同源性分析確定了1株發(fā)酵產(chǎn)氫細(xì)菌，顯示出16S rRNA序列同源性分析快速方便的優(yōu)勢(shì)。系統(tǒng)發(fā)育分析主要分為DNA序列比對(duì)、建立核苷酸替代模型、建立系統(tǒng)發(fā)育樹以及對(duì)樹的評(píng)估4個(gè)步驟[4]。但是，從現(xiàn)代生存物種的大分子中獲得的進(jìn)化歷史信息是不完全的，因此，所推斷出來的系統(tǒng)發(fā)育樹有一定程度的不確定性和假設(shè)性，從同一組數(shù)據(jù)常常推斷出不同的系統(tǒng)發(fā)育樹。目前系統(tǒng)發(fā)育樹的構(gòu)建思想主要包括基于距離的最短距離法、最大簡(jiǎn)約法、最大似然法以及貝葉斯樹估計(jì)方法[5]。其中最大似然法是選擇概率最高的建樹方法，也是目前最常用、最為準(zhǔn)確的方法，但其最大的缺點(diǎn)是計(jì)算復(fù)雜性高，不能有效排除基因突變等因素帶來的干擾信息，從而不能處理大規(guī)模的數(shù)據(jù)，與日益凸顯的“大數(shù)據(jù)”時(shí)代相矛盾。后基因組時(shí)代的到來及各種生物信息的增多，對(duì)物種分類的挑戰(zhàn)不斷增大。而主成分分析 (principal component analysis, PCA)的中心目的就是數(shù)據(jù)的降維，即用少數(shù)幾個(gè)主因子代替龐大的原始數(shù)據(jù)而不丟失有用信息，排除相互重疊的冗余信息[6]。目前主成分分析已經(jīng)廣泛應(yīng)用于人口統(tǒng)計(jì)學(xué)[7]、化學(xué)計(jì)量學(xué)[8]、數(shù)量地理學(xué)[9]等領(lǐng)域，在生物學(xué)領(lǐng)域[10-13]近年來也得到廣泛應(yīng)用。將該方法應(yīng)用到利用16S rRNA的分類尚未見文獻(xiàn)報(bào)道，本文擬利用16S rRNA基因序列轉(zhuǎn)換而成的向量，結(jié)合主成分特征投影方法，進(jìn)行物種的識(shí)別以及分類發(fā)育的初步探討。

1 材料與方法

1.1 主成分分析

PCA是一種能夠解決“維數(shù)災(zāi)難”的有效方法，既能降低數(shù)據(jù)復(fù)雜度，又能更好地分析理解數(shù)據(jù)。在計(jì)算機(jī)出現(xiàn)以前，PCA由于其計(jì)算繁瑣并沒有得到廣泛應(yīng)用，但是在計(jì)算機(jī)出現(xiàn)以后，其發(fā)展非常迅速，成為多元統(tǒng)計(jì)分析中至關(guān)重要的一種方法。它是通過正交變換，將原始數(shù)據(jù)向數(shù)據(jù)變化最大方向作投影，將數(shù)據(jù)分解成得分矩陣T和載荷矩陣PT，即X=TPT，這里的得分矩陣T即為原始變量的線性組合。實(shí)際應(yīng)用中，前幾個(gè)對(duì)數(shù)據(jù)分布影響最大的主成分用于分析，使高維數(shù)據(jù)降低到二維或者三維等低維空間，同時(shí)又不丟失有用信息。

主成分特征投影可通過PCA降維，實(shí)現(xiàn)由前幾個(gè)主成分重構(gòu)的PCA模型捕獲數(shù)據(jù)的大部分信息，還原原始基因序列的信息。相似度較高的對(duì)象在得分圖中傾向于聚集在一起，而相似度較低的對(duì)象將傾向于相互遠(yuǎn)離，實(shí)現(xiàn)聚類分析。一般地，前3個(gè)主成分所代表的信息占到原始數(shù)據(jù)變量的85%以上，這足以表示基因的分類信息。將4種堿基字母形式轉(zhuǎn)換成數(shù)字形式，構(gòu)成向量，一個(gè)列向量代表一條基因序列，采用主成分分析方法基于變量協(xié)方差矩陣對(duì)信息進(jìn)行處理、壓縮和抽提，再將數(shù)據(jù)進(jìn)行三維可視化主成分特征投影作圖。不同物種之間基因序列中堿基的差異明顯，同一物種間差異較小，從而導(dǎo)致在三維視圖中，不同的物種基因序列所處位置不一致，因此該方法可用于物種的分類鑒別。

1.2 奇異值分解

奇異值分解[14]是一種常用的主成分分析方式，可將實(shí)數(shù)矩陣分解成為3個(gè)矩陣的積X=USVT，其中主成分分析中的得分T=US，載荷為PT=VT，采用前3個(gè)得分矢量(PC1、PC2、PC3)作圖，即可得到可視化的三維主成分特征投影圖。

1.3 實(shí)例

以雙歧桿菌(Bifidobacterium)和腸球菌(Enterococcus)為例，研究該方法的可行性。實(shí)驗(yàn)所用基因序列來自GenBank中收錄的16S rRNA基因序列，從中篩選出序列長(zhǎng)度>1 400 bp的菌株(http://www.ncbi.nlm.nih.gov/pubmed)，計(jì)算平臺(tái)為MATLAB2014a，步驟如下：

1：將下載的基因序列轉(zhuǎn)變成數(shù)字形式，c改為1，a改為2，t改為3，g改為4，構(gòu)成列向量，并將數(shù)據(jù)對(duì)齊；

2：將數(shù)據(jù)進(jìn)行歸一化處理，其目的是歸納統(tǒng)一樣本的統(tǒng)計(jì)分布性；

3：將歸一化后的數(shù)據(jù)進(jìn)行奇異值分解，其模型為X=USVT；

4：求取主成分，PC=US；

5：作出可視化的三維主成分特征投影視圖。

同時(shí)，利用MEGA6.0軟件繪制系統(tǒng)發(fā)育樹。

2 結(jié)果與分析

分別采用1 300個(gè)基因和100個(gè)基因建立分類模型，圖1為取序列前1 300個(gè)基因時(shí)主因子變化示意圖，從圖1中可以看出第一主成分的貢獻(xiàn)最大，已經(jīng)代表大部分信息，分別采用第一、第二、第三主成分作圖可以將基因序列中幾乎所有信息抽提出來，用前3個(gè)主因子所做三維可視化圖(圖2)，可明顯看出2種細(xì)菌分屬2類，通過三維圖中距離的遠(yuǎn)近即可判別出親緣關(guān)系的遠(yuǎn)近，因此采用主成分分析方法可以將不同物種區(qū)別開來，并能夠判別親緣關(guān)系的遠(yuǎn)近。

同樣，采用主成分分析方法使用前100個(gè)基因進(jìn)行聚類分析(圖3)，發(fā)現(xiàn)也可以將物種明顯區(qū)別開來。

采用主成分特征投影法得到的分類視圖與MEGA6.0進(jìn)行發(fā)育樹結(jié)果(圖4)基本吻合，但是在實(shí)際應(yīng)用過程中發(fā)現(xiàn)，隨著數(shù)據(jù)量的增多，MEGA6.0發(fā)育樹計(jì)算變慢、時(shí)間長(zhǎng)，并且當(dāng)數(shù)據(jù)超過一定量時(shí)也會(huì)增加結(jié)果顯示的難度，而本方法采用數(shù)據(jù)降維的方式，對(duì)數(shù)據(jù)進(jìn)行抽提、壓縮，能夠有效解決上述問題。隨著大數(shù)據(jù)時(shí)代的到來，對(duì)數(shù)據(jù)進(jìn)行抽提、壓縮，使其更加直觀，勢(shì)必是今后發(fā)展的趨勢(shì)。

圖1 所取主因子數(shù)與特征值大小Fig.1 The figure of matrix's eigenvalue

圖2 1 300個(gè)基因PCA分類效果圖Fig.2 The figure of principal component projection when in 1 300 genes

圖3 100個(gè)基因PCA分類圖Fig.3 The figure of principal component projection when in 100 genes

圖4 利用N-J法構(gòu)建系統(tǒng)發(fā)育樹Fig.4 Neighbour-joining phylogenetic tree based on 16S rRNA gene sequences

3 討 論

本研究采用主成分分析方法使用前100個(gè)基因進(jìn)行聚類分析，發(fā)現(xiàn)也可以將物種明顯區(qū)別開來，但是與采用1 300個(gè)基因分類不同，腸球菌在三維視圖中分布相對(duì)集中，雙歧桿菌在PC1、PC2的二維視圖中也是相對(duì)集中。這可能是因?yàn)樵谛蛄袦y(cè)定過程中，100個(gè)基因測(cè)序誤差相對(duì)較小，也可能是前100個(gè)基因同類之間差異相對(duì)較小，因此同一類樣品的數(shù)據(jù)分布方向相對(duì)集中，不同類之間數(shù)據(jù)分布方向變化較大，三維視圖中同一類細(xì)菌分布集中，不同類之間差異明顯。

因此，該方法對(duì)基因測(cè)序準(zhǔn)確性以及對(duì)齊與否有一定要求。這在主成分分析的原理中亦能尋得答案，主成分分析是將原始數(shù)據(jù)向數(shù)據(jù)分布最大方向作圖，如果錯(cuò)誤信息過多，數(shù)據(jù)將不再投影到差異信息的方向，而是投影到錯(cuò)誤信息的方向，因此在接下來的研究中應(yīng)該研究異常值的剔除以及對(duì)齊基因的尋找優(yōu)化。

本文探究了采用主成分特征投影法對(duì)雙歧桿菌和腸球菌進(jìn)行分類，發(fā)現(xiàn)該方法可以明顯的區(qū)分這兩種物種，為采用多元統(tǒng)計(jì)方法進(jìn)行物種的歸類奠定了理論基礎(chǔ)，但是對(duì)于親緣關(guān)系更為接近的物種分類，有待于基因?qū)R算法的研究。

[1] 楊霞,陳陸,王川慶. 16S rRNA基因序列分析技術(shù)在細(xì)菌分類中應(yīng)用的研究進(jìn)展[J]. 西北農(nóng)林科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,02:55-60.

[2] Worse C R,Kandler O,Wheelis M L. Towards a natural system of organisma:Proposal for the domains archaea,bacteria, and eucaya[J]. Proc Natl Acad Sci,1990,87：5476-4579.

[3] Li Y F,Ren N Q,Yang C P,et al. Biohydrogen production behaviour and molecular characterization of a new species of anaerobic bacterium [A]. Anaerobic Digestion 2004- Proceedings of the 10th World Anaerobic Conference [C]. Montreal,Canada,2004.

[4] 唐曉嗣. 系統(tǒng)發(fā)育樹構(gòu)建中的EM算法[D].廣州：暨南大學(xué),2006.

[5] 馮思玲. 系統(tǒng)發(fā)育樹構(gòu)建方法研究[J]. 信息技術(shù),2009,(6):38-40，44.

[6] 王芳. 主成分分析與因子分析的異同比較及應(yīng)用[J]. 統(tǒng)計(jì)教育, 2003, 5(5): 14-17.

[7] 袁俊, 吳殿廷, 吳錚爭(zhēng). 中國(guó)農(nóng)村人口老齡化的空間差異及其影響因素分析[J]. 中國(guó)人口科學(xué), 2007, 3: 41-47.

[8] 許祿. 化學(xué)計(jì)量學(xué)方法[M].北京：科學(xué)出版社, 1995.

[9] 沈澤昊, 張新時(shí). 中國(guó)亞熱帶地區(qū)植物區(qū)系地理成分及其空間析局的數(shù)量分析[J]. 植物分類學(xué)報(bào), 2000, 38(4): 366-368.

[10]林麗, 李以康, 張法偉, 等. 青藏高原高寒矮嵩草草甸退化演替主成分分析[J]. 中國(guó)草地學(xué)報(bào), 2012, 34(1): 24-30.

[11]張麗英, 張正斌, 徐萍, 等. 黃淮小麥農(nóng)藝性狀進(jìn)化及對(duì)產(chǎn)量性狀調(diào)控機(jī)理的分析[J]. 中國(guó)農(nóng)業(yè)科學(xué), 2013, 47(5): 1013-1028.

[12]Frank D N, Amand A L S, Feldman R A, et al. Molecular-phylogenetic characterization of microbial community imbalances in human inflammatory bowel diseases[J]. Proceedings of the National Academy of Sciences, 2007, 104(34): 13780-13785.

[13]Makeig S, Jung T P, Bell A J, et al. Blind separation of auditory event-related brain responses into independent components[J]. Proceedings of the National Academy of Sciences, 1997, 94(20): 10979-10984.

[14]Wall M E, Rechtsteiner A, Rocha L M. Singular value decomposition and principal component analysis[M].A practical approach to microarray data analysis. Springer US, 2003: 91-109.

Methodology Establishment Using 16S rRNA Gene Sequences to Carry out Species Taxonomy Based on Principal Component Characteristics Projection Method

REN Ying1, BI Chun-xia2, QIN Jiang-nan1, GUO Shuang-shuang1,WANG Bin1, SONG Xu-xia1, YAN Zhi-yong1

(1.Teach. &Res.Div.ofMicrobiol.,Med.Coll.,QingdaoUni.; 2.QingdaoMuni.Hosp.,Qingdao266071)

This paper proposed a new species taxonomy method that differs from the phylogenetic tree based on 16S rRNA gene sequences. Firstly, the letter pattern of bases was converted into digital pattern, to build a multi-dimensional vector. Then according to the method of principal component analysis (PCA) the vector was projected into the maximum direction of database distribution, express the original data into "principal component" linearity, without loss of information of the original data, adopting the display function of principal component to make projection view of the principal component characteristics in three dimentions, to meet the goal of taxonomy. It had been fairly well applied in taxonomical recognition ofBifidobacteriumandEnterococcus.

species taxonomy; base sequence; vector; principal component analysis

山東省優(yōu)秀中青年科學(xué)家科研獎(jiǎng)勵(lì)基金項(xiàng)目(BS2011SW005)；山東省科技公關(guān)基金項(xiàng)目(2007GG3WZ05009)

任瑩 女，碩士研究生。研究方向?yàn)椴≡⑸飳W(xué)。E-mail:qdrenying@126.com

* 通訊作者。男，博士，副教授，碩士生導(dǎo)師。研究方向?yàn)椴≡⑸飳W(xué)。Tel:0532-83780059,E-mail:yanzhiyong@qdu.edu.cn

2015-04-27；

2015-05-24

Q78

A

1005-7021(2015)06-0105-04

10.3969/j.issn.1005-7021.2015.06.021