[發明專利]一種識別癌癥驅動通路的方法

有效

• 申請號: 202011185104.6
• 申請日: 2020-10-30
• 公開（公告）號: CN112270952B
• 公開（公告）日: 2022-04-05
• 發明（設計）人: 朱凱;吳璟莉;李高仕
• 申請（專利權）人: 廣西師范大學
• 主分類號: G16B20/10
• 分類號: G16B20/10;G16B20/50
• 代理公司: 桂林市華杰專利商標事務所有限責任公司 45112
• 代理人: 楊雪梅

1.一種識別癌癥驅動通路的方法，其特征在于，包括如下步驟：

1)構造加權的非二進制突變矩陣：

現有某個癌癥的體細胞突變矩陣拷貝數變異矩陣和基因表達矩陣在體細胞突變矩陣拷貝數變異矩陣和基因表達矩陣三個矩陣中行表示該癌癥的相同樣本集p，列分別表示基因集G_S、G_C和G_E，在矩陣中，s_ij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|G_S|)，i樣本中j基因突變，s_ij值為1，反之值為0；矩陣中每個元素c_ij∈{-2,-1,0,1,2}(i＝1,2,…,|p|,j＝1,2,…,|G_C|)，表示i樣本中j基因拷貝數變異值；在矩陣中e_ij∈R(i＝1,2,…,|p|,j＝1,2,…,|G_E|)，表示i樣本中j基因表達量；令矩陣中的基因集為G_A＝G_S∪G_C，樣本集為p，令a_ij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|G_A|)，其中為突變矩陣，當s_ij取值為1或i樣本中j基因處于統計顯著變異區域時，a_ij值為1，反之值為0，為了進一步整合突變矩陣和表達矩陣在突變矩陣和表達矩陣中取基因集G＝G_A∩G_E，重新得到兩個矩陣A_|p|×|G|和E_|p|×|G|，對于基因表達數據，存在正常樣本表達矩陣N_|n|×|G|，n表示正常樣本，在矩陣N_|n|×|G|中，n_ij∈R(i＝1,2,…,|p|,j＝1,2,…,|G|)，表示i樣本中j基因表達量，令差異倍數矩陣D_|p|×|G|，d_ij∈R(i＝1,2,…,|p|,j＝1,2,…,|G|)，表示i樣本中j基因表達量相比j基因在正常樣本中表達量的差異倍數用表示，其中則d_ij值為否則d_ij值為0，處理好差異倍數矩陣D_|p|×|G|，進一步對突變矩陣A_|p|×|G|進行加權處理，整合成加權突變矩陣，對于A_|p|×|G|，如果a_ij＝1，并且d_ij≥λ₁，則a_ij＝1.5，如果a_ij＝0，并且d_ij≥λ₂，則a_ij＝(2·l)^-1·d_ij，其中λ₁和λ₂是截取差異倍數的閾值，l是j基因對應所有樣本中差異倍數的最大值，針對突變基因，λ₁取較低值，使a_ij∈{1,1.5}，以提高該突變基因的突變可信值；針對不突變基因，λ₂取較高值，使a_ij∈[0,0.5]，以提高該不突變基因的突變可信值，使其可能成為潛在基因，經過加權重新得到加權突變矩陣A_|p|×|G|，a_ij∈[0,1.5](i＝1,2,…,|p|,j＝1,2,…,|G|)；

2)設定識別模型：

針對加權突變矩陣A_|p|×|G|，基于高覆蓋和高互斥兩個特性，重新構建新的整合模型，假設M_|p|×k為矩陣A_|p|×|G|的任一子矩陣，令Γ(m)＝{m_i|m_i＝max{a_im|m∈M},i＝1,2,…,|p|}記錄矩陣M_|p|×k每行中最大權值，令矩陣M_|p|×k的覆蓋度對于矩陣M_|p|×k中一行的互斥度，考慮這一行的離散程度，用變異系數計算每行的互斥度，每行互斥度之和為整個M_|p|×k的互斥度，具體表示如公式(1)所示：

其中當趨近于0值時，對于變異系數值影響很大，所以如果M_|p|×k中一行的最大權值m_i≤0.5，則令該行的互斥度為放縮該行互斥度，避免該行互斥度過高對通路識別造成影響，根據和公式(1)，對整合數據后的最大權重子矩陣問題重新定義模型：給定突變加權矩陣A_|p|×|G|和正整數k(k＜|G|)，在矩陣A_|p|×|G|中確定矩陣M_|p|×k，使函數值W(M)最大，如公式(2)所示：

W(M)＝α(M)+ω(M) (2)，

其中α(M)表示矩陣M_|p|×k的覆蓋度，α(M)由矩陣M_|p|×k中每行最大突變權值相加所得，α(M)越大表示覆蓋樣本越多，并且突變可信值也越大；ω(M)表示矩陣M_|p|×k的互斥度，ω(M)由M中每行變異系數值相加所得，變異系數越大，則離散程度越高，其互斥度就越大；

3)設定適應度函數：

每個染色體對應一個問題解，因此需要對該解進行評估，給定染色體X_i(i＝1,2,…P)，P為種群大小，適應度函數Fitness(X_i)的定義如公式(3)所示：

其中，表示染色體X_i對應的子矩陣；

4)設定交叉算子：

按照染色體適應度從大到小，給染色體X_i一個排名R_i，則每個染色體被選中的概率如公式(4)所示：

為了保證染色體的可行性，采用輪盤賭隨機從父種群選取兩條染色體，重復基因分別給子代的兩條染色體，剩余基因放在一個集合，采用均勻交叉的方式，對于剩余基因的集合中連續的每對基因，隨機生成一個二值數據，如果該二值數據為1，則這對基因的第一個基因放入第一條子染色體，第二個基因放入第二條子染色體，反之，第一個基因放入第二條子染色體，第二個基因放入第一條子染色體，經過一次交叉，生成兩個子染色體；

5)設定變異算子：

給定一條子染色體X＝{x₁,x₂,…,x_k}(x_i＝1,2,…,|G|)，確定候選基因集合從子染色體中隨機刪除一個基因，得到基因集X′，將H_X中基因順序打亂，遍歷前個基因，選出基因g，使適應值Fitness(M_X′∪{g})最大，對應于子矩陣M_X′∪{g}的基因集X′∪{g}為新子染色體，即X＝X′∪{g}；

6)設定合作策略：

采用種群間相互合作策略，在種群交叉、變異和選擇操作后，比較兩個種群適應度最好的染色體和對方適應度最差的染色體，如果最好適應度高于對方最差適應度，則將該條染色體替換對方適應度最差的染色體；

7)設定參數：

輸入加權的非二進制突變矩陣A_|p|×|G|和公式(2)中的模型，參數k用于限制找到的驅動通路大小，然后輸入CGA-MWS相關參數：種群規模P、變異概率P_m、最大演化代數maxstep、最優值保持恒定的閾值maxt；

8)構造初始種群：

染色體編碼采用十進制編碼方式，一個解為k個基因構成的集合，即X＝{x₁,x₂,…,x_k}(x_i＝1,2,…,|G|)，將|G|個基因順序隨機打亂，然后取前k個基因構成初始染色體，生成兩個初始種群pop₀和pop₁，每個種群大小為P/2，計算兩個種群各自染色體的適應值，將pop₀和pop₁中各自最優的染色體相比較，保存最好的個體到變量best中，初始迭代次數step＝0，最優值保持恒定的代數t＝0；

9)執行迭代操作：

(1)若step＞maxstep或t＞maxt，轉入步驟9)的(4)，得到大小為k的驅動通路，否則轉入步驟9)的(2)；

(2)種群pop₀和pop₁各自通過基于排名的概率，采用輪盤賭隨機選兩條父代染色體，通過交叉算子進行交叉，生成兩條子染色體，并各自放入子種群pop₀′和pop₁′中，重復P/4次，針對子種群pop₀′和pop₁′，對每條染色體隨機一個突變概率P_m′，如果P_m′＜P_m，則對該條染色體進行變異操作，將變異操作中得到的適應值最高的染色體替換該條染色體，將pop₀和pop₀′里所有染色體按適應值從高到低排序，取前P/2條染色體放入下一代種群pop_step+1中，pop₁和pop₁′進行相同操作得到下一代種群pop_step+2；

(3)對種群pop_step+1和pop_step+2進行合作策略操作，比較兩個種群的最優適應值，取兩者中適應值最高的染色體，若該染色體適應值大于best染色體的適應值，則更新best染色體，t＝0；否則t＝t+1，step＝step+1，返回步驟9)的(1)；

(4)將best染色體轉換為基因集，由此得到子矩陣M，并子矩陣M其輸出，輸出的子矩陣M即為大小為k的驅動通路。