專利名稱：用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法。更具體地，本發(fā)明針對(duì)一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法，包括(i)選擇不同病毒的堿基序列作為參考序列，這些病毒的基因型或亞型是已知的，并且通過在所述參考序列的多重比對(duì)中計(jì)算序列之間的距離而獲得距離矩陣；以及(ii)開發(fā)一種判別方程，該判別方程可以對(duì)這些參考序列進(jìn)行分類，這是通過對(duì)通過該距離矩陣的多維定標(biāo)對(duì)所述參考序列成簇而獲得的聚簇執(zhí)行判別分析來實(shí)現(xiàn)的，接著根據(jù)所述判別方程對(duì)查詢序列的基因型與亞型進(jìn)行分類。
背景技術(shù)：
在理解趨異病毒的進(jìn)化方面，精確的基因分型(或分亞型)是關(guān)鍵。近來，公共數(shù)據(jù)庫里的病毒序列的數(shù)量的迅速增長被注意到。例如，NCBI基因庫(NCBI GenBank)擁有的 HIV-I與HCV序列條目幾乎每三年翻一番。這些病毒還顯示出非常好的基因型多樣性并且因此已經(jīng)被分類成組，被稱作基因型與亞型(Robertson等人,2000 ；Simmonds等人,2005)。因此，基于這些病毒株的序列相似性對(duì)它們進(jìn)行基因分型(或分亞型)，在理解它們的進(jìn)化、流行病學(xué)以及研發(fā)抗病毒療法或疫苗方面已經(jīng)成為最基本的步驟之一。傳統(tǒng)的分亞型方法包括以下(I)最近鄰法，尋找該查詢序列與被稱作參考的每一亞型的代表的最佳匹配；(2)系統(tǒng)發(fā)育方法，尋找該查詢序列分支至其上的單系群。由于這些亞型原本已經(jīng)被定義為單獨(dú)的聚簇群，所以這些直觀上合理的方法已經(jīng)得以廣泛使用并且對(duì)于許多案例而言十分成功。然而，隨著序列數(shù)目的漸增，觀察到不能被確切地分亞型的離群值或?qū)ζ涠赃@些方法不適宜的離群值。最近一份將這些不同的自動(dòng)分亞型方法與HIV-I序列作比較的報(bào)告顯示，除了亞型B與C之外，它們之中的相符性低于50% (Gifford R、de Oliveira T、Rambaut A、Myers RE、GaleCVΛ Dunn D、Shafer R、Vandamme AM、Kellam P、Pillay D UKCollaborative Group on HIV Drug Resistance:Assessmentof automatedgenotyping protocols as tools for surveillance of HIV-Igeneticdiversity. AIDS2006, 20:1521-1529)。該不相符性的原因之一要?dú)w結(jié)于由于重組而引起的增加的趨異性與復(fù)雜性。還應(yīng)注意到，在那些方法中，緊密關(guān)聯(lián)的亞型(B與D)或分享共同起源的亞型(A和CRF01_AE)顯示出較差的一致性。本發(fā)明人認(rèn)為，這一問題的根本是每一亞型的參考序列的數(shù)目太少。這些方法使用兩至四種手選的參考序列。它們是由各專家在高質(zhì)量的全基因組序列中仔細(xì)挑選的，是要盡量覆蓋每一亞型的多樣性。然而，利用每一亞型的本質(zhì)上小數(shù)目的參考序列，它們不能解決亞型預(yù)測(cè)的可信性；低E值的雙序列比對(duì)或高系統(tǒng)發(fā)育樹的高引導(dǎo)值(bootstrapvalue)表明單元操作的可靠性，但是就整體而言并不必然保證一個(gè)可信的亞型分類。對(duì)缺少統(tǒng)計(jì)置信測(cè)度這一問題的認(rèn)識(shí)帶來了 STAR的引入，這是一種基于特定位點(diǎn)打分矩陣的統(tǒng)計(jì)模型的方法，該特定位點(diǎn)打分矩陣是從每一亞型的多重序列比對(duì)(MSA)建立的。然而，其當(dāng)前的實(shí)施有一些限制它僅適用于HIV-I氨基酸序列，以小數(shù)目的參考(總共11個(gè)亞型的141種)為基礎(chǔ)，并且利用少于1000種序列進(jìn)行了測(cè)試。最近，已經(jīng)引入了新穎的基于核苷酸組成字符串的基因分型(或分亞型)方法。它的獨(dú)特在于它繞過了多重序列比對(duì)并且仍舊達(dá)到高精確度。然而，它也僅使用了 42種參考序列并且已經(jīng)用1156種序列進(jìn)行了測(cè)試?？紤]到這些病毒序列數(shù)目的爆炸式增長，這些傳統(tǒng)方法的測(cè)試案例非常少，最多萬分之一。因此，本發(fā)明的目的是要提供一種新穎的用于對(duì)公知的查詢序列的基因型或亞型進(jìn)行分類的方法。關(guān)鍵是在試圖對(duì)一種查詢序列進(jìn)行分類之前，評(píng)估每一亞型群的聚簇程度如何?？紤]這樣一個(gè)案例，其中這些參考序列大部分都被亞型很好地分開了，除了兩種或更多種亞型至少部分地重疊依賴少數(shù)參考的這些方法可能沒有注意到這一問題并且可能將高分分配給一種明顯的亞型。由于序列范圍內(nèi)的不同突變率，所以每一基因片段的系統(tǒng)發(fā)育動(dòng)力(phylogenetic power)也可能不同。這對(duì)于相對(duì)短的部分序列來說尤為關(guān)鍵。換言之，如果在基因分型(或分亞型)中僅考慮序列區(qū)域的一部分，那么即使這些本應(yīng)區(qū)別成 簇的、具有很好特征的參考也不能被分辨出。這些最近鄰法不能評(píng)估該背景分類模型的這種有效性，因?yàn)樗鼈儍H關(guān)注查詢與參考之間的比對(duì)，而不是參考與參考之間。REGA，基于樹的方法之一，關(guān)注該查詢是在由一組參考形成的聚族的內(nèi)部還是外部(deOliveira TDeforche K、Cassol S、Salminen Μ、Paraskevis D、SeebregtsC、Snoeck J>van Rensburg EJ>ffensing AM、van de Vijver DA、BoucherCA、Camacho R>Vandamme AM An automated genotyping system foranalysis ofHIV-Iand other microbial sequences. Bioinformatics 2005、21:3797-3800)。然而，就本發(fā)明人所知曉的，沒有工具定量地報(bào)道這樣一種測(cè)量。所以，本發(fā)明人提出一種方法，該方法基于這些參考序列之間的距離開發(fā)了這些背景分類模型，重新評(píng)估了它們對(duì)于每一查詢的有效性，并且就后驗(yàn)概率報(bào)告了基因型(或亞型)賦值的統(tǒng)計(jì)顯著性。如此，本發(fā)明的方法適合于其中許多參考序列可用的案例。本發(fā)明通過將主坐標(biāo)分析(PCoA)與線性判別分析(LDA)(兩者是使用生物科學(xué)中普遍的應(yīng)用能很好建立的統(tǒng)計(jì)工具)結(jié)合起來而實(shí)現(xiàn)這些目標(biāo)。PCoA (也稱為經(jīng)典多維定標(biāo)(MDS))，將這些序列標(biāo)繪在高維主坐標(biāo)空間，同時(shí)盡可能地盡力保持它們之間的距離關(guān)系。PCoA已經(jīng)廣泛地應(yīng)用于探索序列集中的全球趨勢(shì)，在系統(tǒng)發(fā)育分析方面對(duì)基于樹的方法進(jìn)行了補(bǔ)充。因?yàn)閬喰鸵呀?jīng)被定義為系統(tǒng)發(fā)育樹中的不同單系類群，所以如果選擇一種適當(dāng)?shù)母呔S，每一亞型應(yīng)該在MDS空間里形成良好分離的聚簇。在此類案例中，可以發(fā)現(xiàn)一組將這些聚簇分開的超平面并且與這些超平面相關(guān)的查詢可以得到分類。為了這一目的，本發(fā)明將LDA (—種直接的并且強(qiáng)大的分類方法)應(yīng)用于MDS坐標(biāo)并且將一種查詢分配給顯示出最聞的關(guān)系后驗(yàn)概率的基因型(或亞型)。這種概率在檢測(cè)任何需要仔細(xì)檢驗(yàn)的模糊案例時(shí)是有用的。本發(fā)明的方法通過留一法交叉驗(yàn)證(L00CV)來測(cè)試這些LDA模型，該驗(yàn)證可以用以通過檢測(cè)誤分類率來估測(cè)模型有效性。由于這些序列是由坐標(biāo)來表示的，因此還可以開發(fā)一種簡單的措施用以檢測(cè)基因型(或亞型)離群值。本發(fā)明人實(shí)質(zhì)上已經(jīng)利用所有來自NCBI基因庫(核苷酸)與GenPept (蛋白質(zhì))的HIV-I和HCV序列對(duì)本發(fā)明進(jìn)行了測(cè)試。披露內(nèi)容技術(shù)問題本發(fā)明的主要目的是提供一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法，包括(i)選擇不同病毒的堿基序列作為參考序列，這些病毒的基因型或亞型是已知的，并且通過在所述參考序列的多重對(duì)比中計(jì)算序列之間的距離而獲得距離矩陣；以及(ii)開發(fā)一種判別方程，該判別方程可以對(duì)這些參考序列進(jìn)行分類，這是通過對(duì)通過該距離矩陣的多維定標(biāo)對(duì)所述參考序列成簇而獲得的聚簇執(zhí)行判別分析而實(shí)現(xiàn)的，接著根據(jù)所述判別方程對(duì)一種查詢序列的基因型與亞型進(jìn)行分類。技術(shù)解決方案本發(fā)明的上述主要目的可以通過提供一種用于對(duì)查詢序列的基因型與亞型進(jìn)行 分類的方法來達(dá)到，包括(i)選擇不同病毒的堿基序列作為參考序列，這些病毒的基因型或亞型是已知的，并且通過在所述參考序列的多重對(duì)比中計(jì)算序列之間的距離而獲得距離矩陣；以及(ii)開發(fā)一種判別方程，該判別方程可以對(duì)這些參考序列進(jìn)行分類，這是通過對(duì)通過該距離矩陣的多維定標(biāo)對(duì)所述參考序列成簇而獲得的聚簇執(zhí)行判別分析而實(shí)現(xiàn)的，接著根據(jù)所述判別方程對(duì)一種查詢序列的基因型與亞型進(jìn)行分類。本發(fā)明的方法的步驟(i)可以進(jìn)一步包括從所述多重比對(duì)中除去插入缺失。另外，本發(fā)明的方法的步驟(ii)的多維定標(biāo)優(yōu)選地是一種主坐標(biāo)分析。此外，本發(fā)明的方法的步驟(ii)的判別分析可以選自不同的方法，比如線性判別分析、二次判別分析、最近鄰距離法、支持向量機(jī)或線性分類。有利效果本發(fā)明的方法可以被有效地用于通過分析快速進(jìn)化的病毒(比如HIV-I與HCV)的序列而對(duì)病毒的基因型或亞型進(jìn)行精確分類。另外，本發(fā)明的方法對(duì)核苷酸和蛋白質(zhì)(多肽)序列都適用。而且，可以應(yīng)用本發(fā)明的方法根據(jù)多態(tài)性標(biāo)記(比如SNP)的距離矩陣將個(gè)別受試對(duì)象分類成群組。附圖簡要說明圖I示出了根據(jù)本發(fā)明的用于對(duì)病毒的基因型(或亞型)分析進(jìn)行分類的方法的示意圖。這些球形表示已知被成簇為四種群簇A-D的序列，并且這些組群的分界面由隔離圓圈表示。每一群簇里的實(shí)心球形分別地表示參考序列，并且查詢序列由星形表示。由于查詢序列位于群簇B與D之間的分界面內(nèi)，所以難以查明該查詢序列的基因型(或亞型)。另一方面，可以通過最近鄰法來將查詢序列分配給最鄰近參考序列并且這種情況發(fā)生在群簇D中。根據(jù)最鄰近參考序列的距離，而不考慮已知分類方法的序列的聚簇模式，就該參考序列的選擇而言，這些結(jié)果可以并不穩(wěn)健(robust)。圖2示出了沿第一(VI)、第二(V2)以及第三(V3)主坐標(biāo)軸的HIV-I序列的示例性MDS示意圖。這些參考序列被示出為根據(jù)其亞型進(jìn)行了顏色編碼的小圓圈。為了清楚起見，沒有對(duì)亞型F-K進(jìn)行標(biāo)記。該查詢位于亞型B的中間(‘ + ’)。圖3示出了對(duì)每一基因片段而言通過MDS維數(shù)K示出的LOOCV錯(cuò)誤率。對(duì)于(a)HIV-I核苷酸、(b)HIV-I蛋白質(zhì)、(c)HCV核苷酸以及(d)HCV蛋白質(zhì)序列的每一基因片段而言，參考序列的預(yù)測(cè)基因型(或亞型)的LOOCV錯(cuò)誤率是通過使該MDS維數(shù)K從I到50進(jìn)行變化來進(jìn)行測(cè)量。一些顯示出與眾不同的較高錯(cuò)誤率的基因片段被標(biāo)記。與序列類型無關(guān)，這些錯(cuò)誤率在k=10后都達(dá)到穩(wěn)定期，這些錯(cuò)誤率在隨后的分析中被使用。圖4示出了沿基因片段的LOOCV錯(cuò)誤率的代表性滑動(dòng)窗口繪圖。這些LOOCV錯(cuò)誤率是沿(a)HIV-Ienv核苷酸與(b)HCV e2蛋白質(zhì)序列的基因片段在滑動(dòng)窗口中繪制的。對(duì)兩種情況而言，該MDS維數(shù)是設(shè)置在k=10。總類表示出于圖8與圖9中。圖5示出了用于HIV-1“主要”分析的離群值O的密度分布。在測(cè)試的161，440個(gè)案例中，根據(jù)本發(fā)明的方法的159，261個(gè)預(yù)測(cè)與LANL亞型信息(實(shí)線)相一致，而剩下的則不一致(虛線)。圖5是利用在R統(tǒng)計(jì)包中執(zhí)行的核密度估計(jì)函數(shù)來產(chǎn)生的。通過0>2濾出的部分標(biāo)為陰影。在過濾出很大部分的不一致案例的同時(shí)，一致性案例的丟失被最小化。圖6示出了 HIV-I超變異序列的離群度值的盒形圖。離群度(O)參數(shù)的盒形圖是針對(duì)由先前研究(Janini M、Rogers M、Birx DR、McCutchan FE Human immunodeficiency virus type IDNA sequencesgenetically damaged by hypermutationare often abundant in patient peripheralbolld mononuclear cells and maybe generated during near-simultaneousinfection and activation of CD4 (+)T cells.J Virol2001，75 (17) :7973-7986 ;Gandhi SK、Siliciano JD、Bailey JR、Siliciano RF、Blankson JN Role ofAP0BEC3G/F_mediated hypermutation in thecontrol of humanimmunodeficiency virus type Iin elite suppressors.J Virol2008，82 (6) : 3125-3130 ;Land AM、Ball TB、Luo M、PilonRm、Sandstrom P、Embree JE、Wachihi C、Kimani J、Plummer FA Human immunodeficiency virus(HIV)typelproviralhypermutation correlates with CD4count in HIV-infectedwomen from Kenya.JVirol2008,82(16) :8172-8182)報(bào)道的561種無功能性與1，519種功能性序列繪制的，這些研究明確地標(biāo)記出每一序列是否為“無功能性的”。圖7示出了本發(fā)明對(duì)HIV-I進(jìn)行分亞型的網(wǎng)頁服務(wù)器屏幕截圖。圖7(a)示出了輸入屏并且圖7(b)-(d)分別示出了輸出的第一頁到最后一頁。圖8在滑動(dòng)窗口中示出了針對(duì)HIV-I核苷酸(上圖)與蛋白質(zhì)(下圖)序列((a) env、(b)gag、(c) nef、(d)pol、(e) vif > (f) vpu)的 L00CV 錯(cuò)誤率。圖9在滑動(dòng)窗口中示出了針對(duì)HCV核苷酸(上圖)與蛋白質(zhì)(下圖)序列((a) utr、(b)arfp、(c)core、(d)el、(e)e2、(f)ns2> (g)ns3、(h)ns4a、(i)ns4b、(j)ns5a、(k)ns5b、(l)okamoto、(m)p7)的 L00CV 錯(cuò)誤率。

圖10示出了針對(duì)該HIV-I “主要”分析的離群度值的柱狀圖與L00CV錯(cuò)誤率。對(duì)于基于本發(fā)明的預(yù)測(cè)與LANL —致的離群度值的分布示出了以大約I. O為中心的尖峰(a)，而那些不一致的則示出了直到10.0的很長的尾巴(b)。在過濾掉低可信度的案例(離群度〈2. O)之后，對(duì)于不一致性預(yù)測(cè)(d)仍留下比一致性預(yù)測(cè)(c)相對(duì)更多的具有較高錯(cuò)誤率的案例。然而，它們的比例不大并且任何基于這些值的過濾方案都沒有被執(zhí)行過。最佳模式在下文中，將參考以下實(shí)例與附圖更詳細(xì)地描述本發(fā)明。這些實(shí)例與附圖僅給出用于說明本發(fā)明而不在于限制本發(fā)明?？傮w過程
本發(fā)明的方法通過創(chuàng)建該查詢與參考序列的多重序列比對(duì)(MSA)來開始該過程。不像常規(guī)的方法，本發(fā)明要求大量的參考，它們應(yīng)該具有高質(zhì)量并且具有謹(jǐn)慎指定的基因型(或亞型)。洛斯阿拉莫斯國家實(shí)驗(yàn)室(Los Alamos National Laboratory (LANL))數(shù)據(jù)庫分配 HIV-1 (http://www.hiv. Ianl. gov/) and HCV (http ://hcv. lanl. gov/)序列的這樣的MSA。LANL還提供有關(guān)該MSA中每一序列的亞型信息。在2007年發(fā)布的HIV-IMSAs中包括總共3，591種核苷酸與3，478種蛋白質(zhì)序列，而在HCV MSAs中總共有3，093種核苷酸與3，077種蛋白質(zhì)序列。應(yīng)該注意，對(duì)一些亞型而言，在該MSA中發(fā)現(xiàn)超過100種序列，同時(shí)有極少亞型僅包括少數(shù)參考序列。該樣品大小失衡是一個(gè)嚴(yán)重問題，但是本發(fā)明提出一種基于全局方差(global variance)的相當(dāng)具有啟發(fā)性的解決方案。為了與其他方法公平比較，本發(fā)明人決定將該查詢與已經(jīng)可從公共數(shù)據(jù)庫中得到的參考序列的MSA進(jìn)行比對(duì)，而不是自己創(chuàng)建MSA，由而對(duì)該參考MSA表示尊重。這樣做具有節(jié)省執(zhí)行時(shí)間的優(yōu)點(diǎn)，這對(duì)網(wǎng)絡(luò)服務(wù)器應(yīng)用程序很關(guān)鍵。對(duì)于這一步驟，使用hmmbuild、hmmcalibrate、andhmmalign(http://hmmer. janeIia. org/)這套程序。在使用一種PERL腳本去除該MSA中的插入缺失之后，使用具有 Jukes-Cantor 修正的 EMBOSS 程序包(http://emboss, sourceforge. net/)的distmat來計(jì)算這些序列間的配對(duì)距離矩陣。 下一步驟是所謂的主坐標(biāo)分析(PCoA)，它將該距離矩陣轉(zhuǎn)變?yōu)槠錁?gòu)成與所搜索的坐標(biāo)的內(nèi)積相等的矩陣。通過所得到的矩陣的奇異值分解，獲得直到指定的較低維的一組特征向量以及相關(guān)特征值。然后將配對(duì)歐氏距離近似于這些原始距離的那些序列的多維坐標(biāo)從包括這些特征向量與特征值的簡單矩陣運(yùn)算中恢復(fù)。每一特征值是沿由相應(yīng)特征向量定義的軸而獲得的方差量，也稱作主坐標(biāo)(PC)。為了方便，這些特征值按降序排列并且通過采用最高的少數(shù)幾個(gè)來達(dá)到維數(shù)降低。如果組內(nèi)變異是忽略不計(jì)的，則最高PCs的數(shù)目或該MDS維數(shù)k應(yīng)該最多是N-1，其中N是參考組的數(shù)目。然而，根據(jù)所考慮的序列區(qū)域，一種亞型可能顯示出復(fù)雜的聚簇模式，分為一個(gè)以上的群簇，比如亞-亞型。因此，本發(fā)明人采用一種經(jīng)驗(yàn)性方法，該方法針對(duì)從I至50范圍的k來調(diào)查這些參考序列的交叉驗(yàn)證誤差。這一步驟是利用R統(tǒng)計(jì)系統(tǒng)(http://www. r-project. org/)中的cmdscale來實(shí)現(xiàn)的(圖2示出了該MDS結(jié)果的一個(gè)示意圖)。然后，本發(fā)明的下一步驟是開發(fā)判別模型，這些判別模型根據(jù)他們的亞型對(duì)參考進(jìn)行最佳地分類并且根據(jù)這些模型給該查詢分配亞型成員(membership)。在此,可以想象應(yīng)用除其他以外的不同分類方法，比如K-最近鄰近法(K-NN)、支持向量機(jī)(SVM)、線性分類器。如果這種MDS步驟真正有效，則這些參考應(yīng)該根據(jù)其亞型成員而被很好地聚簇，并且因此諸如線性判別分析(LDA)或二次判別分析(QDA)這些最簡單的方法應(yīng)該有效。這兩者通過使高斯分布函數(shù)適于每一組的中心而起作用，兩者之間的不同之處在于是使用全局協(xié)方差(LDA)還是使用組協(xié)方差(QDA)。由于可以預(yù)計(jì)組內(nèi)偏差可能組與組之間不同，因此QDA可能更合適。然而，以上提到的樣品大小失衡問題阻礙了應(yīng)用QDA，因?yàn)閷?duì)于一些基因型(或亞型)而言，在小量參考情況下它變得不穩(wěn)定。另一方面，LDA通常應(yīng)用全局協(xié)方差至所有這些亞型并且因此針對(duì)這一問題可以更穩(wěn)健(robust)。盡管它不如QDA嚴(yán)謹(jǐn)，但是只要這些組偏差彼此之間不是過于不同，則這種啟發(fā)式方法運(yùn)行地相當(dāng)好。一旦基于這些參考序列計(jì)算這些線性判別，則屬于特定組的后驗(yàn)概率是作為從該查詢至該組中心的所謂的馬氏距離函數(shù)而給出的。對(duì)于該查詢，之后分配后驗(yàn)(MAP)估計(jì)的最大值，也就是，具有最大可能性的亞型。該后驗(yàn)概率是通過與每一亞型的參考數(shù)目成比例的前者來進(jìn)行衡量的。這一步驟是利用R統(tǒng)計(jì)系統(tǒng)(http://WWW.r-project. org/)中的MASS程序包的Ida來實(shí)現(xiàn)的。預(yù)測(cè)模型的交叉驗(yàn)證這些線性判別模型的有效性是通過這些參考序列的基因型(或亞型)成員的LOOCV來進(jìn)行評(píng)估。對(duì)于這些參考中的每一個(gè)而言，其基因型(或亞型)是通過從這些參考中的其余參考產(chǎn)生的模型來進(jìn)行預(yù)測(cè)的。誤分類錯(cuò)誤率(它是誤分類參考的數(shù)量與參與驗(yàn)證的參考總數(shù)量的比)是對(duì)該背景分類能力的一種敏感量。公共數(shù)據(jù)庫中的許多病毒序列并不是全基因組，而只覆蓋了一些基因或一個(gè)基因的一部分，并且因此它們的系統(tǒng)發(fā)育信號(hào)可以不同。因此，本發(fā)明人利用LOOCV重新評(píng)估了每一預(yù)測(cè)的分類能力。如果在針對(duì)一種給定查詢的MDS空間里這些參考序列得不到很好的辨析，則在LOOCV中會(huì)很明顯，導(dǎo)致高誤分類率。離群值檢測(cè) 即使通過亞型使這些參考以低LOOCV失誤率得以很好地分開，該查詢序列本身還是可能異常它可以是兩種或更多種亞型的復(fù)合，位于數(shù)種亞型的中間(一種重組體情況);它可以僅接近一種亞型群簇(針對(duì)這種亞型具有接近I的P值)但是遠(yuǎn)在該群簇邊界之外(一種趨異情況)。在多變量分析的領(lǐng)域內(nèi)，習(xí)慣是通過計(jì)算自樣品中心的馬氏距離并且通過將其與卡方分布進(jìn)行比較來檢測(cè)離群值。由于該馬氏距離已經(jīng)結(jié)合在LDA后驗(yàn)概率的計(jì)算中，因此本發(fā)明人提出一種有些不同的量，即，離群度0，它是從該查詢至與屬于沿該方向
的那種亞型的參考的最大趨異值有關(guān)的群簇中心的歐氏距離
j j 2'rrrrrrrr^ R\-\S
(Eq. I)其中XQ、Xe以及Xc分別是該查詢、這些參考之一以及該參考組S的中心的MDS向量。該組S包含所有屬于已經(jīng)將該查詢分類給其的基因型(或亞型)的所有參考序列。如果O小于1.0，則該查詢是很好地在該群簇內(nèi)部，否則就在外部。本發(fā)明人基于此開發(fā)了一種簡單的啟發(fā)式過濾器例如，可以將閾值設(shè)置在2. O以容許一些偏差。REGA還通過檢查樹形拓?fù)鋪韴?zhí)行離群值檢測(cè)方案以查看該查詢是在由參考序列組形成的群簇的內(nèi)部還是外部。重組體檢測(cè)的套合分析(Nested Analysis)用于表征重組體病毒株的標(biāo)準(zhǔn)過程包括沿該序列的靴掃描(bootscanning)以定位該重組點(diǎn)。它僅適用于長序列并且對(duì)于依賴于大樣品量的工具(比如本發(fā)明的方法)而言，要實(shí)用地通過互聯(lián)網(wǎng)而服務(wù)，它花費(fèi)時(shí)間太多，除非采用具有數(shù)百CPU的群簇場(cluster farm)。與其執(zhí)行靴掃描，本發(fā)明人通過以下途徑解決了該重組問題(a)對(duì)于包括多于一個(gè)基因的查詢而言，逐基因預(yù)測(cè)亞型；(b)以一種包括重組參考序列的“套合”方式對(duì)該分析進(jìn)行再迭代。HIV-I與HCV包含順序的10個(gè)基因并且因此對(duì)整個(gè)基因組序列進(jìn)行逐基因分析不會(huì)花費(fèi)比單個(gè)基因分析長10倍的時(shí)間。如果不同亞型被以高可信度分配給了一種查詢的不同基因構(gòu)成部分，則暗示了一種重組情況。對(duì)于一些重組體，斷點(diǎn)可以發(fā)生在一個(gè)基因的中間。在此類情況中，有可能的是，分類的后驗(yàn)概率不是僅受一種亞型支配，但是第二個(gè)左右會(huì)具有一個(gè)不可忽略的P值。本發(fā)明人通過對(duì)具有大于O. Ol的P值的亞型以及相關(guān)的重組體亞型予以注意，以一種“套合”方式對(duì)該預(yù)測(cè)過程進(jìn)行再迭代。例如，如果A組或G組的P值大于O. 01，則這些參考包括CRF02_AG組。網(wǎng)頁服務(wù)器開發(fā)已經(jīng)研發(fā)出接受核苷酸序列作為一種查詢并且預(yù)測(cè)該查詢的每一基因片段的基因型(或亞型)的阿帕奇(Apache)網(wǎng)頁服務(wù)器，每個(gè)HIV-I與HCV有一個(gè)網(wǎng)絡(luò)服務(wù)器。接受氨基酸序列作為一種查詢的相應(yīng)蛋白質(zhì)版本也已經(jīng)得以開發(fā)。這些可以在http://WWW.muldas. org/MuLDAS/上免費(fèi)取得。以PERL編寫的每一 CGI程序封裝了已從HMMER、EMB0SS以及R的各自發(fā)布網(wǎng)站上下載的組件程序。由于距離矩陣的運(yùn)算耗費(fèi)許多運(yùn)行時(shí)間，因此本發(fā)明人將該任務(wù)分割為數(shù)個(gè)(典型地是四個(gè))計(jì)算節(jié)點(diǎn)，其中每一個(gè)計(jì)算節(jié)點(diǎn)并行地計(jì)算這些行的多個(gè)部分，并且這些結(jié)果通過主節(jié)點(diǎn)進(jìn)行整合。在英特爾至強(qiáng)CPU Linux盒(IntelXeon CPU Linux box)上，對(duì)一段1000-bp的HIV-I核苷酸序列的典型亞型預(yù)測(cè)要花費(fèi)大約20秒。這些網(wǎng)頁服務(wù)器報(bào)告該查詢的MAP基因型(或亞型)以及每一亞型的后驗(yàn)概率 (posterior P)、這些預(yù)測(cè)模型的留一法交叉驗(yàn)證結(jié)果、以及離群值檢測(cè)結(jié)果(圖7的屏幕截圖)。該查詢的3D示意圖與前三個(gè)PC中的參考是以PNG格式給出并且描述該查詢的所有PC以及這些參考的XML文件可以下載,用于隨后利用GGobi (http://www. ggobi. org/)的動(dòng)態(tài)互動(dòng)可視化(Fig. 2)。這對(duì)于可視地檢查聚簇的質(zhì)量以及對(duì)于確定可以導(dǎo)致識(shí)別出潛在的新型或重組體的離群檢測(cè)結(jié)果來說尤其有用。對(duì)于HIV-1，以上描述的“套合”分析被進(jìn)行再迭代并且該結(jié)果也被報(bào)告。該網(wǎng)站還運(yùn)行存儲(chǔ)了 HIV-I亞型與HCV基因型的預(yù)計(jì)算結(jié)果的數(shù)據(jù)服務(wù)器，這些結(jié)果是利用與這些預(yù)測(cè)服務(wù)器完全一樣的方法預(yù)測(cè)的。定期地(典型地是每天)下載NCBI基因庫與GenP^t中HIV-I或HCV的所有新條目，并且預(yù)測(cè)它們的基因型(或亞型)并存儲(chǔ)在這些數(shù)據(jù)庫中?？梢酝ㄟ^NCBI GI編號(hào)或主入藏號(hào)(primary accessions)檢索這些結(jié)果。還以利用由諸如后驗(yàn)概率、L00CV率、離群度、基因型(或亞型)、或基因片段這些系統(tǒng)計(jì)算的性質(zhì)來查詢這些條目。該檢索的結(jié)果包括從LANL數(shù)據(jù)庫里讀取的基因型(或亞型)信息，如果有的話。結(jié)果本發(fā)明的方法是利用從NCBI基因庫與GenP印t下載的HIV-I與HCV的序列數(shù)據(jù)集進(jìn)行測(cè)試的。針對(duì)還沒有用作參考序列的158，834種HIV-I序列(包括8，832種重組體)以及48，720種HCV序列，從LANL網(wǎng)站上檢索核苷酸序列的亞型信息并且將這些亞型信息用于探尋出源自該核苷酸序列的蛋白質(zhì)序列的亞型信息。對(duì)于一些序列而言，這些基因型/亞型是由最初提交者給出的或由LANL分配。這些測(cè)試數(shù)據(jù)集的基因型(或亞型)命名法HIV-I序列被分組為M (主要(main))組、N (非主要(non-main))組、U (未經(jīng)分類(unclassified))組、O (外類群(outgroup))組。多數(shù)可用的序列屬于M組。由于N組與O組距離M組非常遠(yuǎn)，因此M組的亞型在包括這些遠(yuǎn)離組的MDS示意圖中不能得到很好的解析。因此，本發(fā)明人集中于將M組序列分類為亞型A-D、F-H、J以及K。在M組的亞型中，有時(shí)將A與F進(jìn)一步分別地分開為亞-亞型Al與A2以及Fl與F2。
然而，在LANL數(shù)據(jù)庫中仍有一些新序列在亞型等級(jí)上被報(bào)道。甚至對(duì)于包括在由LANL產(chǎn)生的MSA中的序列也是這種情況。利用本發(fā)明針對(duì)相對(duì)短的序列解析亞-亞型要求一種僅使用相關(guān)亞型序列的“套合”分析。由于這些原因，本發(fā)明人沒有試圖去區(qū)別亞-亞型并且在亞型等級(jí)上對(duì)它們進(jìn)行分類。M組序列的不同亞型可以重組來形成一種新株。如果在三個(gè)以上流行病學(xué)上獨(dú)立的病人中發(fā)現(xiàn)這些株，則稱它們?yōu)榱餍兄亟M形式((circulating recombinant forms) CRFs)。在這些 CRF 中，CRF01_AE 由 A 與現(xiàn)在已滅絕的E株重組形成，并且構(gòu)成一個(gè)與A亞型不同的大家族。M組與CRF01_AE亞型已被稱為“主要”亞型并且本發(fā)明的方法針對(duì)它們作為“主要”分析來進(jìn)行。表I列出了按照亞型以及所有已經(jīng)被LANL分類為“主要”組的測(cè)試序列的基因片段統(tǒng)計(jì)的分項(xiàng)數(shù)據(jù)(相應(yīng)的蛋白質(zhì)序列參考表2)。該分布遠(yuǎn)不一致，代表了研究偏差屬于亞型H、J以及K的序列稀少；特別對(duì)于諸如vif與vpr的輔助蛋白而言，非B株過于稀少，以至于不能精確評(píng)估該分類。 表I. HIV-IM組以及CRF01_AE核苷酸序列的基準(zhǔn)測(cè)試的總結(jié)性統(tǒng)計(jì)(a)過濾之前每一亞型的基因片段的數(shù)目
權(quán)利要求
1.一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法，包括 (i)選擇不同病毒的堿基序列作為參考序列，這些病毒的基因型或亞型是已知的，并且通過在所述參考序列的多重比對(duì)中計(jì)算序列之間的距離而獲得一種距離矩陣；以及 (ii)開發(fā)一種判別方程，該判別方程可以對(duì)這些參考序列進(jìn)行分類，這是通過對(duì)通過該距離矩陣的多維定標(biāo)對(duì)所述參考序列成簇而獲得的聚簇執(zhí)行判別分析而實(shí)現(xiàn)的，接著根據(jù)所述判別方程對(duì)一種查詢序列的基因型與亞型進(jìn)行分類。
2.根據(jù)權(quán)利要求I所述的方法，其中所述步驟(i)進(jìn)一步包括從所述多重比對(duì)中除去插入缺失。
3.根據(jù)權(quán)利要求I所述的方法，其中所述步驟(ii)的所述多維定標(biāo)是一種主坐標(biāo)分析。
4.根據(jù)權(quán)利要求I所述的方法，其中所述步驟(ii)的所述判別分析是選自包括線性判別分析、二次判別分析、最近鄰點(diǎn)距離法、支撐向量機(jī)以及線性分類器的組。
全文摘要
本發(fā)明涉及一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法。更具體地，本發(fā)明針對(duì)一種用于對(duì)查詢序列的基因型與亞型進(jìn)行分類的方法，包括(i)選擇不同病毒的堿基序列作為參考序列，這些病毒的基因型或亞型是已知的，并且通過在所述參考序列的多重比對(duì)中計(jì)算序列之間的距離而獲得一種距離矩陣；以及(ii)開發(fā)一種判別方程，該判別方程可以對(duì)這些參考序列進(jìn)行分類，這是通過對(duì)通過該距離矩陣的多維定標(biāo)對(duì)所述參考序列成簇而獲得的聚簇執(zhí)行判別分析而實(shí)現(xiàn)的，接著根據(jù)所述判別方程對(duì)一種查詢序列的基因型與亞型進(jìn)行分類。
文檔編號(hào)G06F19/24GK102884203SQ201080066436
公開日2013年1月16日 申請(qǐng)日期2010年8月13日 優(yōu)先權(quán)日2010年2月26日
發(fā)明者金尚洙 申請(qǐng)人:崇實(shí)大學(xué)校產(chǎn)學(xué)協(xié)力團(tuán)