本發(fā)明涉及信息技術(shù)領(lǐng)域和生物領(lǐng)域，具體來(lái)說(shuō)涉及一種基于不同尺度tuple詞頻的微生物高通量測(cè)序數(shù)據(jù)分析協(xié)議

背景技術(shù)：

微生物群落是地球上生物多樣性最為豐富的資源，廣泛存在于各種自然環(huán)境中，如土壤、人體皮膚和消化系統(tǒng)中。環(huán)境中的微生物傳統(tǒng)意義上包括細(xì)菌、真菌、病毒以及一些古生菌，不同的環(huán)境中微生物的物種豐度組成多樣性以及微生物功能多樣性都存在著巨大的差異。為了更好地洞悉微生物在不同的微生物環(huán)境中所起的功能作用以及更好地理解微生物和環(huán)境之間的關(guān)系，對(duì)環(huán)境中所有的微生物基因組研究是極有必要的。

傳統(tǒng)的測(cè)序方法獲取的微生物數(shù)量很少，不能從整體上去描述微生物群落之間的結(jié)構(gòu)差異。而高通量測(cè)序技術(shù)可以獲得更完整、更準(zhǔn)確的微生物群落結(jié)構(gòu)，因此高通量測(cè)序技術(shù)逐漸成為研究學(xué)者們對(duì)微生物群落的比較研究的一個(gè)強(qiáng)有力的工具，通過高通量測(cè)序技術(shù)我們可以直接從環(huán)境中獲取大量的微生物測(cè)序樣本，基于這些樣本，微生物群落的比較方法大量被提取，其中主要包括基于16S rRNA的方法、基于配準(zhǔn)的序列比較方法，如Smith-Waterman算法和Blast算法、基于k-tuple的頻度統(tǒng)計(jì)方法。然而基于16S rRNA的方法，在微生物群落的分析比較上，存在很大的局限性，能得到的微生物群落構(gòu)成信息和物種分布范圍都很有限?；诟咄繙y(cè)序技術(shù)獲取的微生物測(cè)序數(shù)據(jù)，許多微生物的基因是未知的，現(xiàn)如今的微生物參考數(shù)據(jù)庫(kù)是極其不完備的，而基于配準(zhǔn)的方法高度依賴已知數(shù)據(jù)庫(kù)或已知基因，這就使得配準(zhǔn)的準(zhǔn)確性和完整性大大降低。

相比于基于配準(zhǔn)的方法，基于無(wú)需比對(duì)的方法克服了高度依賴參考數(shù)據(jù)庫(kù)的不足，為基因間的比較提供了更好的選擇。k-tuple方法是最具代表性的無(wú)需比對(duì)方法，基于k-tuple的頻度統(tǒng)計(jì)方法在微生物群落分析比較方面主要集中在長(zhǎng)度較短的tuple層面上(2-10bp)，結(jié)合概率背景統(tǒng)計(jì)模型和相異度度量方法，通過非監(jiān)督的聚類方法，在衡量微生物群落的差異性方面具有優(yōu)良的性能。然而目前基于這種短k-tuple的方法，只能建立tuple分布的總體統(tǒng)計(jì)模型，找出群落與群落之間的關(guān)系，度量其整體相異度。但是具體來(lái)說(shuō)是哪些特征序列、哪些微生物/基因序列造成群落間的這種差異和樣本類別的分組是k-tuple統(tǒng)計(jì)模型無(wú)法解決的問題。所以通過無(wú)監(jiān)督的聚類方法在對(duì)微生物群落進(jìn)行比較上就顯得不那么完整，而針對(duì)無(wú)監(jiān)督聚類獲取的樣本類別，通過有監(jiān)督的模式分類能夠進(jìn)一步的識(shí)別出不同類別高通量測(cè)序數(shù)據(jù)的特異性tuple，能夠?yàn)榭坍嫴煌悇e的微生物群落具體差異以及尋找生物標(biāo)記物提供重要的參考信息。

眾所周知，生物樣本是由A、C、G、T四種堿基組成的基因序列，k-tuple是指長(zhǎng)度為k的連續(xù)字符串序列。所以一個(gè)測(cè)序樣本的k-tuple頻度特征向量的維度就是4^k維。之前有研究表明，來(lái)自同一個(gè)基因組的k-tuple頻度相近，能夠建立tuple分布的總體統(tǒng)計(jì)模型，但不同基因組的k-tuple頻度有很大區(qū)別。在基于k-tuple頻度的無(wú)需比對(duì)的研究方法中，集中于短tuple(2-10bp)上，且應(yīng)用于無(wú)監(jiān)督的樣本聚類上表現(xiàn)比較優(yōu)越。因此，基于短k-tuple頻度的相異度距離度量方法D₂被提出用來(lái)評(píng)估比較兩個(gè)微生物群落樣本之間的相異度。此后，在D₂基礎(chǔ)上又衍生出和為了更好地應(yīng)用到高通量測(cè)序數(shù)據(jù)上，通過歸一化處理改進(jìn)的和相繼被提出用于比較樣本之間的相異度。

用和計(jì)算距離時(shí)需要結(jié)合合適的背景模型來(lái)進(jìn)行建模。在之前的研究中，用到的是定階次馬爾科夫模型和基于變階次的馬爾科夫模型。然而由于微生物群落是由各種各樣不同種類的微生物基因組混合組成，很難用幾個(gè)確定的階次模擬背景模型，而且需要手動(dòng)設(shè)定模型的階次，然后去集中對(duì)不同階次模型對(duì)聚類結(jié)果的優(yōu)良效果做評(píng)估，工作量和計(jì)算成本都非常高。對(duì)于定階次馬爾科夫模型，階次越高模型越準(zhǔn)確，然而階次越高，需要的數(shù)據(jù)量也越多，一般情況下，我們獲取的數(shù)據(jù)量是很難滿足需求的。而且基于變階次的馬爾科夫模型在對(duì)選擇模型階次時(shí)，對(duì)其構(gòu)建的前綴樹進(jìn)行減枝的過程中，需要手動(dòng)設(shè)定閾值，大大提高了模型的不精確性和計(jì)算的復(fù)雜性。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的主要目的在于克服現(xiàn)有技術(shù)中的上述缺陷，提出一種通過變尺度tuple頻度來(lái)區(qū)分微生物的群落類別，且基于獲取的群落類別，能夠找出區(qū)分群落類別的特異性信息的基于不同尺度tuple詞頻的微生物高通量測(cè)序數(shù)據(jù)分析協(xié)議。

本發(fā)明采用如下技術(shù)方案：

基于不同尺度tuple詞頻的微生物高通量測(cè)序數(shù)據(jù)分析協(xié)議，其特征在于，包括如下步驟：

步驟1：獲取宏基因組樣本的2-10bp的短tuple高通量測(cè)序數(shù)據(jù)，采用插值上下文馬爾科夫模型進(jìn)行建模微生物群落的背景基因組，再采用無(wú)監(jiān)督的聚類方法來(lái)比較宏基因組樣本，得出宏基因組樣本的類別信息；

步驟2：基于步驟1)中聚類得出的類別信息，將≥30bp的長(zhǎng)tuple作為特征，采用有監(jiān)督的樣本分類方法找出描述宏基因組樣本類別的特異性特征長(zhǎng)tuple序列。

優(yōu)選的，所述步驟1)具體包括如下：

步驟1.1:獲取宏基因組樣本的高通量測(cè)序數(shù)據(jù)，生成該宏基因組樣本的短tuple特征頻度向量，對(duì)每個(gè)宏基因組樣本中出現(xiàn)的長(zhǎng)度為2-10bp的tuple的頻度進(jìn)行統(tǒng)計(jì)，并生成相應(yīng)宏基因組樣本的頻度向量；

步驟1.2:采用插值上下文馬爾科夫模型對(duì)微生物群落的背景基因組進(jìn)行建模，估計(jì)頻度向量中每一個(gè)tuple的馬爾科夫概率；

步驟1.3：計(jì)算各個(gè)宏基因組樣本頻度向量間的相異度距離，生成一個(gè)宏基因組樣本間的相異度矩陣；

步驟1.4：根據(jù)相異度矩陣生成一個(gè)聚類樹，用于判斷宏基因組樣本與樣本之間的關(guān)系，找出樣本的類別信息。

優(yōu)選的，所述步驟1.1中，tuple特征定義為宏基因組樣本中可能出現(xiàn)的字符串組合，并選擇長(zhǎng)度為2-10bp的字符串組合作為所述tuple特征。

優(yōu)選的，步驟1.2中，計(jì)算tuple的馬爾科夫概率具體方法如下：

步驟1.2.1：基于宏基因組樣本的頻度向量的互信息量構(gòu)建上下文序列樹；

步驟1.2.2：基于上下文序列樹計(jì)算每個(gè)tuple的馬爾科夫概率。

優(yōu)選的，步驟1.2.1具體如下：基于宏基因組樣本的頻度向量，將k長(zhǎng)度tuple中的每一列的字符放在一個(gè)向量中，形成A₁,A₂,…,A_k k個(gè)向量，分別計(jì)算前面k-1個(gè)向量與最后一個(gè)向量之間的互信息量，互信息量的公式如下：

$I(A_w,B)=\underset{i}{\Σ}\underset{j}{\Σ}P(a_i,b_j)*log\frac{P\left(a_i\right)*P\left(b_j\right)}{P(a_i,b_j)}$

其中,w＝1,2,…,k-1；B＝A_k；a_i,b_j表示向量A,B中的變量；P(a_i,b_j)表示a_i,b_j在對(duì)應(yīng)向量中同時(shí)出現(xiàn)的聯(lián)合概率；P(a_i)表示a_i在對(duì)應(yīng)向量中出現(xiàn)的概率；

找出與B互信息量最大的向量A_w，將此向量對(duì)應(yīng)的下標(biāo)位置作為上下文序列樹的頂點(diǎn)；然后按照該向量中出現(xiàn)的四種不同字符(A,C,G,T)將所有tuple分成四組，；最后對(duì)四組tuple向量矩陣，按照公式中的互信息量公式分別計(jì)算互信息量，找出每組中與B互信息量最大的向量A_s，s＝1,2,…,w-1,w+1,…,k-1，將此向量對(duì)應(yīng)的下標(biāo)位置作為上下文序列樹的對(duì)應(yīng)葉子的子節(jié)點(diǎn)(A,C,G,T)；依次繼續(xù)下去，直到找到最后一個(gè)跟當(dāng)前狀態(tài)向量關(guān)聯(lián)性最大的向量，整個(gè)上下文序列樹構(gòu)建完畢。

優(yōu)選的，在步驟1.2.2中，每個(gè)tuple的馬爾科夫概率公式如下：

P(c₁c₂…c_k)＝P_ICM(c₁)P_ICM(c₂|c₁)…P_ICM(c_k|c₁c₂…c_k-1)

其中，c₁c₂…c_k表示k-tuple序列，P_ICM(c_k|c₁c₂…c_k-1)表示上下文序列c₁c₂…c_k-1轉(zhuǎn)移到當(dāng)前狀態(tài)c_k的ICM轉(zhuǎn)移概率。

優(yōu)選的，對(duì)于計(jì)算上述每個(gè)ICM轉(zhuǎn)移概率，針對(duì)所述k-tuple序列c₁c₂…c_k，應(yīng)用ICM馬爾科夫模型構(gòu)建的上下文序列樹中找出與當(dāng)前狀態(tài)c_k關(guān)聯(lián)性程度從大到小排序的重要位置，重新構(gòu)建其上下文序列，具體如下：要構(gòu)建r階馬爾科夫模型，r≤k-1，從上下文序列樹中找出與當(dāng)前狀態(tài)c_k關(guān)聯(lián)性程度從大到小排序的重要位置對(duì)應(yīng)的狀態(tài)分別為c₃,c₄…,c_r，組成插值上下文序列M_r,那么構(gòu)建ICM的概率模型如下：

P_ICM(c_k|c₁c₂c₃…c_k-1)＝P_ICM,r(c_k|M_r)；

P_ICM,r(c_k|M_r)＝λ_r*P(c_k|M_r)+(1-λ_r)*P_ICM,r-1(c_k|M_r-1)；

$P\left(c_k\right|M_r)=\underset{a\∈A,C,G,T}{\Σ}\frac{N(M_r,A_k=c_k)}{N(M_r,A_k=a)};$

其中，*表示乘積，λ_m表示m階馬爾科夫模型概率所占的權(quán)重系數(shù)；N(M_r,X_k＝c_k)表示所有的k-tuple中插值上下文序列是M_r，第k個(gè)位置是c_k的所有tuple的頻度之和，對(duì)于上述馬爾科夫模型概率所占的權(quán)重系數(shù)的計(jì)算公式為：

其中所述C表示樣本閾值，它由赤池信息量準(zhǔn)則AIC^R(C)確定，具體公式如下：

AIC^R(C)＝-2λ(S；M_k)+2|M_IMM,k,C|；

其中，λ(S；M_k)表示樣本S的偽似然度，計(jì)算公式為：

$\mathrm{\λ}(S;M_k)=\underset{c_1,\mathrm{...}{c}_k\∈(A,C,G,T)}{\Σ}N(c_1,\mathrm{...}{c}_{k-1}){\mathrm{logP}}_{ICM}\left(c_k\right|c_1,\mathrm{...}{c}_{k-1});$

|M_IMM,k,C|表示模型自由參數(shù)的個(gè)數(shù)，當(dāng)AIC^R(C)值最小時(shí)算出的C值作為樣本的閾值；

所述q表示兩個(gè)字符串之間差異度的卡方檢驗(yàn)值，其計(jì)算原理如下：

${\mathrm{\Δ}}_r\left(M_r\right)=\underset{a\∈A,C,G,T}{\Σ}\frac{{\left(N\right(M_r,a)-E(M_r,a\left)\right)}^2}{E(M_r,a)};$

E(M_r,a)＝N(M_r)P_ICM，r(a|M_r)；

其中，N(M_r,a)，E(M_r,a)分別表示字符串頻度的實(shí)際值和理論值，將q＝Δ_r(M_r)作為卡方值，自由度為3，作為卡方檢驗(yàn)的指標(biāo)參數(shù)。

優(yōu)選的，步驟1.3中，應(yīng)用不同的相異度距離度量方法計(jì)算各個(gè)宏基因組樣本頻度向量間的相異度距離，所用到的相異度距離度量方法包括和計(jì)算公式如下：

$D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}};$

$d_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{Y,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}});$

$D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{n_X{p}_{X,i}}\sqrt{n_Y{p}_{Y,i}}}$

$d_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}^2}{n_X{p}_{X,i}}}\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{Y,i}^2}{n_Y{p}_{Y,i}}}});$

其中，和都是計(jì)算兩個(gè)樣本之間的相異度的一種距離度量方法；表示樣本X的頻度向量，表示樣本Y的頻度向量；稱為中心化過程，i＝1,2,…,4^k；C_X,i和C_Y,i分別表示X和Y樣本中第i個(gè)tuple出現(xiàn)的頻度；n_X表示樣本X中tuple個(gè)數(shù)的總和，n_Y表示樣本Y中tuple個(gè)數(shù)的總和；p_X,i和p_Y,i分別表示在插值上下文馬爾科夫背景模型下，樣本X中第i個(gè)tuple的馬爾克夫概率和樣本Y中第i個(gè)tuole的馬爾科夫概率；

如果一個(gè)數(shù)據(jù)集中有n個(gè)樣本，根據(jù)相異度距離度量公式計(jì)算出的兩兩樣本間的相異度，生成一個(gè)n*n維的相異度距離矩陣，這個(gè)矩陣定義如下：

N(n,n)＝(d(x,y))_n×n,d(x,y)＝d(y,x),d(x,x)＝0

其中，d(x,y)是兩個(gè)宏基因組樣本的相異度距離，如果不同樣本間的距離越小，那么d(x,y)的值就越??；d(x,x)表示相同樣本間的距離為0。

優(yōu)選的，在步驟1.4中，在n*n相異度矩陣的基礎(chǔ)上，根據(jù)非加權(quán)平均法層次聚類算法計(jì)算出兩個(gè)群組的相異度距離，其定義如下：

$d(C_i,C_j)=\frac{1}{\left|C_i\right|\·\left|C_j\right|}\underset{x\∈C_i}{\Σ}\underset{y\∈C_j}{\Σ}d(x,y)$

d(x,y)是兩個(gè)宏基因組樣本的相異度距離，|C_i|和|C_j|表示兩個(gè)群組的大小，也就是群組中樣本的個(gè)數(shù)，i,j＝1,2,…,n；由兩兩群組的相異度距離得到聚類樹，從聚類樹中可以直觀的看出群落中各個(gè)樣本之間的結(jié)構(gòu)關(guān)系，得出樣本之間的類別信息。

優(yōu)選的，所述步驟2具體包括以下子步驟：

步驟2.1：對(duì)樣本中出現(xiàn)的長(zhǎng)度為40bp的長(zhǎng)tuple的頻度進(jìn)行統(tǒng)計(jì)，并生成相應(yīng)樣本的頻度向量；

步驟2.2：對(duì)每個(gè)樣本的tuple頻度向量進(jìn)行并行處理，生成所有樣本的長(zhǎng)tuple頻度向量矩陣，然后過濾掉冗余的特征；

步驟2.3：基于所述步驟1中獲取的樣本類別信息，應(yīng)用過濾好的樣本特征對(duì)樣本進(jìn)行有監(jiān)督分類，找到對(duì)分類效果具有很強(qiáng)辨識(shí)性的特異性tuple特征；

步驟2.4：基于步驟2.3獲得的特異性特征，用留一法(LOOCV)來(lái)驗(yàn)證和評(píng)估分類器的準(zhǔn)確率。

優(yōu)選的，在步驟2.2中，將需要分類樣本的tuple頻度向量合并在一起，生成一個(gè)tuple頻度向量矩陣A，A表示為M×N的頻度矩陣，其中N表示樣本數(shù)量，M表示特征維度。

優(yōu)選的，在步驟2.3中，基于步驟1中獲取的樣本類別信息，選定訓(xùn)練集和測(cè)試集樣本，在訓(xùn)練集中選定當(dāng)前類別和目標(biāo)類別，然后根據(jù)對(duì)稱不確定性大于設(shè)定的閾值時(shí)，把冗余的tuple序列特征過濾移除，得到一些類別特異性候選特征，對(duì)稱不確定性定義如下：

其中，N_X表示當(dāng)前類別組成的X樣本集中tuple特征出現(xiàn)的頻度；sum(N_X)表示由當(dāng)前類別組成的X樣本集中特征出現(xiàn)的頻度之和；sum(N_Y)表示目標(biāo)類別組成的Y樣本集中特征出現(xiàn)的頻度之和；n(X)和n(Y)分別表示X和Y樣本集中樣本的個(gè)數(shù)；θ表示X和Y之間對(duì)稱不確定性的閾值；

采用SVM分類器，對(duì)樣本進(jìn)行有監(jiān)督分類，找到能夠描述微生物群落內(nèi)部具有差異性的特異性特征。

優(yōu)選的，在步驟2.4中，基于步驟2.3獲得的特異性特征，用留一法(LOOCV)來(lái)驗(yàn)證和評(píng)估分類器的準(zhǔn)確率P：

$P=\frac{\underset{i=1}{\overset{\left|D\right|}{\Σ}}f\left(g\right(x_i),y_i)}{\left|D\right|}$

其中，P表示分類準(zhǔn)確率，D是由有限個(gè)以(x_i,y_i)形式表示的樣本組合集合x_i是樣本中除y_i以外的屬性列表，y_i表示樣本中類別標(biāo)號(hào)的屬性，g表示分類器模型函數(shù)，輸出結(jié)果為該模型的預(yù)測(cè)結(jié)果，f(g(x_i),y_i)為判別函數(shù)，當(dāng)g(x_i)與y_i相等時(shí)，輸出1，否則，輸出0。

由上述對(duì)本發(fā)明的描述可知，與現(xiàn)有技術(shù)相比，本發(fā)明具有如下有益效果：

1、本發(fā)明基于高通量測(cè)序數(shù)據(jù)，通過變尺度tuple頻度，不但能夠清楚的區(qū)分微生物的群落類別，而且基于獲取的群落類別，能夠找出區(qū)分群落類別的特異性信息。

2、本發(fā)明使用的方法無(wú)需人工選擇馬爾科夫階次，能根據(jù)數(shù)據(jù)本身特點(diǎn)自動(dòng)地選擇馬爾科夫階次，而且馬爾科夫模型中對(duì)應(yīng)的上下文序列可以是連續(xù)的也可以是不連續(xù)的；

3、本發(fā)明對(duì)宏基因組數(shù)據(jù)的聚類效果明顯優(yōu)于定階次馬爾科夫模型的聚類效果；

4、本發(fā)明為了更好更完整的比較分析微生物群落的群落結(jié)構(gòu)和找出微生物群落的類別差異性，針對(duì)不同尺度k-tuple詞頻的微生物高通量測(cè)序數(shù)據(jù)，采用基于無(wú)監(jiān)督聚類和基于聚類得到的樣本類別進(jìn)行有監(jiān)督分類相結(jié)合的微生物群落比較分析方法，把微生物群落的比較分析從統(tǒng)計(jì)分布層面擴(kuò)展到物種和基因分析層面。

5、本發(fā)明為了更好地刻畫微生物群落間的特異性信息，基于無(wú)監(jiān)督聚類獲得的樣本類別，我們首次將≥30bp的長(zhǎng)tuple作為特征，應(yīng)用于微生物高通量測(cè)序數(shù)據(jù)比較分析協(xié)議的有監(jiān)督的樣本分類方法中，用來(lái)找出區(qū)分樣本類別的特異性tuple序列特征。實(shí)例實(shí)驗(yàn)表明，當(dāng)k-tuple長(zhǎng)度等于40bp時(shí)最能代表兩類數(shù)據(jù)存在的差異性。

附圖說(shuō)明

圖1為采用固定階次馬爾科夫模型方法進(jìn)行聚類的結(jié)果；

圖2為插值上下文馬爾科夫模型方法進(jìn)行聚類的結(jié)果。

具體實(shí)施方式

以下通過具體實(shí)施方式對(duì)本發(fā)明作進(jìn)一步的描述。

本發(fā)明提供一種基于不同尺度tuple詞頻的微生物高通道量測(cè)序數(shù)據(jù)分析協(xié)議?；?-10bp的短tuple高通量測(cè)序數(shù)據(jù)，應(yīng)用所述插值上下文馬爾科夫模型進(jìn)行建模微生物群落的背景基因組，來(lái)比較宏基因組樣本，得出宏基因組樣本的類別信息。并且找出樣本分類的特異性特征，所述方法包括以下步驟：

步驟1：獲取宏基因組樣本的2-10bp的短tuple高通量測(cè)序數(shù)據(jù)，采用插值上下文馬爾科夫模型進(jìn)行建模微生物群落的背景基因組，再采用無(wú)監(jiān)督的聚類方法來(lái)比較宏基因組樣本，得出宏基于組樣本的類別信息。具體包括如下

步驟1.1:獲取宏基因組樣本高通量測(cè)序數(shù)據(jù)，生成該宏基因組樣本的短tuple特征頻度向量，對(duì)每個(gè)宏基因組樣本中出現(xiàn)的長(zhǎng)度為2-10bp的tuple的頻度進(jìn)行統(tǒng)計(jì)，并生成相應(yīng)宏基因組樣本的頻度向量；其中tuple特征定義為為宏基因組樣本中可能出現(xiàn)的字符串組合，并選擇長(zhǎng)度為2-10bp的字符串組合作為所述tuple特征。

步驟1.2:基于插值上下文馬爾科夫模型對(duì)微生物群落的背景基因組進(jìn)行建模，估計(jì)頻度向量中每一個(gè)tuple的馬爾科夫概率；計(jì)算tuple的馬爾科夫概率具體方法如下：

步驟1.2.1：基于宏基因組樣本的頻度向量的互信息量構(gòu)建上下文序列樹；。該步驟中，基于樣本的頻度向量，將k長(zhǎng)度tuple中的每一列的字符放在一個(gè)向量中，形成A₁,A₂,…,A_k k個(gè)向量。分別計(jì)算前面k-1個(gè)向量與最后一個(gè)向量(即當(dāng)前狀態(tài)向量)之間的互信息量，互信息量的公式如下：

$I(A_w,B)=\underset{i}{\Σ}\underset{j}{\Σ}P(a_i,b_j)*log\frac{P\left(a_i\right)*P\left(b_j\right)}{P(a_i,b_j)}---\left(1\right)$

其中,w＝1,2,…,k-1；B＝A_k；a_i,b_j表示向量A,B中的變量；P(a_i,b_j)表示a_i,b_j在對(duì)應(yīng)向量中同時(shí)出現(xiàn)的聯(lián)合概率；P(a_i)表示a_i在對(duì)應(yīng)向量中出現(xiàn)的概率。

找出與B互信息量最大的向量A_w，將此向量對(duì)應(yīng)的下標(biāo)位置作為上下文序列樹的頂點(diǎn)；然后按照該向量中出現(xiàn)的四種不同字符(A,C,G,T)將所有tuple分成四組；最后對(duì)四組tuple向量矩陣，按照公式(1)中的互信息量公式分別計(jì)算互信息量，找出每組中與Y互信息量最大的向量A_s，其中s＝1,2,…,w-1,w+1,…,k-1，將此向量對(duì)應(yīng)的下標(biāo)位置作為上下文序列樹的對(duì)應(yīng)葉子(A,C,G,T)的子節(jié)點(diǎn)；依次繼續(xù)下去，直到找到最后一個(gè)跟當(dāng)前狀態(tài)向量關(guān)聯(lián)性最大的向量，整個(gè)上下文序列樹構(gòu)建完畢。對(duì)于上述構(gòu)建的上下文序列樹，每條枝都對(duì)應(yīng)一個(gè)tuple，按照與當(dāng)前狀態(tài)關(guān)聯(lián)性程度從大到小，把這個(gè)tuple中的堿基位置自上而下排列,儲(chǔ)存在這條枝的節(jié)點(diǎn)中。

步驟1.2.2：基于上下文序列樹計(jì)算每個(gè)tuple的馬爾科夫概率。每個(gè)tuple的馬爾科夫概率公式如下：

P(c₁c₂…c_k)＝P_ICM(c₁)P_ICM(c₂|c₁)…P_ICM(c_k|c₁c₂…c_k-1) (2)

其中，c₁c₂…c_k表示k-tuple序列，P_ICM(c_k|c₁c₂…c_k-1)表示上下文序列c₁c₂…c_k-1轉(zhuǎn)移到當(dāng)前狀態(tài)c_k的ICM轉(zhuǎn)移概率。

對(duì)于計(jì)算上述每個(gè)ICM轉(zhuǎn)移概率，針對(duì)上述k-tuple序列c₁c₂…c_k，應(yīng)用ICM馬爾科夫模型構(gòu)建的上下文序列樹中找出與當(dāng)前狀態(tài)c_k關(guān)聯(lián)性程度從大到小排序的重要位置，重新構(gòu)建其上下文序列。例如，要構(gòu)建r階馬爾科夫模型(r≤k-1)，從上下文序列樹中找出與當(dāng)前狀態(tài)b_k關(guān)聯(lián)性程度從大到小排序的重要位置對(duì)應(yīng)的狀態(tài)分別為c₃,c₄,…,c_r，組成插值上下文序列M_r,那么構(gòu)建ICM的概率模型如下：

P_ICM(c_k|c₁c₂c₃…c_k-1)＝P_ICM,r(c_k|M_r) (3)

P_ICM,r(c_k|M_r)＝λ_r*P(c_k|M_r)+(1-λ_r)*P_{ICM，r-1}(c_k|M_r-1) (4)

$P\left(c_k\right|M_r)=\underset{a\∈A,C,G,T}{\Σ}\frac{N(M_r,A_k=c_k)}{N(M_r,A_k=a)}---\left(5\right)$

其中，λ_m表示m階馬爾科夫模型概率所占的權(quán)重系數(shù)；N(M_r,A_k＝c_k)表示所有的k-tuple中插值上下文序列是M_r，第k個(gè)位置是c_k的所有tuple的頻度之和。對(duì)于上述馬爾科夫模型概率所占的權(quán)重系數(shù)的計(jì)算公式為：

其中所述C表示樣本閾值，它由赤池信息量準(zhǔn)則AIC^R(C)確定，具體公式如下：

AIC^R(C)＝-2λ(S；M_k)+2|M_IMM,k,C| (7)

其中，λ(S；M_k)表示樣本S的偽似然度，計(jì)算公式為：

$\mathrm{\λ}(S;M_k)=\underset{c_1,\mathrm{...}{c}_k\∈(A,C,G,T)}{\Σ}N(c_1,\mathrm{...}{c}_{k-1}){\mathrm{logP}}_{ICM}\left(c_k\right|c_1,\mathrm{...}{c}_{k-1})---\left(8\right)$

|M_IMM,k,C|表示模型自由參數(shù)的個(gè)數(shù)。當(dāng)AIC^R(C)值最小時(shí)算出的C值作為樣本的閾值。

所述q表示兩個(gè)字符串之間差異度的卡方檢驗(yàn)值，其計(jì)算原理如下：

${\mathrm{\Δ}}_r\left(M_r\right)=\underset{a\∈A,C,G,T}{\Σ}\frac{{\left(N\right(M_r,a)-E(M_r,a\left)\right)}^2}{E(M_r,a)}---\left(9\right)$

E(M_r,a)＝N(M_r)P_ICM,r(a|M_r) (10)

其中，N(M_r,a)，E(M_r,a)分別表示字符串頻度的實(shí)際值和理論值。將q＝Δ_r(M_r)作為卡方值，自由度為3，作為卡方檢驗(yàn)的指標(biāo)參數(shù)。

步驟1.3：計(jì)算各個(gè)樣本頻度向量間的相異度距離，生成一個(gè)宏基因組樣本間的相異度矩陣。應(yīng)用不同的相異度距離度量方法計(jì)算各個(gè)宏基因組樣本頻度向量間的相異度距離，所用到的相異度距離度量方法包括和計(jì)算公式如下：

$D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}---\left(11\right)$

$d_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{Y,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}})---\left(12\right)$

$D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{n_X{p}_{X,i}}\sqrt{n_Y{p}_{Y,i}}}---\left(13\right)$

$d_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{X,i}^2}{n_X{p}_{X,i}}}\sqrt{\underset{i=1}{\overset{4^k}{\Σ}}\frac{{\tilde{C}}_{Y,i}^2}{n_Y{p}_{Y,i}}}})---\left(14\right)$

其中，和都是計(jì)算兩個(gè)樣本之間的相異度的一種距離度量方法；表示樣本X的頻度向量，表示樣本Y的頻度向量；稱為中心化過程；C_X,i和C_Y,i分別表示X和Y樣本中第i個(gè)tuple出現(xiàn)的頻度；n_X表示樣本X中tuple個(gè)數(shù)的總和，n_Y表示樣本Y中tuple個(gè)數(shù)的總和；p_X,i和p_Y,i分別表示在插值上下文馬爾科夫背景模型下，樣本X中第i個(gè)tuple的馬爾克夫概率和樣本Y中第i個(gè)tuple的馬爾科夫概率。

如果一個(gè)數(shù)據(jù)集中有n個(gè)樣本，根據(jù)相異度距離度量公式計(jì)算出的兩兩樣本間的相異度，生成一個(gè)n*n維的相異度距離矩陣，這個(gè)矩陣定義如下：

N(n,n)＝(d(x,y))_n×n,d(x,y)＝d(y,x),d(x,x)＝0 (15)

其中，d(x,y)是兩個(gè)宏基因組樣本的相異度距離，如果不同樣本間的距離越小，那么d(x,y)的值就越?。籨(x,x)表示相同樣本間的距離為0。

步驟1.4：根據(jù)n*n相異度矩陣生成一個(gè)聚類樹。由此判斷宏基因組樣本與樣本之間的關(guān)系。在相異度矩陣的基礎(chǔ)上，根據(jù)非加權(quán)平均法層次聚類算法計(jì)算出兩個(gè)群組的相異度距離。其定義如下：

$d(C_i,C_j)=\frac{1}{\left|C_i\right|\·\left|C_j\right|}\underset{x\∈C_i}{\Σ}\underset{y\∈C_j}{\Σ}d(x,y)---\left(16\right)$

d(x,y)是兩個(gè)樣本的相異度距離，|C_i|和|C_j|表示兩個(gè)群組的大小，也就是群組中樣本的個(gè)數(shù)，i,j＝1,2,…,n。由兩兩群組的相異度距離得到聚類樹，從聚類樹中可以直觀的看出群落中各個(gè)樣本之間的結(jié)構(gòu)關(guān)系，得出樣本之間的類別信息。

步驟2：基于步驟1)中聚類得出的類別信息，將≥30bp的長(zhǎng)tuple作為特征，采用有監(jiān)督的樣本分類方法找出描述樣本類別的特異性特征tuple序列。具體包括以下子步驟：

步驟2.1：對(duì)宏基因組樣本中出現(xiàn)的長(zhǎng)度為40bp的長(zhǎng)tuple的頻度進(jìn)行統(tǒng)計(jì)，并生成相應(yīng)樣本的頻度向量。該步驟可參照步驟1.1，只是在統(tǒng)計(jì)過程中將tuple的長(zhǎng)度變長(zhǎng)為40bp。

步驟2.2：對(duì)每個(gè)宏基因組樣本的tuple頻度向量進(jìn)行并行處理，生成所有樣本的長(zhǎng)tuple頻度向量矩陣，然后過濾掉冗余的特征。將需要分類樣本的tuple頻度向量合并在一起，生成一個(gè)tuple頻度向量矩陣A，A表示為M×N的頻度矩陣，其中N表示樣本數(shù)量，M表示特征維度。

步驟2.3：基于步驟1中獲取的樣本類別信息，應(yīng)用過濾好的樣本特征對(duì)樣本進(jìn)行有監(jiān)督分類，找到對(duì)分類效果具有很強(qiáng)辨識(shí)性的特異性tuple特征。具體的：基于步驟1中獲取的樣本類別信息，選定訓(xùn)練集和測(cè)試集樣本，在訓(xùn)練集中選定當(dāng)前類別和目標(biāo)類別。然后根據(jù)對(duì)稱不確定性大于設(shè)定的閾值時(shí)，把冗余的tuple序列特征過濾移除。得到一些類別特異性候選特征。對(duì)稱不確定性定義如下：

其中，N_X表示當(dāng)前類別組成的X樣本集中tuple特征出現(xiàn)的頻度；sum(N_X)表示由當(dāng)前類別組成的X樣本集中特征出現(xiàn)的頻度之和；sum(N_Y)表示目標(biāo)類別組成的Y樣本集中特征出現(xiàn)的頻度之和；n(X)和n(Y)分別表示X和Y樣本集中樣本的個(gè)數(shù)；θ表示X和Y之間對(duì)稱不確定性的閾值。

采用SVM分類器，對(duì)樣本進(jìn)行有監(jiān)督分類，找到能夠描述微生物群落內(nèi)部具有差異性的特異性特征。

步驟2.4：基于步驟2.3獲得的特異性特征，用留一法(LOOCV)來(lái)驗(yàn)證和評(píng)估分類器的準(zhǔn)確率P：

$P=\frac{\underset{i=1}{\overset{\left|D\right|}{\Σ}}f\left(g\right(x_i),y_i)}{\left|D\right|}---\left(18\right)$

其中，P表示分類準(zhǔn)確率，D是由有限個(gè)以(x_i,y_i)形式表示的樣本組合集合x_i是樣本中除y_i以外的屬性列表，y_i表示樣本中類別標(biāo)號(hào)的屬性，g表示分類器模型函數(shù)，輸出結(jié)果為該模型的預(yù)測(cè)結(jié)果，f(g(x_i),y_i)為判別函數(shù)，當(dāng)g(x_i)與y_i相等時(shí)，輸出1，否則，輸出0。

本發(fā)明針對(duì)由高通量測(cè)序獲得的宏基因組樣本，對(duì)微生物群落進(jìn)行比較和分析。接下來(lái)詳細(xì)描述本發(fā)明的方法的實(shí)施過程。雖然在以下內(nèi)容中展示了執(zhí)行步驟的邏輯過程，但在某些情況下，可以不同此處的順序執(zhí)行。

首先執(zhí)行步驟1中的步驟1.1,獲取宏基因組樣本的k-tuple頻度向量。k-tuple是指長(zhǎng)度為k的連續(xù)字符串。在本發(fā)明中，通過統(tǒng)計(jì)這些字符串在樣本中出現(xiàn)的頻度，并將這些頻度組合成一個(gè)k-tuple頻度向量，以此來(lái)代表整個(gè)樣本的特征。在本發(fā)明中，先選擇的tuple尺度標(biāo)準(zhǔn)為長(zhǎng)度為2-10bp的字符串作為k-tuple的tuple特征。

為了計(jì)算tuple序列特征的插值上下文馬爾科夫概率，在本實(shí)施例中需要執(zhí)行步驟1.2。在步驟1.2中，首先執(zhí)行步驟1.2.1：根據(jù)樣本的所有tuple建立一個(gè)上下文序列樹，在構(gòu)建的過程中需根據(jù)互信息量最大的準(zhǔn)則，依次找到與當(dāng)前狀態(tài)關(guān)聯(lián)性最大的點(diǎn)，然后添加到上下文序列的節(jié)點(diǎn)中，每個(gè)子節(jié)點(diǎn)再按A,C,G,T作為葉子向下分枝，在每個(gè)分枝下，再依照互信息量最大的準(zhǔn)則，向下添加子節(jié)點(diǎn)。上下文序列樹中，父節(jié)點(diǎn)表示的tuple字符位置包含在子節(jié)點(diǎn)表示的tuple字符位置中。

在步驟1.2中，當(dāng)整個(gè)上下文序列樹構(gòu)建好之后，這時(shí)候每個(gè)tuple的上下文序列按照和當(dāng)前狀態(tài)的關(guān)聯(lián)性大小從大到小都會(huì)依次存儲(chǔ)在樹的每個(gè)節(jié)點(diǎn)中。進(jìn)而進(jìn)行步驟1.2中的1.2.2子步驟操作。由構(gòu)建上下文序列樹的過程可知，原先tuple的上下文序列轉(zhuǎn)移到下一個(gè)狀態(tài)的概率可以用經(jīng)過排好序的上下文序列轉(zhuǎn)移到下一個(gè)狀態(tài)的轉(zhuǎn)移概率替代。根據(jù)這一原則，可以估算出每一個(gè)tuple的馬爾科夫概率。

為了計(jì)算實(shí)施例中k-tuple向量之間的距離，接下來(lái)實(shí)施步驟1.3。對(duì)k-tuple向量分別采取和的相異度方法計(jì)算距離。其中距離度量方法中用到的tuple的插值上下文馬爾科夫概率，在步驟1.2中已求得。

實(shí)施例步驟1.3可以得到一個(gè)相異度矩陣，由這個(gè)相異度矩陣進(jìn)行步驟1.4，即進(jìn)行無(wú)監(jiān)督的層次聚類分析，最終可得到一個(gè)聚類樹。通過觀察聚類樹，可以判斷聚類情況的好壞，找出樣本的類別信息。

實(shí)施例步驟2.1和步驟1.1類似。通過實(shí)施例中步驟2.2對(duì)得到的tuple頻度向量進(jìn)行合并，得到tuple頻度向量矩陣。然后實(shí)施該步驟中的特征過濾，將樣本中的tuple特征出現(xiàn)頻度不為零的特征頻度歸一化為1，然后利用對(duì)稱不確定性來(lái)計(jì)算當(dāng)前類別和目標(biāo)類別的相關(guān)性熵，將相關(guān)性熵大于某一設(shè)定的閾值的特征留下來(lái)，這些特征就是類別特異性候選特征。

利用這些類別特異性候選特征對(duì)基于步驟1所獲得的類別信息進(jìn)行有監(jiān)督的樣本分類，需要執(zhí)行步驟2.3，應(yīng)用SVM分類器，選定訓(xùn)練集和測(cè)試集樣本，在訓(xùn)練集中選定當(dāng)前類別和目標(biāo)類別，通過學(xué)習(xí)建立分類模型，找出能分開當(dāng)前類別和目標(biāo)類別的特異性tuple特征。為刻畫不同類別的微生物群落具體差異以及尋找生物標(biāo)記物提供重要的參考信息。最后執(zhí)行步驟2.4，應(yīng)用留一法對(duì)分類器分類準(zhǔn)確性進(jìn)行評(píng)估。為刻畫不同類別的微生物群落具體差異以及尋找生物標(biāo)記物提供重要的參考信息。

我們選取24個(gè)人體皮膚微生物群落樣本(NCBI基因數(shù)據(jù)庫(kù)http://www.ncbi.nlm.nih.gov/)進(jìn)行無(wú)監(jiān)督聚類實(shí)驗(yàn)，分別采用了固定階次馬爾科夫模型和插值上下文馬爾科夫模型的方法，結(jié)果顯示樣本按人體左右兩個(gè)不同位置分別進(jìn)行聚類，而且插值上下文馬爾科夫模型的方法(參見圖1，)得出的結(jié)果要好于固定階次馬爾科夫模型(參見圖2，)。

在有監(jiān)督聚類過程中，選取了99個(gè)健康成年人和25個(gè)患有腸胃炎(IBD)病人的糞便樣本(Qin,J.,et al.,A human gut microbial gene catalogue established by metagenomic sequencing.Nature,2010.464(7285):p.59-65.)，將25個(gè)IBD病人樣本和25個(gè)健康人樣本作為訓(xùn)練集，用SVM分類器建立分類模型；將剩下的74個(gè)健康人樣本作為測(cè)試集，并進(jìn)行LOOCV實(shí)驗(yàn)來(lái)評(píng)估分類器性能，最終結(jié)果顯示，分別將40-tuple和k-tuple(k＝2-10)作為特征，評(píng)估分類器分類性能，基于40-tuple作為特征構(gòu)建的分類器只要用一個(gè)特征就能平均獲得100％的準(zhǔn)確率，而基于2-10tuple作為特征時(shí)，其最好的分類準(zhǔn)確率為88％(k＝7)，且需要200個(gè)特征。實(shí)驗(yàn)表明，長(zhǎng)tuple序列比短tuple包含更多的顯著的類別特異性信息。一個(gè)最直觀的表現(xiàn)就是分類性能的準(zhǔn)確性。參見表1

表1

上述僅為本發(fā)明的具體實(shí)施方式，但本發(fā)明的設(shè)計(jì)構(gòu)思并不局限于此，凡利用此構(gòu)思對(duì)本發(fā)明進(jìn)行非實(shí)質(zhì)性的改動(dòng)，均應(yīng)屬于侵犯本發(fā)明保護(hù)范圍的行為。