基于miRNA的乳頭狀腎細(xì)胞癌預(yù)后指標(biāo)的生物信息學(xué)研究

王松松,黃慧珍,林 凡,b*,蘇陳穎,李曉冉,林 夢,林曉暉,b

(福建中醫(yī)藥大學(xué)a.中西醫(yī)結(jié)合學(xué)院;b.中西醫(yī)結(jié)合慢性病研究福建省高校重點(diǎn)實(shí)驗(yàn)室,中國福建 福州 350122)

腎細(xì)胞癌(renal cell carcinoma,RCC)是一種異質(zhì)性疾病,具有不同的組織學(xué)和細(xì)胞遺傳學(xué)異常,占成人腎臟惡性腫瘤的90%,其發(fā)病率每年以2.5%的速度上升[1]。其中,乳頭狀腎細(xì)胞癌(papillary renal cell carcinoma,PRCC)是第二常見的腎臟惡性實(shí)質(zhì)腫瘤,占腎臟腫瘤的7%～14%[2～3]?，F(xiàn)有PRCC治療方法主要為手術(shù)治療。不同手術(shù)方案各有優(yōu)缺點(diǎn),如:根治性腎切除術(shù)可有效切除腫瘤,然而可能導(dǎo)致腎功能下降,增加慢性腎功能不全和透析的風(fēng)險(xiǎn),對患者生活質(zhì)量、心血管功能和總生存率都有不利影響;消融治療具有創(chuàng)傷小、療效確切的優(yōu)點(diǎn),但復(fù)發(fā)風(fēng)險(xiǎn)高于腎切除術(shù),遠(yuǎn)期療效需進(jìn)一步觀察[4]。因此,PRCC的治療方法仍需不斷優(yōu)化,而研究其預(yù)后標(biāo)志物則有助于優(yōu)化治療方法以及提高患者預(yù)后判斷的精確性[5]。目前,已經(jīng)有研究報(bào)道PRCC的預(yù)后標(biāo)志物,如Huang等[1]發(fā)現(xiàn) PVT1、PR11-401P9.4 等 6 個(gè)長鏈非編碼RNA(long noncoding RNA,lncRNA)可作為PRCC相關(guān)性預(yù)后生物標(biāo)志物。但此類研究報(bào)道甚少,且研究的深度和廣度都亟待加強(qiáng)。

微RNA(microRNA,miRNA)可降解靶基因mRNA,抑制靶蛋白質(zhì)翻譯,是一類基因表達(dá)調(diào)控的關(guān)鍵分子。其表達(dá)失調(diào)是導(dǎo)致疾病發(fā)生發(fā)展的重要因素之一[6]。同時(shí),由于miRNA相對分質(zhì)子量小,擁有較為穩(wěn)定的二級結(jié)構(gòu),檢測技術(shù)靈敏且成熟,所以其也是一類極具潛力的疾病預(yù)后生物標(biāo)志物。本研究以miRNA為分析對象,探究PRCC的預(yù)后生物指標(biāo),以期為臨床患者治療方案的選擇及預(yù)后的判斷提供更多依據(jù)。

1 數(shù)據(jù)與方法

1.1 miRNA與基因芯片表達(dá)數(shù)據(jù)下載

參照文獻(xiàn)[7]所述方法,利用編程語言R(version 4.0.0;https://www.r-project.org/)的RTCGA包[8]下載癌癥基因組圖譜數(shù)據(jù)庫(The Cancer Genome Atlas,TCGA;https://www.cancer.gov/)中PRCC相關(guān)的miRNA-seq數(shù)據(jù)KIRP.miRNASeq和臨床數(shù)據(jù)KIRP.clinical。該數(shù)據(jù)包含35個(gè)癌旁非癌腎組織和307個(gè)PRCC樣本,以及樣本的編號(ID)、種族(race)、年齡(age)、性別(gender)、病理階段(pathologic stage)、死亡天數(shù)/隨訪天數(shù)(days_to_death/days_to_last_followup)、死亡/存活(dead/alive)等臨床資料。

在GEO(Gene Expression Omnibus;https://www.ncbi.nlm.nih.gov/geo/)[9]中選取腫瘤實(shí)驗(yàn)組與非腫瘤對照組的組織樣本數(shù)均大于或等于3的芯片GSE7023[10]。此芯片包含12個(gè)癌旁非癌腎組織樣本和35個(gè)PRCC樣本。隨后,下載GSE7023中的PRCC基因芯片數(shù)據(jù)和GPL4866基因探針注釋。

1.2 miRNA表達(dá)的差異分析

對于KIRP.miRNASeq數(shù)據(jù),首先,以miRNA在超過10個(gè)樣本中的表達(dá)值大于1為標(biāo)準(zhǔn)[11],去除含缺失值和表達(dá)極其微量的miRNA數(shù)據(jù);其次,將數(shù)據(jù)進(jìn)行l(wèi)og2轉(zhuǎn)化處理;然后,應(yīng)用編程語言R的edgeR包[12]進(jìn)行差異分析,篩選出差異表達(dá)miRNA,篩選標(biāo)準(zhǔn):|log2(fold change)|≥1 且 P＜0.05;最后,利用R包“ggplot2”繪制差異表達(dá)miRNA的火山圖[13]。

1.3 生存分析

為分析上述差異表達(dá)miRNA與患者生存之間的相關(guān)性,采用survival包、survminer包進(jìn)行批量Kaplan-Meier生存分析,找出與PRCC患者生存顯著相關(guān)的miRNA,篩選標(biāo)準(zhǔn)為P＜0.01。

1.4 LASSO及Cox回歸分析

為了篩選出與PRCC患者預(yù)后強(qiáng)相關(guān)的miRNA,通過survival包對與生存顯著相關(guān)的miRNA進(jìn)行單因素Cox回歸分析,篩選標(biāo)準(zhǔn)為P＜0.05。為了避免過度擬合,運(yùn)用glmnet包[14]對單因素Cox回歸分析的結(jié)果進(jìn)行LASSO回歸分析。具體而言:通過隨機(jī)調(diào)整參數(shù)λ,擬合出50個(gè)系數(shù)(coefficient)不同的模型,再進(jìn)行交叉驗(yàn)證擬合,計(jì)算出lambda.min值;選擇λ為lambda.min時(shí)所對應(yīng)的模型,將該模型下系數(shù)不為0的miRNA篩選出來[15～16]。為了排除性別、年齡、腫瘤的發(fā)展階段等對預(yù)后的影響,把LASSO回歸分析篩選得到的miRNA進(jìn)行多因素Cox回歸分析,篩選條件為P＜0.05,以評估m(xù)iRNA作為獨(dú)立預(yù)后因素在患者生存中的作用,得到獨(dú)立影響預(yù)后的miRNA[17],并構(gòu)建生存曲線。

1.5 獨(dú)立影響預(yù)后的miRNA靶基因預(yù)測及表達(dá)分析

miRWalk3.0(http://mirwalk.umm.uni-heidelberg.de/)是分析、預(yù)測miRNA靶基因的重要生物信息學(xué)工具。它采用tarpmir算法對miRNA靶點(diǎn)進(jìn)行預(yù)測,通過隨機(jī)森林法計(jì)算出score。所得score顯示了miRNA與其靶點(diǎn)相互作用“有效”的概率,即score越接近1,預(yù)測的準(zhǔn)確性越大。本研究以miRWalk3.0網(wǎng)站推薦值score≥0.95為篩選條件[18～19],從miRWalk3.0的miRTarBase和 miRDB數(shù)據(jù)庫中獲取獨(dú)立影響預(yù)后的miRNA的靶基因。

利用PRCC基因芯片數(shù)據(jù)分析所得靶基因的差異表達(dá)。具體而言:首先,應(yīng)用編程語言R的dplyr包對GSE7023中的PRCC基因芯片數(shù)據(jù)進(jìn)行探針注釋轉(zhuǎn)換與重復(fù)值去除;然后,利用limma包[20]對其進(jìn)行差異表達(dá)分析,得到差異表達(dá)基因[|log2(fold change)|≥1 且 P＜0.05];最后,取差異表達(dá)基因與miRNA靶基因的交集,獲得差異表達(dá)的miRNA靶基因。

2 結(jié)果

2.1 miRNA表達(dá)的差異分析

所得PRCC miRNA-seq數(shù)據(jù)的差異表達(dá)分析結(jié)果表明,與癌旁非癌腎組織樣本相比,PRCC組織樣本中302個(gè)miRNA表達(dá)無差異,231個(gè)miRNA存在顯著差異表達(dá)[|log2(fold change)|≥1且P＜0.05]。在顯著差異表達(dá)的miRNA中,117個(gè)上調(diào)(占 50.6%),114 個(gè)下調(diào)(占 49.4%)(圖1)。

圖1 PRCC組織中miRNA表達(dá)情況的火山圖相對于癌旁非癌腎組織,黑色節(jié)點(diǎn)表示表達(dá)無差異的miRNA(356個(gè)),紅色節(jié)點(diǎn)表示表達(dá)上調(diào)的miRNA(117個(gè)),綠色節(jié)點(diǎn)表示表達(dá)下調(diào)的miRNA(114個(gè))。Fig.1 Volcano map of miRNA expression in PRCC tissuesCompared with normal tissues,black nodes represent 356 miRNAs without differential expression,red nodes represent 117 miRNAs with up-regulated expression,and green nodes represent 114 miRNAs with down-regulated expression.

2.2 生存分析

利用R語言對上述差異表達(dá)的miRNA進(jìn)行生存分析,發(fā)現(xiàn)有20個(gè)miRNA與PRCC患者的生存顯著相關(guān)(P＜0.01)(表1)。

表1 與PRCC患者生存相關(guān)的20個(gè)miRNATable 1 Twenty miRNAs related to the survival of PRCC patients

2.3 LASSO與Cox回歸分析獲取獨(dú)立影響預(yù)后的miRNA

首先,對與患者生存相關(guān)的20個(gè)miRNA進(jìn)行單因素Cox回歸分析,得到18個(gè)miRNA(P＜0.01)(表2)。然后,通過LASSO回歸對上述18個(gè)miRNA進(jìn)行再篩選(圖2),得到7個(gè)miRNA:mir-1293、mir-551a、mir-937、mir-543、mir-34a、mir-519a-1、mir-517c。進(jìn)一步對LASSO回歸分析結(jié)果進(jìn)行多因素Cox回歸分析,得到兩個(gè)獨(dú)立影響預(yù)后的miRNA,即mir-1293和mir-937(P＜0.05)(表3),它們的生存曲線見圖3。從圖3可知,在PRCC樣本中,mir-1293和mir-937高表達(dá)組的生存率均低于低表達(dá)組。

表2 基于單變量Cox回歸分析獲得的預(yù)后相關(guān)miRNATable 2 miRNAs related to prognosis based on univariate Cox regression analysis

圖2 18個(gè)關(guān)鍵miRNA的LASSO回歸圖每一條用不同顏色字母標(biāo)記的曲線代表一個(gè)miRNA系數(shù)的變化軌跡,縱坐標(biāo)是系數(shù)的值,下橫坐標(biāo)是logλ,上橫坐標(biāo)是此時(shí)模型中非零系數(shù)的個(gè)數(shù)。虛線指示了一個(gè)特殊的λ值:lambda.min,即為得到最小交叉驗(yàn)證均值誤差的λ值。a:mir-551a,b:mir-543,c:mir-937,d:mir-1293,e:mir-130b,f:mir-337,g:mir-517c,h:mir-802,i:mir-379,j:mir-1224,k:mir-323,l:mir-516a-1,m:mir-432,n:mir-519a-1,o:mir-134,p:mir-551b,q:mir-651,r:mir-34a。Fig.2 LASSO regression diagrams of 18 key miRNAsEach curve marked with a different color letter represents the trajectory of a miRNA coefficient.The vertical axis is the value of the coefficient.The lower horizontal axis is logλ,and the upper horizontal axis is the number of non-zero coefficients in the model at this time.The dotted line indicates a special lambda value(lambda.min),which gives the minimum cross-validation mean error.a:mir-551a,b:mir-543,c:mir-937,d:mir-1293,e:mir-130b,f:mir-337,g:mir-517c,h:mir-802,i:mir-379,j:mir-1224,k:mir-323,l:mir-516a-1,m:mir-432,n:mir-519a-1,o:mir-134,p:mir-551b,q:mir-651,r:mir-34a.

表3 7個(gè)miRNA的多因素Cox回歸分析結(jié)果Table 3 Multivariate Cox regression analysis of 7 miRNAs

圖3 獨(dú)立影響預(yù)后的miRNA的生存曲線圖以該miRNA在所有樣本中表達(dá)量中位數(shù)為基準(zhǔn),高于中位數(shù)為高表達(dá),低于中位數(shù)為低表達(dá)。Fig.3 Survival curves of miRNAs independently affecting prognosisBased on the median expression level of the miRNA in all samples,a level higher than the median is considered as high expression,while lower than the medium is considered as low expression.

2.4 獨(dú)立影響預(yù)后的miRNA靶基因及其表達(dá)分析

為了探究mir-1293、mir-937與PRCC的關(guān)系,從miRWalk3.0數(shù)據(jù)庫中,以score≥0.95為篩選條件,分別下載mir-1293的靶基因199個(gè)和mir-937的靶基因132個(gè)。

從GEO數(shù)據(jù)庫中獲取PRCC基因芯片數(shù)據(jù)GSE7023,使用limma包進(jìn)行差異表達(dá)分析。分析結(jié)果顯示,與癌旁非癌腎組織樣本相比,903個(gè)基因的表達(dá)存在顯著性差異[|log2(fold change)|≥1且P＜0.05],包括 638個(gè)(占70.7%)下調(diào)基因和 265個(gè)(占29.3%)上調(diào)基因。其中,mir-1293有CLMN(calmin)、DCN(decorin)、PODXL(podocalyxin-like)、CYP4A11(cytochrome P450 4A11)、KCNJ10 五個(gè)靶基因低表達(dá);mir-937有ST6GAL1(β-galactoside α2,6-sialyltranferase 1)、ATP6V1A(V-ATPase subunit A)、VAV3(vav guanine nucleotide exchange factor 3)三個(gè)靶基因低表達(dá)(圖4)。

圖4 PRCC組織中基因表達(dá)情況的火山圖相對于癌旁非腎癌組織,黑色節(jié)點(diǎn)為表達(dá)無差異的基因,右邊紅色節(jié)點(diǎn)為高表達(dá)的基因,左邊紅色節(jié)點(diǎn)為低表達(dá)的基因,CLMN、DCN、PODXL、CYP4A11、KCNJ10為 mir-1293靶基因,ST6GAL1、ATP6V1A、VAV3 為 mir-937靶基因。Fig.4 Volcano map of gene expression in PRCC tissuesCompared with normal tissues,black nodes represent genes without differential expression,red nodes on the right represent genes with up-regulated expression,and red nodes on the left represent genes with down-regulated expression.CLMN,DCN,PODXL,CYP4A11 and KCNJ10 are the target genes of mir-1293.ST6GAL1,ATP6V1A and VAV3 are the target genes of mir-937.

3 討論

PRCC是一種常見的腎臟惡性腫瘤,現(xiàn)在還沒有標(biāo)準(zhǔn)治療方案,且國內(nèi)研究也不多[21]。現(xiàn)有研究多為通過分析PRCC的臨床病理特征來研究其診療及預(yù)后因素等,在分子水平上開展預(yù)后標(biāo)志物的研究報(bào)道甚少。有研究發(fā)現(xiàn),miRNA與腎癌相關(guān),可調(diào)控多個(gè)癌基因表達(dá),從而影響腎癌的發(fā)生發(fā)展[22～23]。本研究通過miRNA的差異表達(dá)分析和生存分析發(fā)現(xiàn),mir-517c、mir-543等20個(gè)miRNA與PRCC病人總體存活率具有顯著相關(guān)性(P＜0.01)。其中,mir-1293和mir-937可能是PRCC預(yù)后的獨(dú)立影響因素。

mir-1293是一個(gè)新近引人關(guān)注的腫瘤相關(guān)miRNA,然而相關(guān)研究尚存在爭議。例如:Liu等[23]的最新研究報(bào)道,mir-1293可以通過靶向HAO2(hydroxy acid oxidase 2)促進(jìn)腎癌細(xì)胞的活力、侵襲和遷移,其高表達(dá)與預(yù)后不利相關(guān)。而另一些研究則表明,mir-1293具有抑制轉(zhuǎn)移性腎透明細(xì)胞癌的作用,其靶基因TIMP-1表達(dá)增加與腎透明細(xì)胞癌患者的預(yù)后較差有關(guān)[24];mir-1293能夠通過抑制BRD4(bromodomain-containing protein 4)和DNA修復(fù)基因,抑制異種移植小鼠模型中的體內(nèi)腫瘤生長,其可作為基于miRNA的抑瘤候選物[25]。

本研究提示mir-1293高表達(dá)與PRCC不良預(yù)后相關(guān),其中CLMN和DCN可能是mir-1293影響PRCC預(yù)后的靶基因。CLMN可通過下調(diào)細(xì)胞周期蛋白D1表達(dá),影響細(xì)胞G1/S的轉(zhuǎn)變,從而抑制細(xì)胞周期進(jìn)程,是潛在的腫瘤抑制基因[26],其高表達(dá)可抑制神經(jīng)母細(xì)胞瘤細(xì)胞增殖[27],延長乳腺癌患者的無復(fù)發(fā)生存期[28]。DCN可影響細(xì)胞增殖、擴(kuò)散、遷移和分化,其表達(dá)在多種癌組織中降低,并與腫瘤大小明顯相關(guān)。相關(guān)研究報(bào)道,DCN過表達(dá)可抑制腎癌細(xì)胞的增殖和轉(zhuǎn)移,提示其是一種具有潛在臨床價(jià)值的腎細(xì)胞癌抑制因子[29～30]。

諸多研究表明,mir-937具有促癌、抑癌雙重性,具體作用與所處組織環(huán)境有關(guān)。在肺癌組織中,mir-937上調(diào),其通過INPP4B(inositol polyphosphate 4-phosphatase type Ⅱ)促進(jìn)細(xì)胞增殖[31]。臨床分析表明,mir-937表達(dá)與結(jié)腸癌腫瘤淋巴結(jié)轉(zhuǎn)移和TNM分期有關(guān),mir-937高表達(dá)的患者預(yù)測總體生存率低[32～33]。此外,有研究在乳腺癌組織中檢測到mir-937的高表達(dá),且mir-937高表達(dá)與癌癥患者存活率降低相關(guān)[34]。也有研究發(fā)現(xiàn),mir-937是腫瘤的抑制因子。例如:mir-937的過表達(dá)可通過靶向FOXL2(forkhead box L2)使PI3KAkt信號通路失活,從而抑制胃癌細(xì)胞的增殖和轉(zhuǎn)移[35];在乳腺癌中,mir-937可通過抑制FOXQ1(forkhead box Q1)表達(dá)起到抑癌作用[36]。

本研究的分析結(jié)果提示,mir-937在PRCC組織中高表達(dá)與患者不良預(yù)后呈顯著相關(guān)性,可作為PRCC預(yù)后的生物標(biāo)志物。基因表達(dá)芯片數(shù)據(jù)的驗(yàn)證結(jié)果顯示,mir-937的潛在靶基因VAV3、ST6GAL1和ATP6V1A在PRCC中下調(diào)。上述3個(gè)mir-937潛在靶基因均與腫瘤關(guān)系密切。其中,VAV3是鳥嘌呤核苷酸交換因子,其高表達(dá)與促腫瘤生長、侵襲、轉(zhuǎn)移和血管生成有關(guān)[37～38];ST6GAL1被認(rèn)為對膀胱癌具有抑癌作用[39],但其上調(diào)可誘導(dǎo)胃癌、結(jié)腸癌、肝癌、前列腺癌、卵巢癌、乳腺癌和宮頸癌,被認(rèn)為具有促進(jìn)和抑制癌癥進(jìn)展的雙重作用[40];ATP6V1A的高表達(dá)有助于胃癌的良好預(yù)后[41],但對大腸癌的預(yù)后存在不利影響[42],其對腫瘤的作用機(jī)制依然存在爭議。

綜上所述,本研究挖掘出了PRCC中差異表達(dá)的miRNA,并獲得預(yù)后的獨(dú)立影響因素 mir-1293與mir-937。PRCC組織中的mir-1293與mir-937高表達(dá)可能分別通過抑制其靶基因表達(dá),促進(jìn)腫瘤發(fā)生發(fā)展,影響治療預(yù)后,故有望作為PRCC預(yù)后不良的生物標(biāo)志物,供臨床上對患者的治療和預(yù)后進(jìn)行判斷。