用于疾病亞型問(wèn)題的基于網(wǎng)絡(luò)的聚類方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明是關(guān)于逆向研究疾病亞型領(lǐng)域，特別涉及用于疾病亞型問(wèn)題的基于網(wǎng)絡(luò)的 聚類方法。
【背景技術(shù)】
[0002] 對(duì)于由基因變異導(dǎo)致的疾病的研究一直是一個(gè)非常熱門的議題。在這些疾病中， 很多疾病都有不同的亞型。所謂亞型（subtype)，是同一個(gè)疾病下的不同的種型。它們可能 由不同的病因引起，并且有不同的臨床表征。例如HIV有1型和2型，腫瘤亞型有卵巢癌、 肺癌、子宮癌等。
[0003] 對(duì)于疾病亞型的很多研究，現(xiàn)階段還是集中在對(duì)于其病理的正向研究。而運(yùn)用逆 向工程技術(shù)（reverse engineering)，逆向研究疾病亞型也逐漸成為一個(gè)熱門的話題。"逆 向工程技術(shù)"是一個(gè)研究主體系統(tǒng)的過(guò)程。它通過(guò)研究主體系統(tǒng)來(lái)鑒定系統(tǒng)的各個(gè)成分 以及它們之間的相互關(guān)聯(lián)，并用另一種方式在更加抽象，更加上層的水平上對(duì)該系統(tǒng)進(jìn)行 代表。逆向工程技術(shù)在疾病亞型鑒定與分類方面研究上的一個(gè)非常重要的應(yīng)用，就是通 過(guò)已有的臨床信息，例如基因表達(dá)數(shù)據(jù)（gene expression data)等，運(yùn)用包括聚類分析 (cluster analysis)在內(nèi)的各種手段，反向研究并鑒定疾病的亞型。然而，由于基因的表達(dá) 之間并非是獨(dú)立的，而是會(huì)通過(guò)各種關(guān)系相互作用。因此，最終獲取的基因表達(dá)數(shù)據(jù)，也應(yīng) 該是相互關(guān)聯(lián)的。而在以往的研究中，人們只是使用傳統(tǒng)的聚類方法，并沒(méi)有考慮這種基因 表達(dá)數(shù)據(jù)間的相互關(guān)聯(lián)。因此，將有關(guān)基因的作用關(guān)系的信息加入聚類分析中是一個(gè)自然、 新穎的想法并且值得一試。

【發(fā)明內(nèi)容】

[0004] 本發(fā)明的主要目的在于克服現(xiàn)有技術(shù)中的不足，提供能更好的將疾病亞型進(jìn)行分 類，更好的還原真實(shí)的疾病亞型的基于網(wǎng)絡(luò)的聚類方法。為解決上述技術(shù)問(wèn)題，本發(fā)明的解 決方案是：
[0005] 提供用于疾病亞型問(wèn)題的基于網(wǎng)絡(luò)的聚類方法，具體包括下述過(guò)程：
[0006] (1)獲得O-G矩陣以及基因調(diào)控網(wǎng)絡(luò)；
[0007] (2)選取適用于具體問(wèn)題的基于網(wǎng)絡(luò)的距離定義，構(gòu)建距離矩陣；
[0008] (3)運(yùn)用k-medoids算法對(duì)O-G矩陣進(jìn)行聚類分析；聚類時(shí)距離的選擇用基于網(wǎng) 絡(luò)的距離；
[0009] (4)得出最終關(guān)于疾病亞型的分類；
[0010] 所述過(guò)程（1)具體包括下述步驟：
[0011] 步驟A :根據(jù)基因調(diào)控網(wǎng)絡(luò)（即基因-蛋白質(zhì)調(diào)控網(wǎng)絡(luò)，是一個(gè)細(xì)胞中DNA片段集 合通過(guò)相互間的各種非間接作用，比如RNA作用以及蛋白質(zhì)表達(dá)作用，來(lái)影響其mRNA以及 蛋白質(zhì)表達(dá)水平的相互關(guān)系）的特性（例如網(wǎng)絡(luò)的平均出度、入度等參數(shù)），構(gòu)建隨機(jī)的有 向圖來(lái)代表基因調(diào)控網(wǎng)絡(luò)G(V，E);其中每個(gè)頂點(diǎn)i e V代表基因i及其產(chǎn)生的IiiRNA1和蛋 白質(zhì)i (轉(zhuǎn)錄因子i);每條有向邊e]ie E代表著"轉(zhuǎn)錄因子j調(diào)控基因i的轉(zhuǎn)錄"這種調(diào) 控關(guān)系；
[0012] 步驟B :根據(jù)產(chǎn)生的基因調(diào)控網(wǎng)絡(luò)G(V，E)，對(duì)每個(gè)基因i建立激活函數(shù)仁（·），具 體建立方式為：
[0013] 對(duì)于任意的基因i e V，i = l，2，K，n，我們從G(V，E)中找出所有與i相鄰且以i 為有向邊終點(diǎn)的點(diǎn)，構(gòu)成影響因子集合{qpqytqj ;其中，Q1表示與i相鄰且以i為有向 邊終點(diǎn)的某基因中對(duì)基因i起影響作用的因子，q2表示與i相鄰且以i為有向邊終點(diǎn)的某 基因中對(duì)基因i起影響作用的因子，q ln表示與i相鄰且以i為有向邊終點(diǎn)的某基因中對(duì)基 因i起影響作用的因子，η表示基因調(diào)控網(wǎng)絡(luò)中基因的數(shù)量；
[0014] 確定解離常數(shù)Iclj，且Iclj從定義在[0. 01，1]區(qū)間上的均勻分布中選??；
[0015] 確定希爾系數(shù)II1 j，且II1 j服從[1，10]區(qū)間中的高斯分布函數(shù)，V2,4);
[0016] 確定相對(duì)活性a i，且α，人定義在[0, 1]區(qū)間上的均勻分布上采樣；
[0017] 步驟C :確定無(wú)噪聲動(dòng)態(tài)基因調(diào)控模型，即確定公式（2. 1)的各個(gè)參數(shù)；
[0018]
[0019] 式（2. 1)中，\表示基因 i的濃度；yi表示蛋白質(zhì)i的濃度；表示HiRNA1的濃 度變化率；表示蛋白質(zhì)i的濃度變化率；叫表示基因i的最大轉(zhuǎn)錄速率；r i表示mRNA i 的翻譯速率；Afy表示HiRNA1的降解速率；表示蛋白質(zhì)i的降解速率；A (·)表示基因 i的激活函數(shù)；
[0020] 確定公式（2. 1)中各個(gè)參數(shù)的具體方式為：mRNA的半衰期:/fKi以及蛋白質(zhì)的半衰 期:?產(chǎn)"（以分鐘為單位）從定義在[5,50]區(qū)間上的高斯分布謂:27.5,56.:25)上采樣；
[0021] 根據(jù)公式（2. 9)，獲得mRNA以及蛋白質(zhì)的降解速率，最大轉(zhuǎn)錄速率叫以及翻譯速 率A服從[0. 01，0. 011]區(qū)間上的均勻分布；
[0022]
[0023] 式（2. 9)中，表示HiRNA1的降解速率；λ"表示蛋白質(zhì)i的降解速率；mRNA的 半衰期及蛋白質(zhì)的半衰期?Γ 5"(以分鐘為單位）；
[0024] 步驟D :在獲得了基因調(diào)控網(wǎng)絡(luò)以及無(wú)噪聲動(dòng)態(tài)基因調(diào)控模型之后，選定mRNA濃 度χ(Χρ χ2, Κ，χη)以及蛋白質(zhì)濃度y(yp y2, K，yn)的初始值（可以令各個(gè)xjp y i服從[0， 1]區(qū)間上的均勻分布，并隨機(jī)選取作為初始值），然后求解公式（2.1)，得到最終的基因表 達(dá)數(shù)據(jù)；
[0025] 所述過(guò)程⑵具體是指：根據(jù)過(guò)程⑴所獲得的基因網(wǎng)絡(luò)的拓?fù)潢P(guān)系G(V，E)，定 義三種基于網(wǎng)絡(luò)的距離，用于比較1 1(111，112，1^111)與12(1 21，122，1(，1211)的差別；其中11(叉 11， x12, K，xln)、x2(x21，x22, K，x2n)分別表示兩個(gè)被試者 PjP P 2的 mRNA 濃度；
[0026] 令G(V，E)代表該基因調(diào)控網(wǎng)絡(luò)，其中每個(gè)頂點(diǎn)i e V代表基因i及其產(chǎn)生的KiRNA1 和蛋白質(zhì)i (轉(zhuǎn)錄因子i);它關(guān)聯(lián)的Xi表示該基因轉(zhuǎn)錄的HiRNAi濃度；令每條有向邊E 代表著"轉(zhuǎn)錄因子j調(diào)控基因i的轉(zhuǎn)錄"這種調(diào)控關(guān)系；記T 1表示與節(jié)點(diǎn)i相連的邊數(shù)（即 節(jié)點(diǎn)i的度），Ii表示節(jié)點(diǎn)i的入度，〇 i表示節(jié)點(diǎn)i的出度；
[0027] 其中，基于網(wǎng)絡(luò)的Jaccard距離定義為：
[0028]
公式（3. 10);
[0029] 其中，令G(V，E)代表該基因調(diào)控網(wǎng)絡(luò)，其中每個(gè)頂點(diǎn)i e V代表基因i及其產(chǎn)生 的mRNAjP蛋白質(zhì)i (轉(zhuǎn)錄因子i);它關(guān)聯(lián)的X ;表示該基因轉(zhuǎn)錄的mRNA ;濃度；T ;表示與節(jié) 點(diǎn)i相連的邊數(shù)（即節(jié)點(diǎn)i的度），Ii表示節(jié)點(diǎn)i的入度，0 i表示節(jié)點(diǎn)i的出度；X H指被試 者Pl的IiiRNA1濃度；X 21指被試者P2的mRNA i濃度；η表示基因調(diào)控網(wǎng)絡(luò)中基因的數(shù)量；
[0030] 基于網(wǎng)絡(luò)的Euclidean距離：
[0031] ；
[0032]
[0033] 其中，X11指被試者Pl的mRNA i濃度；X 21指被試者P2的mRNA i濃度；X i;指被試者 Pl的HiRNAj濃度；X 2j指被試者P2的mRNA j濃度；η表示基因調(diào)控網(wǎng)絡(luò)中基因的數(shù)量；
[0034] 基于網(wǎng)絡(luò)的Pearson距離：
[0035]
[0036] 其中，知指被試者Pl的mRNA i濃度；X 21指被試者P2的mRNA i濃度；η表示基因調(diào) 控網(wǎng)絡(luò)中基因的數(shù)量； CN 105160208 A 說(shuō)明書 4/7 頁(yè)
[0037] 1廣示節(jié)點(diǎn)i的入度；
這里的心指被試者Pi的 InRNA1濃度；這里的X 12指被試者Pi的mRNA 2濃度；
[0038] 所述過(guò)程⑶具體是指：將過(guò)程（2)中定義的距離引入聚類分析中，使用 k-medoids聚類分析方法，對(duì)過(guò)程（1)所獲得的基因表達(dá)數(shù)據(jù)進(jìn)行聚類；
[0039] 假設(shè)有η個(gè)被試者，我們將η個(gè)被試者劃分為k類，K-medoids聚類算法是，基于 網(wǎng)絡(luò)的Pearson距離具體的算法具體方法如下：
[0040] (a)從η個(gè)數(shù)據(jù)對(duì)象中任意選取k個(gè)數(shù)據(jù)對(duì)象作為medoids-聚類的中心，
[0041] (b)選定基于網(wǎng)絡(luò)的Person距離，即：
[0042]
[0043] 然后分別計(jì)算余下的數(shù)據(jù)對(duì)象到各個(gè)聚類中心的距離，并將余下的數(shù)據(jù)對(duì)象分配 到離自己最近的聚類中，最終得到k組劃分，G1, G2, ...，Gk;
[0044] (c)數(shù)據(jù)對(duì)象分配完成后，順序選取一個(gè)數(shù)據(jù)對(duì)象來(lái)代替原來(lái)的聚類中心，并計(jì)算 代替后的優(yōu)化目標(biāo)函數(shù)
[0045] 其中，d(X1, X2)定義如下：
[0046]
[0047] 同理定義(!(Xi, xj和丨其中，
為從X1, x2, . . .，xA選取的k 個(gè)聚類中心；表示Xj e G1;
[0048] 再選擇f最小的數(shù)據(jù)對(duì)象來(lái)代替聚類中心，這樣K個(gè)mediods就改變了；
[0049] (d)與前一次的聚類中心相比較，如果發(fā)生變化轉(zhuǎn)到方法（b)，如果不發(fā)生變化轉(zhuǎn) 到方法（e);
[0050] (e)將聚類的結(jié)果輸出；
[0051] 所述過(guò)程（4)具體是指：根據(jù)過(guò)程（3)的聚類結(jié)果，得出最終關(guān)于疾病亞型的分 類。
[0052] 與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：
[0053] 對(duì)于特定的基因網(wǎng)絡(luò)，基于網(wǎng)絡(luò)的聚類方法將有更好的組間相似性，更有效地還 原三種亞型。此外，當(dāng)有大量的基因需要測(cè)定其表達(dá)數(shù)據(jù)時(shí)，現(xiàn)有的方法可能無(wú)法同時(shí)對(duì)所 有的基因進(jìn)行精確的測(cè)量。此時(shí)，我們提出的"基于網(wǎng)絡(luò)的聚類"法使得我們通過(guò)優(yōu)先精確 測(cè)量信息基因的表達(dá)數(shù)據(jù)，并不會(huì)大大地削弱對(duì)于疾病亞型的鑒定效果。
【附圖說(shuō)明】
[0054] 圖1為本發(fā)明的操作流程圖。
【具體實(shí)施方式】
[0055] 下面結(jié)合附圖與【具體實(shí)施方式】對(duì)本發(fā)明作進(jìn)一步詳細(xì)描述：
[0056] 現(xiàn)在，我們假定一共有32位被試者P2.....P32，其中被試者P 2.....P8為正 常狀態(tài)病人，被試者P9、P1C.....P16為患有基因疾病亞型D i的病人，被試者P 17、P18.....P24 患有基因疾病亞型D2，被試者P25、P26.....P32患有基因疾病亞型D 3。D2、D3中的每一種 亞型都代表著某些基因表達(dá)的失常。為了模擬這個(gè)表達(dá)失常過(guò)程，對(duì)于某一種亞型，我們從 整個(gè)基因調(diào)控網(wǎng)絡(luò)中隨機(jī)選取一定的節(jié)點(diǎn)（也就是基因），對(duì)其最大轉(zhuǎn)錄速率Hl1進(jìn)行擾動(dòng)。 對(duì)于不同的亞型，我們選取不同的基因進(jìn)行擾動(dòng)。我們希望做的是通過(guò)對(duì)32位被試者最后 的mRNA濃度向量進(jìn)行聚類分析，試圖分出對(duì)照組與三種疾病亞型。
[0057] 步驟A :我們根據(jù)基因調(diào)控網(wǎng)絡(luò)的某些特性（例如網(wǎng)絡(luò)的平均出度、入度等參數(shù)） 來(lái)構(gòu)建隨機(jī)的有向圖來(lái)代表基因調(diào)控網(wǎng)絡(luò)構(gòu)建基因調(diào)控網(wǎng)絡(luò)。假設(shè)我們要產(chǎn)生由η個(gè)基因 組成的基因調(diào)控網(wǎng)絡(luò)，根據(jù)基因調(diào)控網(wǎng)絡(luò)的特性，我們將產(chǎn)生一張平均入度為2,且分布滿 足冪定理分布（power law distribution)的隨機(jī)有向網(wǎng)絡(luò)G(V，Ε)，其中I V| = η。此外， 圖中不允許有自環(huán)的出現(xiàn)。
[0058] 步驟B :根據(jù)我們產(chǎn)生的基因調(diào)控網(wǎng)絡(luò)G(V，E)