基于miRNA的乳頭狀腎細胞癌預后指標的生物信息學研究

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

(福建中醫藥大學a.中西醫結合學院;b.中西醫結合慢性病研究福建省高校重點實驗室,中國福建 福州 350122)

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

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

1 數據與方法

1.1 miRNA與基因芯片表達數據下載

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

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

1.2 miRNA表達的差異分析

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

1.3 生存分析

為分析上述差異表達miRNA與患者生存之間的相關性,采用survival包、survminer包進行批量Kaplan-Meier生存分析,找出與PRCC患者生存顯著相關的miRNA,篩選標準為P＜0.01。

1.4 LASSO及Cox回歸分析

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

1.5 獨立影響預后的miRNA靶基因預測及表達分析

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

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

2 結果

2.1 miRNA表達的差異分析

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

圖1 PRCC組織中miRNA表達情況的火山圖相對于癌旁非癌腎組織,黑色節點表示表達無差異的miRNA(356個),紅色節點表示表達上調的miRNA(117個),綠色節點表示表達下調的miRNA(114個)。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語言對上述差異表達的miRNA進行生存分析,發現有20個miRNA與PRCC患者的生存顯著相關(P＜0.01)(表1)。

表1 與PRCC患者生存相關的20個miRNATable 1 Twenty miRNAs related to the survival of PRCC patients

2.3 LASSO與Cox回歸分析獲取獨立影響預后的miRNA

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

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

圖2 18個關鍵miRNA的LASSO回歸圖每一條用不同顏色字母標記的曲線代表一個miRNA系數的變化軌跡,縱坐標是系數的值,下橫坐標是logλ,上橫坐標是此時模型中非零系數的個數。虛線指示了一個特殊的λ值:lambda.min,即為得到最小交叉驗證均值誤差的λ值。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個miRNA的多因素Cox回歸分析結果Table 3 Multivariate Cox regression analysis of 7 miRNAs

圖3 獨立影響預后的miRNA的生存曲線圖以該miRNA在所有樣本中表達量中位數為基準,高于中位數為高表達,低于中位數為低表達。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 獨立影響預后的miRNA靶基因及其表達分析

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

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

圖4 PRCC組織中基因表達情況的火山圖相對于癌旁非腎癌組織,黑色節點為表達無差異的基因,右邊紅色節點為高表達的基因,左邊紅色節點為低表達的基因,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是一種常見的腎臟惡性腫瘤,現在還沒有標準治療方案,且國內研究也不多[21]。現有研究多為通過分析PRCC的臨床病理特征來研究其診療及預后因素等,在分子水平上開展預后標志物的研究報道甚少。有研究發現,miRNA與腎癌相關,可調控多個癌基因表達,從而影響腎癌的發生發展[22～23]。本研究通過miRNA的差異表達分析和生存分析發現,mir-517c、mir-543等20個miRNA與PRCC病人總體存活率具有顯著相關性(P＜0.01)。其中,mir-1293和mir-937可能是PRCC預后的獨立影響因素。

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

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

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

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

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