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