基于生物信息学分析LncRNA,AC005479.2在甲状腺乳头状癌中的表达及相关分析

梁新科，王鹏飞，霍胜男，李建英，侯庆田

(河北省邯郸市邯钢医院甲状腺科，河北 邯郸 056000)

甲状腺乳头状癌(papillary thyroid carcinoma，PTC)是最常见的甲状腺癌类型，占甲状腺恶性肿瘤的70%～85%，并且全球发病率逐年增加[1]。虽然甲状腺乳头状癌普遍预后良好，但其肿瘤的复发和转移亦严重影响患者的生存率。近年来，基因检测得以广泛普及，可以研究基因表达与疾病之间的关系，为肿瘤机制的研究提供了明确的方向。本研究通过生物信息学方法筛选分析PTC与健康人群的甲状腺组织的差异表达基因，寻找其发挥作用的相关信号通路，探索PTC的发病机制，为PTC诊断及治疗提供新的研究方向。

1.1一般资料 选取TCGA数据库中2000年1月—2020年12月收录并具有完整随访信息的甲状腺乳头状癌样本510例与癌旁组织样本58例。纳入标准：①肿瘤组织样本来自甲状腺乳头状癌患者，对照组为正常组织；
②样本须同时含有基因表达谱数据与临床信息。

1.2TCGA数据的分析 使用Rstudio中的“Limma”软件包对TCGA数据进行差异分析，寻找甲状腺乳头状癌组织与正常组织之间是否存在差异长链非编码RNA(long non-coding RNA, lncRNA)，并绘制火山图和热图。差异lncRNA筛选标准为：|log2FC|>1，FDR<0.05。

1.3WGCNA构建及临床相关性分析 使用“AnnotationDbi”R包将该数据中的探针ID转换成基因符号以进行进一步分析。使用R软件中的“FlsahClust”软件包进行样本层次聚类(implementation of optimal hierarchical clustering)，过滤变异最大的25%基因。利用WGCNA软件包中的“Pick-Soft-Threshold”函数调节参数β的权重。使用WGCNA软件包将数据中有相关性和相邻关系的lncRNA计算成为拓扑重叠矩阵(topological overlap matrix，TOM) ，然后计算相应的相异度(1-TOM) 。使用1-TOM作为距离度量，使用动态树切割法进行分层聚类以识别模块，将高度相似的模块通过聚类标记并合并。最小基因模块>100，合并相似模块的阈值=0.1用于搜索在心肌甲状腺乳头状癌中起重要作用的模块。获得基因模块后根据模块化特征将甲状腺乳头状癌患者的临床特征和模块内的基因联合分析，探究模块内基因的生物学意义。使用相关性热图探究基因与模块的相关性，并进行聚类分析。此外，WGCNA算法还对在同一模块中基因进行相互作用网络的预测。Cytoscape是一种分析软件，可以为生物学家提供生物网络分析和二维可视化。

1.4模块中差异lncRNA的筛选 Funrich是一种生物学分析软件。本研究中，使用维恩图，将blue模块中的lncRNA与差异lncRNA取交集，得到blue模块中差异lncRNA。

1.5统计学方法 本研究在Rstudio中使用“survival”包进行单因素和多因素Cox回归分析，进行模型变量筛选。受试者工作特征(receiver operating characteristic curve， ROC)曲线用于评估基因的诊断价值。P<0.05为差异有统计学意义。

1.6功能富集分析 选取甲状腺乳头状癌样本，将其分为AC005479.2低表达组和AC005479.2高表达组，随后导入GSEA中进行基因本体论(Gene Ontology，GO)分析和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Geomes，KEGG)分析。

2.1差异lncRNA的筛选 经TCGA数据库筛选，共获得差异表达lncRNA58个，其中包括高表达基因36个与低表达基因22个(P<0.05)，见图1A。LncRNA表达情况见图1B。

图1 甲状腺乳头状癌基因表达差异A.甲状腺乳头状癌基因表达差异火山图;B.甲状腺乳头状癌LncRNA表达情况Figure 1 Differences in gene expression of thyroid papillary carcinoma

2.2加权基因共表达网络构建及临床相关性分析 “pickSoftThreshold”函数的结果显示，当权重参数β=12时，log(k)与log[p(k)]之间相关系数的平方>0.9。选取软阈值β=12构建WGCNA，见图2A。使用动态树切割法识别基因模块，设定模块内最低基因个数为100，最终得到11个相应的模块，见图2B。结果显示，模块blue与甲状腺乳头状癌显著相关，见图2C。基因与模块的相关性分别展现在图2D与图2E中。预测出的Blue模块中基因之间的相互作用见图2F。

图2 加权基因共表达网络构建A.构建WGCNA;B.基因模块;C.蓝色模块与甲状腺乳头状癌显著相关;D.基因与模块的相关性;E.基因与模块的相关性;F.蓝色模块中基因之间的相互作用Figure 2 Construction of weighted gene co-expression networks

2.3模块中差异基因的筛选 通过维恩图将得到甲状腺乳头状癌组织与正常组织之间的58个差异lncRNA与blue模块中174个模块lncRNA取交集，得到蓝色模块内差异lncRNA共有34个，见图3。

图3 蓝色模块中的差异lncRNAFigure 3 Differential lncRNA in blue module

2.4预后模型的构建与验证 对蓝色模块内34个差异lncRNA进行单因素和多因素Cox回归分析(图4A、B)，得到AC005479.2。随后对AC005479.2绘制ROC曲线，曲线下面积(area under curve，AUC)值为0.838(图4C)。

图4 差异lncRNA的回归分析A.单因素Cox回归分析;B.多因素Cox回归分析;C.ROC曲线Figure 4 Regression analysis of differential lncRNAs

2.5基因富集分析 通过GO和KEGG分析，揭示AC005479.2的潜在通路与功能，结果显示AC005479.2通过脂肪细胞因子信号途径、丝裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)信号通路、钙黏蛋白结合参与细胞黏附、对凋亡信号通路的正调控等功能参与甲状腺乳头状癌的病理过程，见图5。

图5 富集分析显示AC005479.2的潜在通路与功能A.AC005479.2低表达组的GO分析;B.AC005479.2高表达组的GO分析;C.AC005479.2低表达组的KEGG分析;D.AC005479.2高表达组的KEGG分析Figure 5 Enrichment analysis showing the potential pathways and functions of AC005479.2

甲状腺癌是最常见的内分泌肿瘤，占人类所有恶性肿瘤的0.5%～1%[2]。PTC被认为是复发率最高的甲状腺癌类型，但同时也被认为是具有良好预后的惰性癌症，占已分化甲状腺癌的70%～90%[3]。对于大多数PTC患者，预后良好，然而，在某些情况下会出现侵袭性行为，30%～90%的患者表现出临床或隐匿性颈部淋巴结转移[4-5]。lncRNA已被证实参与PTC的细胞迁移和增殖过程[6]。因此，寻找肿瘤与正常组织之间的差异lncRNA有助于进一步了解甲状腺乳头状癌的发病机制，为PTC的术前诊断提供生物标志物和治疗靶点。

基于生物信息学分析，研究显示，lncRNA AC005479.2在甲状腺乳头状癌中表达与癌旁组织差异有统计学意义(P<0.05)，同时GO和KEGG结果显示AC005479.2可能通过脂肪细胞因子信号途径、MAPK信号通路、钙黏蛋白结合参与细胞黏附、对凋亡信号通路的正调控等功能参与甲状腺乳头状癌的病理过程。

近年来研究显示，肥胖患病率的增加与包括乳腺癌、结肠癌和前列腺癌在内的各种人类癌症的发病率、进展和预后较差有关[7]。最近，一些研究显示甲状腺癌的发病率随着肥胖率的显著上升而增加，并且越来越多的证据表明肥胖与甲状腺癌风险增加之间存在关联[8-9]。脂肪因子是脂肪细胞分泌的各种酶、激素、细胞因子、生长因子、蛋白质和其他生物活性物质，包括脂联素、瘦素、抵抗素和白细胞介素等。据报道，体重指数越高的患者脂肪因子浓度越高[10]。由于营养过剩或运动不足所致肥胖时，脂肪因子将出现异常分泌，并参与多种病态的发生。脂肪因子参与以下病理和生理过程，如胰岛素敏感性、炎症和增殖，这些在肿瘤发生和发展过程中均发挥重要作用。因此，脂肪因子可能是将肥胖与甲状腺癌联系起来的目标之一。本研究显示，AC005479.2可能通过脂肪因子参与了甲状腺乳头状癌的发生发展过程。

MAPK通路代表一种普遍存在的信号转导通路，在细胞活力、凋亡和侵袭的调控中发挥重要作用。MAPK信号通路能够被受体酪氨酸激酶、G蛋白偶联受体和部分细胞因子受体激活，将细胞外的信号刺激传递给细胞，引起细胞增殖、凋亡、分化等一系列反应，进而促进肿瘤的发生发展。先前的研究表明，MAPK信号通路在多种癌症中被激活，包括胃癌[11]、肺癌[12]、卵巢癌[13]和肝癌[14]。这一途径已被证明在肿瘤发生发展中发挥关键作用，因此，成为包括甲状腺癌在内的恶性肿瘤的重要治疗靶点[15]。MAPK属于丝氨酸/苏氨酸蛋白激酶，MAPK信号级联是肿瘤细胞增殖、分化、凋亡和耐药性的关键途径。此外，MAPK激活可以磷酸化核转录因子和蛋白激酶，调节相关基因的转录并参与各种生理过程[16]。MAPK信号通路由4个子通路组成：extracellular signal-regulated kinase(ERK)、c-Jun N-terminal kinase(JNK)、p38 和big MAPKinase(BMK)。研究表明，ERK信号通路的激活促进细胞生长、分化、迁移和存活[17]。JNK和p38的激活涉及细胞分化、凋亡和细胞存活的变化[18]。研究显示，MAPK通路中，BRAF的突变及RET/ PTC重排能够促进甲状腺滤泡细胞向乳头状甲状腺癌的转化，这种改变被认为是甲状腺癌发生和进展的早期表现[19]。据Nucera等[20]报道，根据甲状腺癌的遗传学变异，有28%～86%的患者伴随BRAF基因突变。其中BRAF V600E突变占甲状腺乳头状癌的45%[19]。Zhang等[21]研究显示BRAF基因激活的非编码RNA(BRAF-activated non-coding RNA，BANCR)通过激活ERK-MAPK信号通路进而对乳头状甲状腺癌细胞增殖和迁移产生影响。

MAPK还可通过调节上皮-间质转化(epithelial mesenchymal transformation EMT)信号通路发挥作用，上皮-间质转化被认为是癌细胞从原发部位向远处组织或淋巴结转移过程中的一个重要步骤，通过该过程上皮细胞逐渐获得间充质特征并增强迁移和侵袭的特性。甲状腺乳头状癌的多灶性病变显示出淋巴结转移倾向[22]，PTC的淋巴结转移导致预后不良的结果，研究显示EMT参与了PTC细胞的局部转移过程。EMT的分子特征包括将极化上皮表型标志物(如E-cadherin)转换为间充质标志物(如N-cadherin、纤连蛋白等)。转录因子(如Snail、Twist、ZEB1和ZEB2等)的激活触发EMT的发生。据报道，MAPK/ERK的激活会增加Snail、Slug和N-钙黏蛋白的水平，并降低E-钙黏蛋白的表达水平，从而导致EMT[23-25]。EMT通路已被广泛证明在肿瘤进展和转移中起关键作用，可作为甲状腺癌的临床治疗靶点[15,26]。本研究显示AC005479.2可能通过MAPK通路，改变细胞黏附状态，形成EMT环境，促进PTC的淋巴结转移，导致PTC的进展。

综上所述，本研究基于TCGA甲状腺乳头状癌及癌旁组织的基因数据，显示AC005479.2在甲状腺乳头状癌和癌旁组织中的表达差异有统计学方法意义，而且通过GO和KEGG分析其可能通过MAPK通路、脂肪因子通路，调节EMT及凋亡的发生，从而参与甲状腺乳头状癌的发生发展过程，为验证AC005479.2的作用提供理论依据。