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

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

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

資料與方法

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

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

3.WGCNA模塊構建及可視化:在R軟件中使用“WGCNA”包,刪除GEO數據集(GSE9750和GSE52903)中基因表達離散程度較小(sd≤0.7)的基因,使用“goodsamplegenes()”函數檢查缺失值。使用“pickSoftThreshold()”函數計算共表達網絡的軟閾值。當軟閾值等于4,共表達網絡更逼近無尺度網絡。構建加權鄰接矩陣,基于拓撲重疊矩陣(TOM)的相異度(1-TOM)進行層次聚類,構建相關模塊。將模塊與臨床特征數據關聯后,繪制模塊身份(module membership,MM)及基因顯著性(gene significance,GS)的散點圖,明確模塊內基因的顯著性。

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

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

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

結 果

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

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

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

圖1 校正后的GSE9750和GSE52903數據集PCA分析

圖2 宮頸癌中差異表達基因的火山圖

表2 宮頸癌中差異表達最顯著前10個基因

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

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

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

圖4 WGCNA網絡構建及顯著性模塊確認A.基因共表達網絡和共表達模塊的層次聚類樹;B.特征模塊與臨床特征之間的關系。紅色為正相關,藍色為負相關

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

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

圖6 宮頸癌中預后相關樞紐基因的表達情況A.GSE9750和GSE52903數據集;B.TCGA和GTEx數據集。*P<0.001

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

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

討 論

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

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

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

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

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