基于WGCNA聯(lián)合LASSO篩選宮頸癌預(yù)后樞紐基因

郭依琳 王 璐 徐 臻 趙 虎 韓少聰 王武亮

宮頸癌是嚴(yán)重威脅女性健康的惡性腫瘤,在全球范圍內(nèi)發(fā)生率和病死率均居女性惡性腫瘤第4位[1]。在我國(guó)宮頸癌的發(fā)生率和病死率分別位居第6位和第8位,盡管近年發(fā)生率有所下降,但發(fā)病年齡逐漸呈年輕化[2, 3]。HPV是宮頸癌發(fā)病的主要原因[4]。早期宮頸癌患者以手術(shù)治療為主,預(yù)后較好。中晚期及復(fù)發(fā)的宮頸癌患者以放化療治療為主,預(yù)后較差[5]。目前可以用來(lái)預(yù)測(cè)宮頸癌患者預(yù)后的分子標(biāo)志物較少。因此,尋找合適的預(yù)后生物學(xué)標(biāo)志物和治療靶點(diǎn)可能有助有提高宮頸癌患者的總生存率。

近年來(lái),隨著基因組學(xué)和生物信息學(xué)的不斷發(fā)展,利用高通量測(cè)序和基因芯片技術(shù)篩選宮頸癌預(yù)后相關(guān)標(biāo)志物具有重要研究?jī)r(jià)值。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis,WGCNA)是一種識(shí)別與表型性狀相關(guān)的基因模塊和關(guān)鍵樞紐基因的方法[6]。最小絕對(duì)值選擇與收縮算子(least absolute selection and shrinkage operator,LASSO)是一種通過(guò)構(gòu)造懲戒函數(shù)壓縮部分低權(quán)重回歸系數(shù),明確兩變量之間的關(guān)聯(lián)程度的回歸算法[7]。WGCNA聯(lián)合LASSO算法可以更好的篩選出相關(guān)性更高的樞紐基因。因此,本研究基于基因表達(dá)匯編(Gene Expression Omnibus,GEO)數(shù)據(jù)庫(kù)、腫瘤基因組圖譜(The Cancer Genome Atlas,TCGA)數(shù)據(jù)庫(kù)和基因型和基因表達(dá)量關(guān)聯(lián)(Genotype Tissue Expression,GTEx)數(shù)據(jù)庫(kù),運(yùn)用WGCNA聯(lián)合LASSO算法篩選出宮頸癌預(yù)后相關(guān)的分子標(biāo)志物,并初步分析其與宮頸癌腫瘤微環(huán)境(tumor microenvironment,TME)中免疫細(xì)胞浸潤(rùn)的相關(guān)性,為宮頸癌發(fā)生、發(fā)展中的分子機(jī)制研究提供依據(jù)。

資料與方法

1.數(shù)據(jù)的獲取和處理:從GEO數(shù)據(jù)庫(kù)選取GSE9750和GSE52903數(shù)據(jù)集。GSE9750數(shù)據(jù)集包括33例宮頸癌組織和24例正常宮頸組織;GSE52903數(shù)據(jù)集包括55例宮頸癌組織和17例正常宮頸組織。從TCGA數(shù)據(jù)庫(kù)選取304例子宮頸癌組織和3例正常宮頸組織。從GTEx數(shù)據(jù)庫(kù)選取10例正常宮頸組織。下載芯片及RNA測(cè)序轉(zhuǎn)錄組數(shù)據(jù)和對(duì)應(yīng)的臨床特征數(shù)據(jù)。在R軟件(版本號(hào)4.1.2)中利用“sva”包中的“combat()”函數(shù)對(duì)GSE9750和GSE52903數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化處理,合并和校正,主成分分析(Principal component analysis,PCA)檢驗(yàn)校正的數(shù)據(jù)。利用“l(fā)imma”包中的“normalizeBetweenArrays()”函數(shù)對(duì)TCGA和GTEx數(shù)據(jù)集進(jìn)行合并和校正。

2.差異表達(dá)基因的篩選和功能富集分析:在R軟件中使用“l(fā)imma”包篩選宮頸癌中差異表達(dá)的基因(differentially expressed gene,DEG),篩選標(biāo)準(zhǔn)P<0.05且|log2差異倍數(shù)(fold change,FC)|>2。繪制DEG的火山圖。在R軟件中使用“cluserProfiler”包對(duì)宮頸癌和正常宮頸中DEG進(jìn)行基因本體論(Gene Ontology,GO)功能學(xué)。錯(cuò)誤發(fā)現(xiàn)率(false discovery rate,FDR)值<0.05被認(rèn)為差異有統(tǒng)計(jì)學(xué)意義。

3.WGCNA模塊構(gòu)建及可視化:在R軟件中使用“WGCNA”包,刪除GEO數(shù)據(jù)集(GSE9750和GSE52903)中基因表達(dá)離散程度較小(sd≤0.7)的基因,使用“goodsamplegenes()”函數(shù)檢查缺失值。使用“pickSoftThreshold()”函數(shù)計(jì)算共表達(dá)網(wǎng)絡(luò)的軟閾值。當(dāng)軟閾值等于4,共表達(dá)網(wǎng)絡(luò)更逼近無(wú)尺度網(wǎng)絡(luò)。構(gòu)建加權(quán)鄰接矩陣,基于拓?fù)渲丿B矩陣(TOM)的相異度(1-TOM)進(jìn)行層次聚類,構(gòu)建相關(guān)模塊。將模塊與臨床特征數(shù)據(jù)關(guān)聯(lián)后,繪制模塊身份(module membership,MM)及基因顯著性(gene significance,GS)的散點(diǎn)圖,明確模塊內(nèi)基因的顯著性。

4.樞紐基因的篩選:WGCNA中按照MM >0.8且GS >0.5篩選模塊內(nèi)的基因。將篩選到的模塊基因與差異表達(dá)基因取交集,并在R軟件中使用“glmnet”包,在α=1的條件下,選擇合適的λ值,用十折交叉驗(yàn)證法實(shí)現(xiàn)內(nèi)部驗(yàn)證,確定最佳樞紐基因。使用TCGA可視化分析工具GEPIA,進(jìn)一步篩選出與預(yù)后相關(guān)的樞紐基因。

5.樞紐基因在宮頸癌和正常宮頸組織的表達(dá):在GEO數(shù)據(jù)集(GSE9750和GSE52903)和TCGA聯(lián)合GTEx數(shù)據(jù)集中,檢測(cè)預(yù)后相關(guān)的樞紐基因在宮頸癌組織與正常宮頸組織之間的表達(dá)情況。人類蛋白圖譜數(shù)據(jù)庫(kù)(human protein atlas,HPA)分析與預(yù)后相關(guān)的樞紐基因在宮頸癌組織與正常宮頸組織之間的蛋白表達(dá)情況。

6.免疫浸潤(rùn)和免疫檢查點(diǎn)基因相關(guān)性分析:利用聚類分析和單樣本基因集富集分析(single sample gene set enrichment analysis,ssGSEA)計(jì)算宮頸癌組織和正常宮頸組織中免疫細(xì)胞浸潤(rùn)程度,繪制小提琴圖。計(jì)算預(yù)后相關(guān)的樞紐基因與28種免疫細(xì)胞以及36個(gè)免疫檢查點(diǎn)基因的相關(guān)性,以P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

結(jié) 果

1.納入患者臨床特征分析:本研究共納入392例宮頸癌組織和54例正常宮頸組織的轉(zhuǎn)錄組表達(dá)數(shù)據(jù)。GSE9750、GSE52903和TCGA數(shù)據(jù)集具有完整的臨床病理數(shù)據(jù)。3組數(shù)據(jù)集在年齡分布和病理類型之間比較,差異無(wú)統(tǒng)計(jì)學(xué)意義(P<0.05),但在臨床分期和生存狀態(tài)之間比較,差異有統(tǒng)計(jì)學(xué)意義(P<0.05,表1)。

表1 392例宮頸癌組織的臨床病理特征[n(%)]

2.宮頸癌中差異表達(dá)的基因和功能富集分析:將GSE9750和GSE52903數(shù)據(jù)集合并且校正后,利用PCA主成分分析檢驗(yàn)校正后的數(shù)據(jù),未發(fā)現(xiàn)存在批次效應(yīng)(圖1)。在88例宮頸癌樣本和41例正常宮頸樣本中,根據(jù)設(shè)置的篩選條件,共發(fā)現(xiàn)109個(gè)DEG,其中表達(dá)上調(diào)38個(gè),表達(dá)下調(diào)71個(gè),并繪制火山圖(圖2)。根據(jù)P值選擇表達(dá)差異最顯著的前10個(gè)基因(表2)。對(duì)109個(gè)DEG進(jìn)行GO功能注釋分析。

圖1 校正后的GSE9750和GSE52903數(shù)據(jù)集PCA分析

圖2 宮頸癌中差異表達(dá)基因的火山圖

表2 宮頸癌中差異表達(dá)最顯著前10個(gè)基因

GO功能結(jié)果注釋發(fā)現(xiàn),差異表達(dá)基因可能參與109條通路的功能,其中生物過(guò)程(biological process,BP)富集度最大的3條通路為皮膚發(fā)育、表皮發(fā)育和表皮細(xì)胞分化;細(xì)胞組分(cellular component,CC)富集度最大的3條通路為角化包膜、中心顆粒體和DNA復(fù)制前起始復(fù)合物;分子功能(molecular function,MF)富集度最大通路為絲氨酸酶活性(圖3)。

圖3 宮頸癌中差異基因的GO功能富集分析

3.WGCNA網(wǎng)絡(luò)構(gòu)建及顯著性模塊確認(rèn):對(duì)GSE9750和GSE52903數(shù)據(jù)集中所有基因構(gòu)建共表達(dá)網(wǎng)絡(luò)。當(dāng)軟閾值等于4時(shí),共表達(dá)網(wǎng)絡(luò)接近為無(wú)尺度網(wǎng)絡(luò)。確定軟閾值后,采用動(dòng)態(tài)剪切法,獲得9個(gè)模塊(圖4A)。進(jìn)一步分析各模塊的特征向量,發(fā)現(xiàn)turquoise模塊(包含126個(gè)基因)與宮頸癌的相關(guān)性最高(r=0.87,P=2×10-40),并且GS和MS最高,確認(rèn)為顯著性塊(圖4B)。

圖4 WGCNA網(wǎng)絡(luò)構(gòu)建及顯著性模塊確認(rèn)A.基因共表達(dá)網(wǎng)絡(luò)和共表達(dá)模塊的層次聚類樹(shù);B.特征模塊與臨床特征之間的關(guān)系。紅色為正相關(guān),藍(lán)色為負(fù)相關(guān)

4.樞紐基因的篩選:將turquoise模塊中126個(gè)基因與差異表達(dá)109個(gè)基因取交集,得到27個(gè)基因。利用LASSO模型進(jìn)行內(nèi)部驗(yàn)證,最終確定了7個(gè)樞紐基因(圖5A)。GEPIA數(shù)據(jù)庫(kù)進(jìn)一步分析得到3個(gè)預(yù)后相關(guān)的樞紐基因,分別為MCM2、APOD和RAD54L(圖5中B～D)。

圖5 宮頸癌中篩選樞紐基因A.LASSO模型建立;B.MCM2表達(dá)水平的生存曲線圖;C.APOD表達(dá)水平的生存曲線圖;D.RAD54L表達(dá)水平的生存曲線圖

圖6 宮頸癌中預(yù)后相關(guān)樞紐基因的表達(dá)情況A.GSE9750和GSE52903數(shù)據(jù)集;B.TCGA和GTEx數(shù)據(jù)集。*P<0.001

5.預(yù)后相關(guān)樞紐基因在宮頸癌和正常宮頸組織表達(dá)情況:在GSE9750和GSE52903數(shù)據(jù)集中,與正常宮頸組織比較,MCM2和RAD54L在宮頸癌組織中高表達(dá),而APOD則低表達(dá)(圖6A)。在TCGA聯(lián)合GTEx數(shù)據(jù)集中,3個(gè)樞紐基因的表達(dá)結(jié)果與GEO數(shù)據(jù)集(GSE9750和GSE52903)結(jié)果一致(圖6B)。為了在蛋白水平驗(yàn)證3個(gè)預(yù)后相關(guān)的樞紐基因的表達(dá)情況,筆者在HPA數(shù)據(jù)庫(kù)中分析了其在宮頸癌組織和正常宮頸組織中免疫組化的表達(dá)情況。結(jié)果發(fā)現(xiàn)MCM2和RAD54L在宮頸癌組織中的表達(dá)高于正常組織,而APOD在宮頸癌組織和正常宮頸組織都屬于低表達(dá)(圖7)。

6.預(yù)后相關(guān)的樞紐基因與免疫相關(guān)性分析:ssGSEA方法計(jì)算宮頸癌組織和正常宮頸組織中免疫細(xì)胞細(xì)胞浸潤(rùn)程度,初步探究宮頸癌TME中免疫細(xì)胞浸潤(rùn)程度。結(jié)果發(fā)現(xiàn)激活的B淋巴細(xì)胞、激活的CD4+T細(xì)胞、激活的樹(shù)突狀細(xì)胞、調(diào)節(jié)性T細(xì)胞和輔助性T細(xì)胞在宮頸癌組織中表達(dá)上調(diào);而肥大細(xì)胞、嗜酸性粒細(xì)胞和中性粒細(xì)胞中則在宮頸癌組織中表達(dá)下調(diào)(圖8A)。進(jìn)一步探究3個(gè)預(yù)后相關(guān)的樞紐基因與28種免疫細(xì)胞以及36個(gè)免疫檢查點(diǎn)基因的相關(guān)性。結(jié)果發(fā)現(xiàn),在28種免疫細(xì)胞中,MCM2基因與CD56bright.NK細(xì)胞呈明顯正相關(guān),與肥大細(xì)胞呈明顯負(fù)相關(guān);APOD基因與漿細(xì)胞樣樹(shù)突細(xì)胞和激活的B細(xì)胞呈明顯正相關(guān);RAD54L基因與大部分免疫細(xì)胞呈負(fù)相關(guān)(圖8B)。在36個(gè)免疫檢查點(diǎn)基因中,MCM2基因與免疫檢查點(diǎn)LGALS9呈明顯正相關(guān),與NRP1呈明顯負(fù)相關(guān);APOD基因與CD70和TNFSF9呈明顯正相關(guān);RAD54L基因與CD44呈明顯正相關(guān)(圖8C)。

討 論

宮頸癌傳統(tǒng)的治療方式主要為手術(shù)切除、放療和化療,不同患者對(duì)放化療的敏感度差異很大,并且即使患者接受上述治療后仍會(huì)出現(xiàn)復(fù)發(fā)和轉(zhuǎn)移的可能[8, 9]。宮頸癌的復(fù)發(fā)和轉(zhuǎn)移是目前造成宮頸癌患者死亡的主要原因,臨床上迫切需要發(fā)現(xiàn)更多分子標(biāo)志物用于宮頸癌的預(yù)后預(yù)測(cè)。因此,本研究基于GEO數(shù)據(jù)庫(kù)中GSE9750和GSE52903數(shù)據(jù)集,利用WGCNA和LASSO篩選宮頸癌樞紐基因,GEPIA數(shù)據(jù)庫(kù)進(jìn)一步篩選預(yù)后相關(guān)的樞紐基因;然后在GEO數(shù)據(jù)集(GSE9750和GSE52903)和TCGA聯(lián)合GTEx數(shù)據(jù)集中比較預(yù)后相關(guān)的樞紐基因在宮頸癌和正常宮頸組織中的表達(dá)情況,并在HPA數(shù)據(jù)庫(kù)中驗(yàn)證;最后利用ssGSEA分析宮頸癌TME免疫細(xì)胞浸潤(rùn)情況,探究3個(gè)預(yù)后相關(guān)的樞紐基因與免疫細(xì)胞浸潤(rùn)和免疫檢查點(diǎn)基因表達(dá)的相關(guān)性。

近年來(lái),隨著計(jì)算機(jī)技術(shù)及人工智能等在生物醫(yī)學(xué)方面的快速發(fā)展,生物信息學(xué)已經(jīng)成為基于大數(shù)據(jù)進(jìn)行分子標(biāo)志物篩選的有力手段之一。方萌等[10]研究利用WGCNA和機(jī)器深度算法確定了5個(gè)miR-141-3p的靶基因FOXA1、DMBX1、TMEM98、RHPN1、SRMS,并構(gòu)建了預(yù)后模型。Liu等[11]研究利用WGCNA方法,在GSE26511數(shù)據(jù)集中,發(fā)現(xiàn)ACKR1基因可能與宮頸癌淋巴結(jié)轉(zhuǎn)移和預(yù)后有關(guān)。同樣,Wang等[12]也利用相同的方法發(fā)現(xiàn)RRM2對(duì)宮頸癌的診斷和預(yù)后具有一定的價(jià)值。然而,運(yùn)用單一的生物信息學(xué)分析方法,可能造成干擾數(shù)據(jù)過(guò)多,對(duì)結(jié)果的準(zhǔn)確性造成影響。因此本研究選了WGCNA和LASSO模型兩種方法聯(lián)合分析增加研究結(jié)果的準(zhǔn)確性和真實(shí)性。

本研究在GSE9750和GSE52903數(shù)據(jù)集中共獲得109個(gè)差異表達(dá)的基因和126個(gè)WGCNA篩選turquoise模塊基因,取交集后得到27個(gè)樞紐基因,利用LASSO模型確定7個(gè)樞紐基因。利用GEPIA數(shù)據(jù)庫(kù)進(jìn)一步分析最終得到3個(gè)預(yù)后相關(guān)的樞紐基因,分別為MCM2、APOD和RAD54L。MCM2基因編碼微小染色體維持蛋白2,是DNA復(fù)制起始的主要調(diào)控因子之一。有研究發(fā)現(xiàn),MCM2的過(guò)表達(dá)與多種惡性腫瘤的發(fā)生、發(fā)展相關(guān)[13]。Das等[14]研究發(fā)現(xiàn),MCM2在宮頸癌組織和宮頸癌細(xì)胞系中高表達(dá),然而MCM2的表達(dá)水平與臨床病理分期無(wú)關(guān),這在Amaro等[15]研究中也被證實(shí)。APOD基因編碼載脂蛋白D,有研究發(fā)現(xiàn),APOD參與多種免疫反應(yīng)、細(xì)胞凋亡和腫瘤的發(fā)生、發(fā)展[16]。然而,APOD在宮頸癌方面的相關(guān)研究較少,僅Wang等[17]通過(guò)生物信息學(xué)技術(shù)確定了包括APOD基因在內(nèi)的10個(gè)免疫相關(guān)性基因。RAD54L基因在同源重組修復(fù)中發(fā)揮重要作用。有研究發(fā)現(xiàn)RAD54L在脈絡(luò)叢癌的發(fā)生、發(fā)展中起到重要作用,但其在宮頸癌中的研究較少。因此,本研究在GEO數(shù)據(jù)集(GSE9750和GSE52903)和TCGA聯(lián)合GTEx數(shù)據(jù)集中比較了3個(gè)預(yù)后相關(guān)的樞紐基因在宮頸癌組織和正常宮頸組織中的表達(dá)情況,發(fā)現(xiàn)與正常宮頸組織比較,MCM2和RAD54L在宮頸癌組織中高表達(dá),而APOD則低表達(dá),為后續(xù)3個(gè)樞紐基因在宮頸癌中的作用機(jī)制提供前期研究的依據(jù)。

有研究發(fā)現(xiàn),腫瘤的進(jìn)展不僅受到腫瘤本身生物學(xué)行為特征的影響,同時(shí)也有受到TME的影響。TME是腫瘤生長(zhǎng)、轉(zhuǎn)移和調(diào)節(jié)腫瘤免疫反應(yīng)的關(guān)鍵因素[18]。本研究采用ssGSEA方法,計(jì)算了宮頸癌組織和正常宮頸組織中28種免疫細(xì)胞浸潤(rùn)程度,發(fā)現(xiàn)大部分免疫細(xì)胞,如激活的B淋巴細(xì)胞、激活的CD4T細(xì)胞、激活的樹(shù)突狀細(xì)胞等在宮頸癌組中表達(dá)高于正常宮頸組織。另外,筆者還比較了3個(gè)預(yù)后相關(guān)的樞紐基因與免疫細(xì)胞和免疫檢查點(diǎn)基因表達(dá)的相關(guān)性,發(fā)現(xiàn)APOD基因與大多數(shù)免疫細(xì)胞的表達(dá)呈明顯負(fù)相關(guān)。特別地,3個(gè)預(yù)后相關(guān)的樞紐基因與部分免疫檢查點(diǎn)基因的表達(dá)呈正相關(guān)。這些結(jié)果提示篩選的3個(gè)預(yù)后分子標(biāo)志物與宮頸癌TME的免疫浸潤(rùn)水平存在關(guān)聯(lián),可能對(duì)宮頸癌的免疫治療可以提供一定的參考依據(jù)。

綜上所述,本研究利用WGCNA聯(lián)合LASSO方法篩選出了3個(gè)與預(yù)后相關(guān)的樞紐基因,MCM2、APOD和RAD54L,探究了其在宮頸癌和正常宮頸組織中的表達(dá),并初步分析了其與宮頸癌TME的免疫浸潤(rùn)水平的關(guān)系,為宮頸癌的預(yù)后預(yù)測(cè)和免疫治療提供指導(dǎo)意義。