范敬炎,韩欣桐,丁家安,王树急,杜晶晶,朱玉峰
(1.锦州医科大学;2.锦州医科大学附属第一医院,辽宁 锦州 121000)
结肠癌(colon cancer)是一种异质性疾病,给患者和医疗系统造成了巨大危害[1]。因此,在分子水平上更深入地了解其调控机制,确定可靠的预后生物标志物对于结肠癌患者的生存预测及个体化治疗具有重要的作用。
长链非编码RNA(LncRNA)的长度超过200个核苷酸,并且在特定组织以及不同类型肿瘤中差异表达,尽管lncRNA不编码蛋白质,但它们仍然具有许多功能,与肿瘤发生发展密切相关[2]。而且lncRNA在体液中易于检测,使其成为理想的临床预后判断的生物标志物[3]。
随着各种生物信息学分析工具以及公共数据库的出现,从高通量数据中鉴定关键基因变得越来越容易。其中,加权基因共表达网络分析(WGCNA)是寻找高度相关基因模块的有力工具,已被广泛用于鉴定候选生物标志物。在这项研究中,我们从公共数据库下载结肠癌患者的RNA表达数据和临床信息,并建立共表达网络以挖掘可以判断结肠癌患者预后的lncRNA。使用独立验证数据集,包括结肠癌样本组织和临床生存信息,验证lncRNA表达水平和结肠癌患者预后的关系。同时深入研究候选lncRNAs在CRC中的表达和功能,以提高我们对结肠癌发生发展的分子机制的理解,并为CRC诊疗提供候选靶标。
结肠癌的基因表达原始数据和临床信息从TCGA数据库中获取,该数据库中包含了467份结肠癌样本,包括268份I~II期和199份III~IV期结肠癌样本。根据临床信息整理相应患者的生存状态和生存时间。本研究以TCGA数据集为训练集,构建加权基因共表达网络,识别结肠癌预后相关基因。使用数据集GSE39582作为测试集进行生存分析来验证我们的结果。数据集GSE39582下载自GEO数据库,该数据集包括566例结肠癌样本,其中TNM分期I期33例,II期264例,III期205例,IV期60例,总生存数据556例,无复发生存数据519例。使用GSE20916数据集验证候选预后分子在结肠癌中的表达水平,该数据集包括101例结肠癌及35例正常结肠组织。使用GSE41568数据集观察候选预后分子在原位结肠癌和转移性结肠癌组织中的表达水平,该数据集包括39例原位结肠癌和94例转移性结肠癌。数据集GSE39582、GSE20916和GSE41568的芯片均在同一平台(GPL570)上进行检测。
1.2.1 筛选差异表达的LncRNA 利用R语言软件的Limma函数包对原始表达式数据进行预处理。从Gencode数据库(https://www.gencodegenes.org)获取人类基因组(hg38)和相关注释文件(31版)。该注释文件用于识别LncRNA。基因类型为“lincrna”、“antisense”、“processed transcript”、“sense_intronic”、“TEC”、“bidirectional promoter lncRNA”、“sense_overlapping”、“macrolncRNA”或“non coding”的分子定义为LncRNA[4]。当多个探针与一个相同的LncRNA匹配时,我们取表达值的平均值[5]。选择FDR<0.05,Log2(FC)>2,P<0.05作为筛选差异基因的标准。
1.2.2 加权基因共表达网络分析 剔除TCGA数据集中异常离群样本,然后逐步进行网络构建和模块聚类。为了构建无尺度基因共表达网络,我们使用了R语言中的WGCNA函数包[6]。首先对所有的基因对进行皮尔逊相关矩阵分析。然后利用幂函数amn=|cmn|β(其中amn是m和n基因之间的邻接,cmn是m和n基因之间的皮尔逊相关性)构造加权邻接矩阵。参数β是软阈值参数,经过计算选择β=14(无尺度r2=0.91)以确保构建无尺度网络[7]。然后,将邻接矩阵转化为拓扑重叠矩阵,由此得到的拓扑重叠是基于两个基因间共表达关系的一种有生物学意义的基因相似性度量[8]。采用动态树切割法进行模块识别,用颜色命名模块,采用Pearson相关分析计算各个模块与TNM分期、总生存期、无复发生存期等临床特征的相关性和P值,相关系数最高的基因模块内的LncRNA作为候选预后分子标志物纳入下一步分析。
1.2.3 筛选预后风险LncRNA及功能预测 使用GSE39582数据集对候选LncRNA进行生存分析,绘制Kaplan-meier生存曲线,所有P值小于0.05的LncRNA作为潜在的预后分子标志物。在此基础上,利用另外的独立数据集GSE20916进行分析,研究正常结肠标本和结肠癌标本中LncRNA表达水平差异。利用GSE41568在原位结肠癌和转移结肠癌样本之间验证LncRNA表达水平差异。为了探索这些LncRNA影响相关临床特征的潜在机制,将候选LncRNA上传到LncACTdb2.0数据库中,提取lncRNA-mRNA共表达系数绝对值大于0.6且P<0.05的网络,并将其中的mRNA定义为该lncRNA的靶基因。对靶基因进行GO功能和KEGG信号通路富集分析,将P<0.01和FDR<0.01设置为筛选标准[9]。
使用R语言的 survival 函数包进行Kaplan-Meier生存分析,生存率比较采用log-rank检验法,P< 0.05 表示差异具有统计学意义。
在TCGA数据库中下载结肠癌数据集,共发现差异表达基因4227个,其中上调基因1899个,下调基因2328个。差异表达基因中,mRNA及LncRNA的数量分别为1817个(上调755个,下调1062个)及2410个(上调1144个,下调1266个)。差异基因的表达数据谱及其临床信息将用于进行接下来的加权基因共表达网络分析。
基于加权基因共表达网络分析将所有差异表达的基因分成12个模块,见图1A。为了推测这些基因的临床意义,将基因模块与结肠癌患者的临床信息相结合。肿瘤组织分期及患者生存资料作为选择功能模块的重要评估指标。结果显示粉色模块(pink module)的特征值与结肠癌患者的总生存期具有较高的相关性(r=0.68,P=0.0001),见图1B。为了探索预后分子标志物,本研究选择粉色模块中的LncRNA进行下一步的生存分析。
A:动态分层聚类图;B:特征基因模块和结肠癌不同临床信息相关性的热图
位于粉色模块中的9个LncRNA(LINC00473,FAM228B,CTC-277H1.6,LOC642846,CHKB-AS1,LRRC75A-AS1,LINC00299,LOC400684,LINC01021)被视为候选的结肠癌预后分子标志物。下载GEO数据库中的GSE39582数据集,利用 R 软件 survival 软件包及GSE39582数据集中的临床资料对9个候选LncRNA进行生存分析。根据每个候选lncRNA的表达值的中位数将结肠癌患者分成两组, Kaplan-Meier生存分析结果显示LINC01021高表达的患者的总体存活率更低(P=0.001),复发风险更高(P=0.001),见图2。其余LncRNA表达水平与结肠癌患者临床预后的关系,差异无统计学意义。
使用另外的独立数据集GSE20916对候选lncRNA在结肠癌中的表达水平进行验证,LINC01021在结肠癌中表达水平高于正常结肠组织,两组LINC01021表达水平比较,差异有统计学意义(t=9.549,P<0.01),见图3A。此外,在数据集GSE41568中,LINC01021在转移结肠癌中的表达水平高于原位结肠癌,差异有统计学意义(t=4.927,P<0.01),见图3B。使用LncACTdb2.0数据库预测LINC01021下游靶基因,共137个基因,见图4。为了进一步推测LINC01021的功能,对下游靶基因进行功能富集分析和信号通路富集分析,GO分析发现lncRNA-mRNA共表达网络中的mRNA主要参与细胞增殖负调控,调节细胞周期素依赖蛋白激酶活性,调节细胞周期等,KEGG分析发现lncRNA-mRNA共表达网络中的mRNA主要参与PI3K/AKT信号通路激活和CTNNB1磷酸化级联反应等,结果见图5。
图5 GO功能富集分析和KEGG信号通路富集分析
A:LINC01021在结肠癌组织及正常组织中的表达水平;B:LINC01021在原位结肠癌组织及转移结肠癌组织中的表达水平
图2 结肠癌患者差异表达的LINC01021生存曲线分析
红色代表上调,绿色代表下调;圆形代表 mRNA,三角形代表 lncRNA,图形越大表示共表达系数绝对值越大
在本研究中,利用TCGA数据库的数据集,我们首先通过加权基因共表达网络分析筛选了与结肠癌预后相关的基因模块,然后利用GEO数据库中的独立数据集对模块中的基因进行了生存分析,确定了长链非编码RNA(LINC01021)与结肠癌预后相关,高表达LINC01021的结肠癌患者预后差。通过一系列的生物信息学分析,发现LINC01021下游的靶基因主要参与结肠癌细胞周期的调控,以及通过激活PI3K/AKT信号通路和参与CTNNB1磷酸化级联反应等过程影响结肠癌的发生发展,表明LINC01021可能在结肠癌的发展中发挥致癌作用。
LINC01021是长度为1112个碱基的长链非编码RNA,定位于人染色体5p14.1。既往研究报道lncRNA LINC01021是p53的靶基因,能够与p53直接结合[10]。使用CRISPR /Cas9方法敲除p53结合位点的启动子序列可以降低结肠癌HCT116细胞系中LINC01021表达,增加结肠癌细胞对多柔比星和5-氟尿嘧啶的化疗敏感性[11]。国内研究报道LINC01021在食管癌组织和食管癌细胞系中的高表达,LINC01021 基因通过影响上皮间质转化促进食管癌细胞侵袭转移[12]。我们的结果发现LINC01021在转移结肠癌中的表达水平高于原位结肠癌,而结肠癌的侵袭转移与上皮间质转化密切相关,这是一个值得探索的研究方向。我们的研究结果还发现LINC01021高表达的结肠癌患者的总体存活率更低,复发风险更高。这些结果说明LINC01021表达的定量可能具有重要的临床预后价值。
我们通过生物信息学分析发现LINC01021下游靶基因发挥的作用机制主要是影响细胞周期调控,以及通过激活PI3K/AKT信号通路和参与CTNNB1磷酸化级联反应。PI3K/AKT信号通路有广泛的生物活性,可促进细胞增殖,在恶性肿瘤内往往呈过度激活的状态,是最主要的抑制细胞凋亡的信号途径[13]。CTNNB1磷酸化的降低往往是由于Wnt信号通路活化,随着信号通路活化和CTNNB1蛋白在核内积累,将导致结肠癌发生[14]。因此,本研究结果对未来探索LINC01021调控结肠癌发生的分子机制奠定了一定的理论基础,具有非常重要的指导意义。当然,需要进一步的实验来验证LINC01021的临床和生物学功能。
总之,我们采用WGCNA等生物信息学方法从TCGA数据库中研究了结肠癌患者的RNA序列和临床数据。我们得出结论是LINC01021与结肠癌患者的预后相关,高表达LINC01021的结肠癌患者预后差。LINC01021有可能成为一种新的预后指标,有助于结肠癌患者个性化治疗及临床预后判断。