基于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预后不良的生物标志物,供临床上对患者的治疗和预后进行判断。