季姜帥,裴頌文,2
1(上海理工大學(xué) 光電信息與計算機工程學(xué)院,上海 200093)
2(復(fù)旦大學(xué) 上海市數(shù)據(jù)科學(xué)重點實驗室,上海 200433)
E-mail:swpei@usst.edu.cn
癌癥是威脅人類健康的主要疾病之一,常見的診斷與治療方法有CT、放療和化療等,雖然有一定的效果,但患者在接受治療過程中受到的輻射會產(chǎn)生一定的副作用.近年來,基因芯片技術(shù)飛速的發(fā)展,帶來了許多蘊含巨大科學(xué)價值的數(shù)據(jù),從基因?qū)用嫔贤诰虬┌Y的致病基因并提出針對性的治療方法是一大趨勢.聚類分析可以有效識別基因大數(shù)據(jù)背后的潛在結(jié)構(gòu)和價值.如利用聚類算法對腫瘤相關(guān)基因表達數(shù)據(jù)進行分型,可以識別出在臨床、病理等方面表現(xiàn)差異的腫瘤亞型,不同的亞型對應(yīng)著不同的治療方法,為臨床科學(xué)研究提供新思路和新方法[1].
聚類分析是一種基于數(shù)據(jù)點之間相似度的無監(jiān)督學(xué)習(xí)方法.目前,聚類分析已廣泛應(yīng)用于圖像處理、模式識別、電力系統(tǒng)、統(tǒng)計學(xué)以及生物信息學(xué)等領(lǐng)域.陳葉旺等人致力于密度聚類算法研究,基于有效半徑和密度熱量,對DBSCAN、MeanShift、DPeak等進行改進,提出了DHeat算法,能自適應(yīng)確定掃描半徑[2];之后進一步使用快速近似KNN算法檢測所有類型相同的核心塊、非核心塊和噪聲塊,設(shè)計了一種快速合并核心塊并將非核心點分配到合適簇以加速聚類過程[3].為了解決基于質(zhì)心的聚類方法無法識別非各向同性的極其復(fù)雜數(shù)據(jù)的問題.陳葉旺等人提出了基于查找密度核的混合分散式算法DCore[4].Qaddoura等人使用基于距離矢量的選擇操作迭代地考慮一個最近的點,提出了具有索引比的最近點算法NPIR,用于發(fā)現(xiàn)非球形簇的聚類[5].黃曉輝等人利用特征加權(quán)方法,提出了KICIC算法,通過在子空間內(nèi)最大化簇中心與其他簇的距離進行聚類[6].
遺傳算法采用自然進化模型,通過對進化過程中的種群,反復(fù)進行選擇、交叉和變異等操作,模擬自然界中種群演變的過程,直到滿足一定的條件才停止[7].雖然在整個進化過程中遺傳操作是隨機的,但它可以有效地根據(jù)已有的信息自適應(yīng)地控制搜索的方向,找出最優(yōu)解.由于遺傳算法本身的結(jié)構(gòu)決定了它可以很好地應(yīng)用在基因表達量數(shù)據(jù)中.
近年來,遺傳算法在生物醫(yī)學(xué)中的應(yīng)用研究越來越廣泛.Toubiana等人通過逐漸增加性狀與基因模塊的子集之間的相關(guān)性來優(yōu)化性狀與基因模塊的關(guān)系,提出了基于遺傳算法的隨機優(yōu)化算法[8].劉鵬等人為了解決醫(yī)學(xué)圖像中的噪聲問題,提出了基于遺傳算法的網(wǎng)絡(luò)進化方法,通過基于經(jīng)驗的貪婪探索策略和轉(zhuǎn)移學(xué)習(xí)來加快進化過程[9].高楠等人將遺傳算法與基因組排序相結(jié)合,以產(chǎn)生一種可以在有限的時間和空間內(nèi)解決DCJ中值問題的新方法[10].
層次聚類算法輸出最主要的部分是聚類樹,其可以為研究分層聚類的結(jié)果提供直觀的表現(xiàn)形式[11],由于層次聚類的邏輯原理易于理解,通用性較好,具有準(zhǔn)確的聚類結(jié)果并且可以根據(jù)一定閾值,自行選擇在所需級別對聚類樹進行“剪枝”操作以獲得合適的聚類等特性,因此廣泛應(yīng)用于大數(shù)據(jù)挖掘中,特別是在生物信息學(xué)中應(yīng)用更為普遍.
對于高通量測序得到的基因表達量數(shù)據(jù),一般具有高噪聲、高維度與非對稱性等特點.雖然聚類是智能算法,但傳統(tǒng)的聚類算法存在著易陷入局部最優(yōu)解、易受噪聲干擾和時空復(fù)雜度高等問題,在應(yīng)用于基因表達量數(shù)據(jù)時不能有效準(zhǔn)確地聚類.
基于以上問題,本文提出了一種基于改進遺傳算法的層次聚類算法(HCIGA:Hierarchical Clustering of Improving Genetic Algorithm),HCIGA融合了精英保留法與輪盤賭的選擇算子,優(yōu)化了適應(yīng)度函數(shù),采用了自適應(yīng)交叉率和變異率,并引入小生境技術(shù)保持種群的多樣性,提高算法的全局搜索能力,并在層次聚類中使用有序表來降低復(fù)雜度.最后分別用兩組異質(zhì)基因數(shù)據(jù)集驗證算法的有效性.
在層次聚類算法中,數(shù)據(jù)被逐層合并或者分解成一棵層次分明的聚類樹.按照樹的生成步驟可分為兩大類:自頂而下的分裂式層次聚類算法和自底而上的凝聚型層次聚類算法.針對分裂式層次聚類算法計算復(fù)雜度高的問題,本文采用凝聚型層次聚類算法為研究對象.
對于基因表達數(shù)據(jù)來說,層次聚類算法的優(yōu)點是通過設(shè)置不同的相關(guān)參數(shù)值,可以得到不同粒度水平上的多層次聚類結(jié)構(gòu);在聚類形狀方面,層次聚類適用于任意形狀的聚類;層次聚類算法可以方便地構(gòu)建出一棵基因聚類樹,充分反映基因之間相似性的關(guān)系.
雖然層次聚類算法在基因表達量數(shù)據(jù)中的應(yīng)用有諸多優(yōu)勢,但層次聚類算法本身還存在著一些不足之處:聚類操作需要使用相似度矩陣,導(dǎo)致算法的時間和空間復(fù)雜度大,所以必須進行優(yōu)化后才能在較大的數(shù)據(jù)集中使用;層次聚類的結(jié)果取決于上一個合并點或分裂點,下一次聚類會在前一次聚類基礎(chǔ)之上繼續(xù)進行合并或分裂,一旦聚類結(jié)果形成,已完成的聚類不能被撤消,簇間也不能交換對象,即層次聚類最大的缺陷就是不可逆性.因此選擇合適的分裂或聚合點是層次聚類的關(guān)鍵.
現(xiàn)有常用聚類算法在某些通用數(shù)據(jù)集中能取得較為理想的聚類結(jié)果,但當(dāng)應(yīng)用于異質(zhì)基因表達量數(shù)據(jù)集時,由于基因數(shù)據(jù)的非對稱性以及分離度不明顯,會導(dǎo)致聚類效果不佳.同時現(xiàn)有常用聚類算法僅基于統(tǒng)計學(xué)而很少考慮異質(zhì)基因內(nèi)部聯(lián)系與患者預(yù)后效果,聚類結(jié)果不一定具備生物學(xué)研究意義.因此本文在遺傳聚類算法中融合生物學(xué)分析,使HCIGA算法在異質(zhì)基因表達量數(shù)據(jù)中具有生物學(xué)意義.
傳統(tǒng)遺傳算法雖然在初期能保持種群的多樣性,但隨著進化的逐步深入,多個個體易集中于某一極值點上,陷入局部最優(yōu)解的陷阱.因此本文引入了小生境技術(shù)[12],融合改進的遺傳算法與凝聚型層次聚類,對異質(zhì)基因表達量數(shù)據(jù)集進行聚類分析,解決基因數(shù)據(jù)的高維度、高噪聲和非對稱性問題,并獲得較高的聚類精確度和較快的收斂速度.HCIGA算法的流程如圖1所示.
圖1 HCIGA流程圖Fig.1 Flow chart of HCIGA
HCIGA先對基因表達量數(shù)據(jù)預(yù)處理,對基質(zhì)細(xì)胞和免疫細(xì)胞計算腫瘤純度得分,并通過Wilcoxon秩和檢驗進行差異分析,對篩選的差異基因進行層次聚類,利用適應(yīng)度函數(shù)F(x)計算種群中所有染色體的適應(yīng)度值,在未達到停止條件時:將適應(yīng)度值降序排列,選擇適應(yīng)度高的優(yōu)良個體組成新種群,對父代進行選擇、自適應(yīng)交叉和變異,通過排擠小生境策略不斷更新種群.當(dāng)種群滿足終止條件時,將最佳染色體種群解碼并利用改進的凝聚型層次聚類算法對樣本進行聚類得到最佳亞型.
2.2.1 編碼方式
基因表達數(shù)據(jù)是一個n×m維的實值矩陣M[xij]n*m,如圖2所示,其中n代表基因的個數(shù),m表示實驗樣本總數(shù),行名一般用基因名稱標(biāo)注,列名為樣本名稱,矩陣中每個元素xij表示第i(1≤i≤n)個基因Gi在第j(1≤j≤m)個樣本Sj中的表達量.
圖2 基因表達矩陣示意圖Fig.2 Schematic diagram of gene expression matrix
對于遺傳算法編碼方式的選擇,綜合考慮染色體基因的選擇、交叉、變異等運算后,本文采用二進制編碼方式,“0”表示未選擇基因,“1”表示選擇基因,因此種群個體可以編碼為長度為n位的“0”和“1”字符串.采用這種編碼方式可以大大提高效率,在解碼時只需將“1”對應(yīng)的基因挑選出來即可.
2.2.2 種群初始化
雖然隨機初始化種群可以保證種群的多樣性,避免種群陷入局部最優(yōu)解,但在隨機初始化種群的過程中,容易產(chǎn)生非法編碼.因此,本文對種群的初始化加以限制:隨機生成正整數(shù)a控制選取的基因數(shù),10/n≤a≤n,其中n為差異表達基因數(shù),a表示選取基因數(shù),即隨機將a個差異基因位置“1”,其余n-a個差異基因位置“0”.并利用凝聚型層次聚類隨機生成k個初始中心.
2.2.3 適應(yīng)度函數(shù)
適應(yīng)度函數(shù)是控制種群遺傳的重要部分,本文根據(jù)基因表達量數(shù)據(jù)的特點,綜合考慮簇的規(guī)模和簇內(nèi)外相似度,優(yōu)化后的適應(yīng)度函數(shù)F(x)如公式(1)所示:
(1)
其中ni為第i個簇內(nèi)個體總數(shù),n為種群內(nèi)個體總數(shù),k為簇的個數(shù),Li為第i個簇內(nèi)相似度,S為簇間相似度,α為平衡系數(shù),用于平衡簇內(nèi)與簇間相似度的比例,不同的數(shù)據(jù)集對應(yīng)不同的α.
簇內(nèi)相似度Li定義如公式(2)所示:
(2)
其中m為樣本總數(shù),xij為基因Gi在樣本Sj中的表達量,Mi是第i簇基因的樣本中心:
(3)
簇間相似度S定義如公式(4)所示:
(4)
其中σ為高斯常數(shù),為了計算的方便,可以取2σ2=1,d2(Mi,Mj),為簇樣本中心Mi和Mj的距離.
在對樣本進行聚類的過程中,若簇內(nèi)相似度越大,且簇間相似度越小,則聚類的效果越理想.
2.2.4 遺傳算子
選擇算子:由于輪盤賭選擇是隨機操作,適應(yīng)度值高的染色體也有被淘汰的可能性,為了解決這個問題,使最優(yōu)個體被保留,本文設(shè)計了融合精英保留法與輪盤賭的選擇算子,對所有染色體的適應(yīng)度降序排列,前5%的優(yōu)秀染色體直接遺傳給下一代,后95%的染色體參與輪盤賭選擇,以此來保持種群中優(yōu)良個體.
交叉概率Pc和變異概率Pm是影響遺傳算法性能的關(guān)鍵,而Pc和Pm難以精確設(shè)定.本文采用了自適應(yīng)策略根據(jù)種群的適應(yīng)度來動態(tài)地調(diào)整Pc和Pm.
交叉算子:由于單點交叉和多點交叉的搜索空間較小,同時還可能會引起位置的偏差,本文采用效果更好的均勻交叉作為交叉算子,可以有效提高算法的局部搜索能力.
主要操作如下:
1)從種群中隨機選取一對交叉的父染色體.
2)隨機生成與染色體等長的“0”和“1”掩碼串,若掩碼當(dāng)前位為“0”,代表由父染色體2提供變量值,若掩碼相應(yīng)位為“1”,則代表由父染色體1提供變量值.
改進后的自適應(yīng)交叉概率Pci如公式(5)所示:
(5)
其中Pci為第i代的交叉概率;Pc0為初始交叉概率;k1為交叉比例系數(shù);Geni為當(dāng)前執(zhí)行第i代遺傳操作;Genmax為最大遺傳代數(shù);Fc為兩個待交叉染色體適應(yīng)度的平均值;Fa為種群適應(yīng)度的平均值.
變異算子:由于本文采用的是二進制編碼方式,在進行基因位變異運算時可以輕松實現(xiàn).具體過程為:從種群中隨機抽取一條染色體,將其以自適應(yīng)變異概率Pmi改變某位值,以此擴大算法搜索空間.
改進后的自適應(yīng)變異概率Pmi計算公式如式(6)所示:
(6)
其中Pmi為第i代的變異概率;Pm0為初始變異概率;k2為變異比例系數(shù);Geni為當(dāng)前執(zhí)行第i代遺傳操作;Genmax為最大遺傳代數(shù);Fm為待變異染色體適應(yīng)度值;Fa為種群適應(yīng)度的平均值.
2.2.5 小生境技術(shù)
為了提高種群的多樣性和全局最優(yōu)解搜索能力,本文在遺傳算法中引入基于排擠機制的小生境技術(shù).在小生境遺傳算法中采用海明距離Hij來表示染色體Ri和Rj之間的相似度,如公式(7)所示:
(7)
其中rik表示染色體Ri的第k位.很明顯,當(dāng)Ri和Rj之間的海明距離Hij越小,兩染色體之間的相似度越高,當(dāng)Hij小于某一閾值γ時(γ為某一很小的正整數(shù)),染色體Ri和Rj即屬于同一小生境,此時比較兩者的適應(yīng)度值,對適應(yīng)度值小的染色體設(shè)置懲罰函數(shù)F′(x)=0.6*F(x),以此降低其適應(yīng)度,減小其在子代中的遺傳概率.經(jīng)過排擠小生境操作后,適應(yīng)度在閾值γ內(nèi)的兩個染色體,僅需保留適應(yīng)度較高者,極大地提高了種群的多樣性.同時在種群個體數(shù)量恒定的情況下,使染色體在整個解空間內(nèi)分散開來,增強了全局搜索能力,從而解決局部最優(yōu)解的問題.
2.2.6 改進的層次聚類算法
由于傳統(tǒng)的層次聚類算法的復(fù)雜度達到了O(n3),在高維數(shù)據(jù)集中效率低下,根本原因是傳統(tǒng)層次聚類算法在每次合并兩個相似度最大的簇后需要重新生成相似度矩陣,這大大消耗了時間與資源.本文利用有序表對各基因相似度進行降序存儲,提高層次聚類算法的效率.
在相似性度量方面,通用的聚類算法采用歐幾里德距離作為樣本間相似性的度量,這在非基因表達量數(shù)據(jù)集是可行的,但由于多數(shù)情況下,基因表達量是實驗樣本與對照組之間的Ratio值,歐式距離無法準(zhǔn)確刻畫基因間的相似度.所以本文選用皮爾遜相關(guān)系數(shù)作為基因間相似性的度量.
設(shè)在n×m維的基因表達量矩陣M中,存在兩個基因組向量Ga=(xa1,xa2,…,xam)和Gb=(xb1,xb2,…,xbm),1≤a≤m,1≤b≤m,xij表示第i(1≤i≤n)個基因Gi在第j(1≤j≤m)個樣本Sj中的表達量.則基因Ga與Gb的相似度可以定義為Sab:
(8)
(9)
(10)
在相似性度量設(shè)計方面,本文選擇均鏈接法,綜合考慮了兩個類間各基因表達值之間的相關(guān)性,從而可以有效避免類連鎖現(xiàn)象.
本文實驗驗證的開發(fā)環(huán)境為PyCharm 2020.2.1,數(shù)據(jù)預(yù)處理軟件為RStudio 1.2.1335.
為了綜合評估本算法的性能,本文設(shè)計了兩類實驗,第1類為使用經(jīng)過實驗驗證并且?guī)?biāo)簽的4組基因表達量數(shù)據(jù)集:腦腫瘤[13]、肺癌[14]、腎臟癌[15]和乳腺癌[16],用于HCIGA精確度的評估,4組數(shù)據(jù)集的屬性描述如表1.第2類為使用從癌癥基因組圖譜官網(wǎng)(TCGA)中獲取不帶標(biāo)簽的HNSC數(shù)據(jù)集進行聚類分析.從TCGA數(shù)據(jù)庫中下載的HNSC數(shù)據(jù)集包括19645個基因,502個腫瘤樣本和44個正常樣本,以及528例臨床樣本.
表1 外部標(biāo)準(zhǔn)數(shù)據(jù)集參數(shù)Table 1 Parameters of external standard data sets
由于基因表達量數(shù)據(jù)的測量條件不同,各樣本之間會存在一些差異,且并不是所有的基因組都與疾病的產(chǎn)生存在影響,所以在實驗之前必須對數(shù)據(jù)進行清洗.
對于存在少量缺失值的樣本,本文采用KNN填補缺失表達量;若缺失值較多或表達量為負(fù)值的基因樣本則直接刪除;使用RMA和分位數(shù)歸一化分別建立這些芯片矩陣數(shù)據(jù)集的基因表達水平,以獲得每個基因和樣品的表達值[17],并對基因表達量數(shù)據(jù)進行標(biāo)準(zhǔn)化.利用ESTIMATE[18]對腫瘤樣本進行打分,并根據(jù)得分的中位值,將腫瘤樣本劃分為高得分組和低得分組.本文使用顯著性分析[19]檢測差異表達基因,利用Wilcoxon秩和檢驗,設(shè)置閾值|logFC|>1且FDR<0.05進行差異基因過濾.
遺傳算法部分重要實驗參數(shù)匯總見表2.
表2 主要實驗參數(shù)Table 2 Main experimental parameters
3.3.1 外部指標(biāo)性能度量
本文使用經(jīng)典層次聚類算法AGNES、遺傳K-Means算法GKA、DPeak算法、譜聚類算法和改進遺傳層次聚類算法HCIGA在4個帶類別標(biāo)簽的標(biāo)準(zhǔn)數(shù)據(jù)集(Brain Cancer、Lung Cancer、Renal Cancer、Breast Cancer)進行聚類分析,經(jīng)過10次重復(fù)獨立運行,得到五種算法在各數(shù)據(jù)集中平均聚類準(zhǔn)確率,如圖3所示.
圖3 5種算法在4個數(shù)據(jù)集中正確率的對比Fig.3 Accuracy of the five algorithms in the four data sets
由圖3可以看出AGNES在4個數(shù)據(jù)集中聚類精度都是最低,GKA和SC算法的聚類精度在Brain Cancer和Breast Cancer數(shù)據(jù)集中效果相當(dāng),但低于HCIGA和DPeak算法,而DPeak算法在Renal Cancer數(shù)據(jù)集中聚類精度不佳,低于GKA、譜聚類算法和HCIGA算法,僅優(yōu)于AGNES算法,HCIGA算法在Brain Cancer等4個數(shù)據(jù)集中的聚類精度均為最高,聚類效果最佳.
根據(jù)第2節(jié)HCIGA算法分析部分,計算腫瘤純度得分的時間可以忽略不計,相應(yīng)的時間復(fù)雜度主要包括遺傳算法和層次聚類兩個階段.設(shè)m為種群大小,n為染色體長度.在遺傳操作階段,初始化種群的時間可忽略不計,適應(yīng)度函數(shù)計算的時間復(fù)雜度為O(mn),在選擇操作中,主要采用的是輪盤賭,在最壞的情況下需要遍歷所有對象,因此時間復(fù)雜度為O(m2),在交叉和變異操作階段的時間復(fù)雜度均為O(m),為了提高算法的全局搜索能力,引入小生境技術(shù)增加了一定的時間復(fù)雜度O(m).在聚類操作階段,采用有序表結(jié)構(gòu)代替每次重復(fù)的計算相似度矩陣,可使聚類操作的時間復(fù)雜度降至O(n2),故算法HCIGA的時間復(fù)雜度為O(m2+n2),稍遜于DPeak等主流聚類算法的時間復(fù)雜度O(n2)和GKA的O(m2+nkt),但優(yōu)于譜聚類的O(n3),且聚類精度對比于另外4種算法有所提升.
圖4 5種算法在4個數(shù)據(jù)集中RAND指數(shù)折線圖Fig.4 Five algorithms in four data sets of RAND index line graph
3.3.2 內(nèi)部指標(biāo)性能度量
TCGA-HNSC數(shù)據(jù)集無模型給定的簇劃分,因此無法使用外部指標(biāo)來衡量算法模型的性能,本文選用了收斂速度與生存分析作為HCIGA有效性的內(nèi)部評價指標(biāo).
圖5 HNSC數(shù)據(jù)集的聚類樹Fig.5 Clustering tree of HNSC data set
對聚類個數(shù)的確定是影響聚類結(jié)果的關(guān)鍵因素.本文將HCIGA在HNSC數(shù)據(jù)集中得到的聚類樹進行剪枝操作[20],如圖5所示,選取聚類數(shù)K為2、3、4、6、11、13這6種方案,并通過Silhouette指數(shù)對這六種方案的聚類數(shù)評估,Silhouette指數(shù)見表3.
表3 各聚類數(shù)的Silhouette指數(shù)值Table 3 Silhouette indices of each clustering number
Silhouette指數(shù)值能很好地反應(yīng)聚類類別的優(yōu)劣,最大值對應(yīng)著最優(yōu)劃分,由表3可以看出,當(dāng)聚類數(shù)K=6時,所對應(yīng)的Silhouette指數(shù)值最大,因此最終確定將聚類結(jié)果劃分為6簇的方案.
在收斂速度評價指標(biāo)方面,HCIGA與GKA的收斂速度曲線如圖6所示,HCIGA在運行到第30代左右時,適應(yīng)度值已基本達到最優(yōu);而GKA在第50代附近才取得最優(yōu)適應(yīng)度值,由此可見,HCIGA的收斂速度明顯優(yōu)于GKA算法.
圖6 HCIGA與GKA收斂速度Fig.6 Convergence rates of HCIGA and GKA
針對TCGA-HNSC的生存分析,如圖7所示,可以分別繪制出HNSC 6種亞型的總體生存曲線,圖示6條生存曲線走勢各不相同,隨著生存時間的推移,病人的存活率逐漸下降,亞型間的p值為6.98e-03小于0.05,這表明6種亞型的生存曲線具有統(tǒng)計學(xué)價值,即本算法聚類后的亞型結(jié)果與生存密切相關(guān),可以很好地區(qū)分HNSC亞型.
圖7 各亞型間生存曲線圖Fig.7 Survival curves among subtypes
隨著高通量測序技術(shù)的蓬勃發(fā)展,基因數(shù)據(jù)井噴式增長,為了有效處理基因表達量數(shù)據(jù),聚類分析作為一種無監(jiān)督數(shù)據(jù)分析方法,成為基因研究的重要手段之一.本文在凝聚型層次聚類算法的基礎(chǔ)上結(jié)合ESTIMATE評分體系、差異分析以及改進遺傳算法,提出面向異質(zhì)基因數(shù)據(jù)的智能層次聚類算法HCIGA.
本文將HCIGA在4個帶類別標(biāo)簽數(shù)據(jù)集上的聚類精度與AGNES、GKA、DPeak以及SC算法進行對比,結(jié)果表明HCIGA的聚類精度在4個數(shù)據(jù)集中均為最高.本文并將HCIGA應(yīng)用于TCGA數(shù)據(jù)庫中無類別標(biāo)簽數(shù)據(jù)集HNSC,通過收斂速度以及生存分析等方法證明,HCIGA較GKA在收斂速度上取得明顯提升,且利用HCIGA在HNSC中識別的6種亞型具有統(tǒng)計學(xué)意義和生物學(xué)價值.