生物信息学预测橘皮素治疗结直肠癌核心基因

宋子健，李建伟，胡和智

作为第三大恶性肿瘤的结直肠癌(colorectal cancer, CRC)具有恶性程度高、病程进展迅速、易复发和转移等特点，对人类健康和生命安全构成重大威胁[1]。目前对于CRC的分子机制和核心基因的不完全了解阻碍了对CRC的各项研究。橘皮素是一类天然黄酮类物质，属于黄酮类化合物，广泛存在于芸香科植物川橘果皮、酸橙果皮和柑橘茎叶中，目前已被证实具有抑制细菌和抗肿瘤等药理作用[2]。因此，该研究利用橘皮素治疗结直肠癌的RNA-seq数据预测核心基因，并利用生存分析对其进行预后分析，为CRC的诊断及治疗药物的研制提供新的作用靶点。

1 材料与方法

1.1 数据资料收集使用DMSO溶剂和橘皮素药物处理CRC的HCT116细胞48 h，使用TRIzol提取实验组和溶剂对照组细胞的总RNA，使用逆转录试剂盒将其逆转为cDNA，进行RNA测序。最后，将细胞样本分为3个橘皮素实验组和3个无橘皮素对照组。

1.2 差异基因筛选方法本研究筛选差异表达使用R语言中的程序包对橘皮素实验组和无橘皮素对照组的RNA-seq数据进行基因差异表达分析，差异基因的筛选标准为|log2FC|>1和P<0.05。

1.3 lncRNA关键基因筛选根据基因种类，从筛选出的差异基因集中分别提取lncRNA差异表达基因集、miRNA差异表达基因集和mRNA差异表达基因集。对lncRNA差异表达基因集采用如下方法筛选得到lncRNA关键基因：①提高筛选标准，以|log2FC|>2和P<0.05为新阈值，分别筛选出差异表达更为显著的lncRNA差异表达基因集和miRNA差异表达基因集，分别记为集合A和集合B，通过归并排序算法对两个集合均按照P值由小到大进行排序；②从StarBase数据库[3]中收集与miRNA差异表达基因集B中miRNA存在调控关系的lncRNAs，记为lncRNA基因集C；③将集合A与集合C取交集，得到的lncRNA基因集记为集合D；④利用GEPIA数据库[4]中的临床数据对集合D中的lncRNAs进行生存分析，得到有显著预后价值的lncRNA关键基因集。显著lncRNA能够通过大数据网站(如GEPIA等)预测它们和临床病理参数的关系，以供后续研究使用。

1.4 lncRNA-miRNA-mRNA调控网络构建首先，利用DIANA网站[5]获取lncRNA-miRNA调控关系数据，然后通过miRDB网站[6]获取miRNA-mRNA调控关系数据。利用Cytoscape3.7.2软件[7]构建lncRNA-miRNA-mRNA调控网络，将其中的mRNAs记为mRNA关键基因，以便后续研究。

1.5 基因功能注释和通路富集分析为了更深层次了解橘皮素在治疗CRC中的调控功能，对mRNA关键基因集进行基因本体论(gene ontology,GO)分析和京都基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)通路分析，从而找到差异基因的分子功能、参与的主要生物过程以及它们所属的代谢通路等。DAVID网站[8]是一个生物信息数据库，它整合了生物学数据和分析工具，主要用于GO富集分析和KEGG通路分析。将mRNA关键基因集导入到DAVID数据库中进行分析，相关阈值设置为P<0.05、Kappa Score=0.5、Min Level=5和Max Level=8。

1.6 mRNA核心基因集筛选通过STRING数据库，收集mRNA关键基因集的蛋白质互作网络(PPI)数据，PPI分数设置为0.4，并使用Cytoscape3.7.2软件构建相应的PPI网络[9]。基于MCODE算法和cytoHubba插件中拓扑分析方法Degree、MNC和EPC分别对PPI网络进行分析。MCODE算法中的阈值设为：Node Density Cutoff=0.1、Node Score Cutoff=0.2、K-Core=2、Max.Depth=100;3个拓扑分析方法的阈值设置为：Degree Score≥3、MNC≥2、EPC≥4.8，最终得到mRNA核心基因集。

1.7 统计学处理通过GEPIA数据库对mRNA核心基因集进行在线生存分析，验证其有效性。数据集限定为COAD和READ数据集，时间轴单位设置为月；基因表达差异采用t检验，在CRC中表达量与预后的关系采用Log-rank检验，以P<0.05表示差异有统计学意义。核心mRNA能够通过大数据网站(如GEPIA等)预测它们和临床病理参数的关系，以供后续研究使用。

2 结果

2.1 差异基因筛选RNA测序后共有21 460条基因数据，根据阈值P<0.05，|log2FC|>1对基因数据进行差异表达分析，共得到2 614个差异表达基因，上调差异表达基因1 711个，下调差异表达基因903个。其中，含有197个lncRNA差异表达基因，128个miRNA差异表达基因，1 938个mRNA差异表达基因。差异表达基因的火山图如图1所示，图中红色为上调基因，绿色为下调基因。

图1 RNA-seq数据的差异表达分析火山图

2.2 lncRNA关键基因集以阈值|log2FC|>2且P<0.05为标准对197个lncRNA差异表达基因再次进行筛选，得到116个lncRNA基因，命名为集合A；以阈值|log2FC|>2且P<0.05为标准对miRNA差异表达基因集进行筛选，得到92个miRNA基因，命名为集合B；从StarBase数据库中共收集到65个与集合B存在调控关系lncRNA基因，命名为集合C；对集合A和集合C取交集，得到32个lncRNA基因，命名为集合D；利用GEPIA数据库进行生存分析，共得到5个具有显著预后价值的lncRNA基因(MALAT1、NEAT1、LINC00342、LINC01133、LINC00662)，记为lncRNA关键基因集，生存曲线如图2所示。由图2可知，5个lncRNA关键基因的LogrankP均<0.05，MALAT1(LogrankP=0.022)、NEAT1(LogrankP=0.013)、LINC00342(LogrankP=0.035)、LINC01133(LogrankP=0.036)和LINC00662(LogrankP=0.045),这证明5个lncRNA关键基因均对患者的生存周期有显著影响。

2.3 lncRNA-miRNA-mRNA调控网络针对5个lncRNA关键基因，从DIANA网站获取651条与它们相关的lncRNA-miRNA调控关系数据，从miRDB数据库获取419条与它们相关的miRNA-mRNA调控关系数据。基于以上两组调控关系数据，利用Cytoscape3.7.2软件的Merge功能构建相应的lncRNA-miRNA-mRNA调控网络。lncRNA-miRNA-mRNA网络由1 057个节点和909条边组成，共有117个mRNA关键基因。调控网络如图3所示。

2.4 mRNA关键基因的GO和KEGG分析结果通过DAVID数据库对117个mRNA关键基因进行GO富集分析，在生物过程方面，mRNA关键基因功能主要集中于RNA聚合酶II启动子转录的调控和TOR信号通路等生物过程，表1给出了生物过程中差异性显著的前5个条目。KEGG通路分析结果表明CRC的发生发展与PI3K-Akt信号通路密切相关。表2给出差异性显著的前5条KEGG通路分析通路。通过GO富集分析和KEGG通路分析证实mRNA关键基因集与CRC密切相关。

图2 5个lncRNA关键基因的生存曲线图

表1 差异性显著的前5个GO条目

表2 差异显著的前5个KEGG通路分析条目

2.5 mRNA核心基因集将117个mRNA关键基因上传至STRING数据库，构建对应的PPI网络，将结果导出并保存为tsv格式。利用Cytoscape3.7.2软件将该PPI网络可视化，由48个节点和45条边组成，如图4所示。利用MCODE算法和cytoHubba插件中Degree、MNC及EPC拓扑分析方法对PPI网络进行分析。MCODE算法分析结果如图5所示，拓扑分析方法前10名结果如表3所示。利用Python语言的numpy程序包对MCODE算法的结果基因集和拓扑分析方法的结果基因集取交集，得到6个mRNA核心基因(FOS、GADD45A、CCND2、MYCN、BACH1和MXD1)。这些核心基因在橘皮素治疗CRC中起到重要的调控作用，有成为生物标志物和药物靶点的潜力。

表3 拓扑分析方法前10名结果

图3 lncRNA关键基因对应的lncRNA-miRNA-mRNA调控网络

图4 mRNA关键基因对应的PPI网络

2.6 mRNA核心基因与患者的预后关系利用GEPIA数据库对6个mRNA核心基因进行生存分析，各自的生存曲线图如图6所示。其中，FOS、CCND2和MXD1表达水平对患者的总生存时间有着显著影响(P<0.05)。而GADD45A、MYCN和BACH1对患者的生存率影响差异无统计学意义。

图5 MCODE分析结果图

3 讨论

CRC的发生与外界环境、行为方式、遗传等多种因素密切相关，尽管目前对CRC的研究已取得了较大进步，但是CRC的预后仍然效果不佳。随着测序技术的飞速发展和生物信息技术的不断突破，CRC的分子机制研究成为了当前的一个热点，寻找CRC诊断及预后的生物标志物和药物靶点为CRC的诊疗提供了新的思路。橘皮素已被证实具有抑制细菌和抗肿瘤等药理作用，但对橘皮素治疗CRC的分子机制却不甚了解。因此，本文以橘皮素治疗CRC的RNA-seq数据为研究对象，通过生物信息学分析方法，筛选出117个mRNA关键基因。GO富集分析和KEGG通路分析结果表明，关键基因均富集在与CRC相关的功能和通路上。其中，KEGG分析结果表明PI3K-Akt信号通路与CRC密切相关。研究[10]表明，SPOCK1在CRC细胞系中过表达，沉默SPOCK1可逆转CRC细胞中的EMT过程，显著减弱了迁移/侵袭，抑制体外增殖和体内肿瘤的生长。敲除SPOCK1明显降低了HCT116细胞中p-PI3K和p-Akt的蛋白表达水平。此外，mRNA关键基因还显著富集在乙型肝炎、肺结核、胰腺癌、甲型流感、前列腺癌等多种疾病的相关信号通路，提示橘皮素对多种疾病具有治疗作用，为今后的相关研究提供了新的思路。

图6 6个mRNA核心基因表达与患者预后的生存曲线

通过生存分析，从mRNA关键基因中筛选出3个与CRC预后密切相关的核心基因(FOS、CCND2和MXD1)。其中，FOS和CCND2已被文献[11-12]证实与CRC有密切关系。有研究[13]表明，MXD1参与了乳腺癌癌细胞的增殖和转移过程，在乳腺癌组织中MXD1表达显著下调，并影响了乳腺癌患者的预后。因此，课题组推断MXD1也极有可能与CRC的发生和发展相关。综上所述，3个核心基因有成为生物标志物和药物靶点的可能，对CRC的发病机制及治疗提供了新的思路，也为CRC药物靶点研究提供重要参考。