基于lncRNA-mRNA共表达网络筛选胃癌分期相关的lncRNA*

哈尔滨医科大学卫生统计学教研室(150081) 丛雨欣 张秋菊 田 伟 张 奇 李 称 赵 敏 刘美娜

【提 要】 目的 构建lncRNA-mRNA共表达网络，探索与胃癌分期相关的lncRNA及其调控关系，为研究胃癌的进展机制及寻找胃癌治疗的潜在靶点提供依据。方法 利用TCGA数据库，收集胃癌RNA-seq数据及临床信息数据；采用表达数量性状位点(eQTL)分析与加权基因共表达网络分析(WGCNA)相结合的方法，构建lncRNA-mRNA共表达网络模块，结合临床信息筛选与胃癌分期相关的模块；采用Kruskal-Wallis秩和检验筛选模块内胃癌不同分期差异表达的lncRNA。结果 获得286例胃癌组织和30例癌旁对照组织样本的RNA-seq数据；eQTL分析得到5118对顺式作用和1 953 109对反式作用lncRNA-mRNA；2 999个lncRNA和3 884个mRNA纳入WGCNA，产生25个共表达模块，其中与胃癌分期高度相关的模块有3个；模块midnightblue、orange内18个枢纽lncRNA中有14个lncRNA在胃癌不同分期差异表达。结论 本研究筛选出14个与胃癌分期相关的lncRNA，这14个lncRNA可能通过调控mRNA的表达影响胃癌的进展，分析其对应的网络调控关系，为研究lncRNA-mRNA调控机制及探索胃癌治疗靶点提供了参考和方向。

lncRNA(long non-coding RNA)是一类转录长度超过200nt的非编码序列，通常由RNA聚合酶II转录形成。由于不编码蛋白质，lncRNA起初被认为是没有功能的转录垃圾。近年来发现，lncRNA在表观遗传学、转录及转录后水平等方面调控基因的表达，如干扰邻近基因的表达、作为共因子调节转录因子的活性、在转录后水平可与mRNA形成双链复合物调控基因表达水平。更值得关注的是，lncRNA在细胞异常调节、诱导肿瘤发生过程中起关键作用[1]。多项研究表明，lncRNA 与胃癌的发生、发展、转移和预后密切相关[2-3]。因此，探索与胃癌分期相关的lncRNA，明确其对靶基因的调控关系对胃癌机制研究及治疗具有重要意义。

加权基因共表达网络分析(weighted gene co-expression network analysis，WGCNA)作为一种高通量数据挖掘算法，已广泛应用于筛选疾病相关基因的研究领域，如利用WGCNA构建与胃癌进展相关的lncRNA网络[4]。本研究在此基础上将表达数量性状位点(expression quantitative trait loci，eQTL)分析与WGCNA相结合，收集胃癌RNA-seq数据及临床信息数据，利用eQTL筛选出的相关程度高的lncRNA-mRNA共表达对，对应的lncRNA和mRNA纳入WGCNA，构建共表达网络模块，获得与胃癌分期相关的lncRNA并分析其调控关系，为胃癌进展及治疗靶点研究奠定基础。

数据与方法

1.数据来源

在TCGA(the cancer genome atlas)数据库中，纳入286例胃癌组织和30例癌旁对照组织样本的RNA-seq reads count数据。收集匹配的临床数据，包括年龄、性别、肿瘤TNM分期等信息。

2.分析方法

(1)预处理方法

RNA-seq数据的预处理：包括提取样本信息、构建基因表达矩阵、将探针名转化为基因名，最终获得行名为样本名、列名为基因名的矩阵；利用edge R进行差异分析，设置log2 fold change界值为1，FDR界值为0.05。

(2)表达数量性状位点分析

表达数量性状位点分析由数量性状位点分析(quantitative trait loci，QTL)发展而来，QTL指定位控制数量性状的基因在基因组中的位置[5]，eQTL是将每个基因的表达水平作为数量性状进行分析[6]。通过全基因组mRNA表达量测序得到特定组织样本的基因表达量，以全部SNP为自变量、每种mRNA表达量为因变量进行线性回归，得到每一个SNP位点和mRNA表达量之间的关系。由于eQTL分析中表型数据不限于离散型数据，因此eQTL可作为lncRNA-mRNA共表达分析方法[7]。本研究将SNP数据替换为lncRNA表达量数据，分析mRNA与lncRNA表达量的相关性，获得mRNA-lncRNA共表达对。

eQTL分析的优点是可以区分顺式作用关系和反式作用关系，有利于具体分析lncRNA对mRNA的调控作用。顺式调控[8]是指lncRNA对染色体临近位置(小于1Mb)的mRNA的表达调控，位于编码蛋白上下游的lncRNA可能与启动子或共表达基因的其他顺式作用元件有交集，从而在转录或者转录后水平对基因的表达进行调控。反式调控是指对不同染色体或者染色体远端位置的mRNA进行表达调控，调控关系不受空间距离的限定，因此trans-eQTL数量庞大，且可能存在假阳性相关。基于此，本研究选取eQTL分析获得的全部顺式作用共表达对和前5%具有显著性的反式作用共表达进行WGCNA分析，构建lncRNA与mRNA共表达网络。

(3)加权基因共表达网络分析

加权基因共表达网络分析是一种在高通量基因表达数据中，利用系统生物学思想，计算基因表达相关性，构建基因共表达模块，进而发现具有生物学意义模块的高通量数据挖掘算法[9]。该算法首先计算各基因间Pearson相关系数，构建加权邻接矩阵：

aij=power(Sij，β)≡|Sij|β

其中aij代表基因i与基因j的邻接系数；Sij代表基因i与基因j的Pearson相关系数。β为软阈值，软阈值的确定应满足使共表达网络服从无尺度网络分布，即出现连接度为k的节点的对数lgk与该节点出现的概率的对数lg[P(k)]呈负相关，且R2应>0.8。

基因间的邻接矩阵转换为拓扑矩阵(topological matrix，TOM)，将某个基因与所有直接、间接相关的基因间的关系联系起来，基于拓扑重叠性计算基因与基因间的相异度，根据相异度利用动态剪切树法划分基因的功能模块，完成共表达网络的基础构建。继而计算模块特征值(即该模块的第一主成分)与临床表型信息的相关系数，筛选出与胃癌分期相关的模块作为关键模块，分析关键模块内的lncRNA及其共表达的mRNA。根据网络模块节点中度的分析，确定模块中的枢纽lncRNA，分析其在胃癌分期是否差异表达。由于目前已知功能的lncRNA微乎其微，这种分析策略有效缩小候选lncRNA的范围，有助于寻找与胃癌分期相关的lncRNA。

(4)统计分析

所有分析通过R 3.6.1和Cytoscape 3.7.2实现。①edgeR包筛选差异表达基因；②筛选出的基因通过MatrixEQTL包进行eQTL分析，获得lncRNA-mRNA共表达对；③利用WGCNA包实现共表达网络的构建并筛选与胃癌分期相关的关键模块；④将模块内基因共表达网络的权重信息导入Cytoscape 3.7.2软件，筛选模块枢纽基因并绘制基因共表达网络图；⑤对枢纽lncRNA进行Kruskal-Wallis秩和检验，分析各枢纽lncRNA在不同胃癌分期是否差异表达；⑥模块内基因映射至在线网站DAVID(http：//david-d.ncifcrf.gov/)中，进行GO和KEGG富集分析。

结 果

1.预处理结果

差异分析获得4 767个差异表达的mRNA，其中2 466个mRNA上调，2 301个mRNA下调；3 542个差异表达的lncRNA，其中2 767个lncRNA上调，775个lncRNA下调。对差异基因进行火山图的可视化，见图1。

图1 差异表达基因火山图

2.构建胃癌lncRNA-mRNA共表达网络

以年龄、性别为协变量对差异基因进行eQTL分析，获得5 118对顺式作用共表达对和1 953 109对反式作用共表达对，取全部的顺式作用共表达对和前5%具有显著性的反式作用共表达对，去除重复基因，最终获得2 999个lncRNA和3 884个mRNA进入下一步的WGCNA分析。

选取软阈值β=4，通过动态剪切树法进行模块初步识别并合并相似模块，设置每个基因网络模块最少的基因数目为50，模块合并阈值为0.3，最终得到25个基因网络模块，灰色模块是无法聚类到其他任何模块的基因集合，见图2。

图2 加权基因共表达网络模块聚类图

根据各模块的特征向量分别计算每个模块与胃癌分期的相关性，绘制模块样本性状相关性热图，见图3。结果显示，midnightblue、orange、yellow、purple、blue、royalblue、green共七个模块与胃癌分期的相关性有统计学意义(P<0.05)。

图3 临床信息与模块相关性热图

3.筛选与胃癌分期相关lncRNA

选择与胃癌分期相关性较强的前三个模块，midnightblue、orange、yellow模块，各模块按度的大小排序前50的基因为枢纽基因，见图4。三个模块中分别包含10个、8个、3个枢纽lncRNA和40个、42个、47个枢纽mRNA。秩和检验结果获得14个与胃癌分期相关lncRNA：midnightblue模块有AP002954.4、AC002331.1、LINC01272、RP11-44K6.4、RP11-638I2.9、AC069363.1、LINC01094，orange模块有RP11-443C10.1、RP11-13P5.2、AC093850.2、RP11-576I22.2、AC007750.5、RP11-95H3.1、RP11-867G23.1；其中AP002954.4在T3、T4期表达水平相较于T1期上调，AC002331.1、RP11-44K6.4在T4期表达水平相较于T1、T2期上调，AC069363.1在T4期表达水平相较于T1期上调，RP11-638I2.9在T4期表达水平相较于T2期上调，LINC01272在T2期表达水平相较T1期上调，LINC01094、RP11-443C10.1、RP11-13P5.2、AC093850.2、RP11-576I22.2、AC007750.5、RP11-95H3.1、RP11-867G23.1在T2、T3、T4表达水平相较于T1期上调，差异均具有统计学意义，见表1。与14个lncRNA共表达的mRNA共有59个，其中有56个mRNA秩和检验有统计学意义。

表1 枢纽lncRNA秩和检验结果

图4 midnightblue(a)、orange(b)、yellow(c)模块的共表达网络图

4.富集分析

lncRNA-mRNA共表达模块中的59个mRNA进行GO和KEGG分析，结果见表2。GO分析结果显示，模块基因主要涉及到免疫应答(immune response)、防御反应(defense response)、趋化作用(chemotaxis)和炎性反应(inflammatory response)等生物过程；KEGG分析结果显示，模块基因主要富集于细胞因子-受体相互作用通路(cytokine-cytokine receptor interaction)、Toll样受体信号通路(toll-like receptor signaling pathway)、趋化因子信号通路(chemokine signaling pathway)、黏着斑通路(focal adhesion)等，这些通路与机体炎性反应有关。

表2 富集分析结果

讨 论

本研究首次将WGCNA与eQTL相结合应用在癌症数据中，筛选与胃癌分期相关的lncRNA。由于目前大部分lncRNA功能未知，通过分析lncRNA与mRNA表达相关性，利用已知功能的mRNA推导lncRNA的功能成为一种重要的研究策略。相较于既有的单纯通过WGCNA构建与胃癌进展相关的lncRNA网络的研究[4]，本研究结合eQTL的优势在于，一方面筛选出相关程度较高的lncRNA与mRNA进行共表达网络分析，提高lncRNA功能推断的准确性；另一方面将lncRNA与mRNA的共表达区分为顺式共表达对与反式共表达对，有利于后续进一步研究生物学上lncRNA与mRNA的调控关系。本研究获得14个与胃癌分期相关的lncRNA，总体上随着胃癌分期的进展，lncRNA表达水平呈上升趋势。与14个lncRNA共表达的59个mRNA中56个mRNA在不同分期表达水平差异有统计学意义，这说明eQTL与WGCNA相结合的方法能够有效筛选出表达相关性高的lncRNA-mRNA，筛选出的lncRNA可能通过调控mRNA的表达在胃癌的发展进程中具有关键作用，可作为潜在的生物标志物。研究发现，midnightblue模块的枢纽lncRNA LINC01272在肺鳞癌与癌旁组织中显著差异表达，且可以作为诊断肺鳞癌早期与晚期的标志物[10]；orange模块的枢纽lncRNA AC093850.2在肺鳞癌、肺腺癌、乳腺癌及胰腺癌中表达水平均高于癌旁组织[11]。目前仍未有研究发现这些枢纽lncRNA与胃癌分期有关，本研究结果为深入研究胃癌分期相关lncRNA提供了参考和依据。

lncRNA通过调控mRNA的转录、稳定性等多种方式影响靶基因的表达，lncRNA-mRNA共表达网络是研究lncRNA功能和调控机制的重要方式。GO分析和KEGG分析结果显示lncRNA-mRNA共表达网络中的mRNA主要参与免疫应答、防御反应、炎性反应等生物过程，可以推断这些枢纽lncRNA可能通过调控mRNA影响机体的炎性反应、免疫应答等，在一定程度影响胃癌的进展过程。已有多项研究发现，非编码RNA可以通过影响肿瘤和免疫细胞中免疫调节分子的表达影响抗肿瘤免疫应答、参与炎性反应过程[12]。midnightblue模块中的枢纽lncRNA LINC01272已被证实与炎性反应有关，Wang Sen等的研究发现LINC01272在炎症性肠病患者中的表达水平显著高于健康对照，是炎症性肠病的潜在诊断标志物[13]。

综上，利用eQTL和WGCNA相结合的方法构建lncRNA-mRNA共表达网络，探索与胃癌分期相关的lncRNA及其调控关系。筛选出的14个与胃癌分期相关的lncRNA可为研究胃癌的进展机制提供依据，其对应的网络调控关系为后续更深入研究lncRNA-mRNA调控机制及探索胃癌治疗靶点提供了参考和方向。