溃疡性结肠炎及其恶性并发症的生物信息学分析和潜在治疗药物筛选

夏守兵，许春杰，蒋春晖，顾 磊，孙隆慈，徐 庆

上海交通大学医学院附属仁济医院胃肠外科，上海 200127

溃疡性结肠炎（ulcerative colitis，UC），又称特发性溃疡性结肠炎，是一种起始于直肠的非特异性慢性肠道炎症性疾病，目前发病原因尚不明确。UC可发病于任何年龄，但多见于青春末期和成年早期。临床症状主要以持续或反复发作的腹泻、黏液脓血便为主，可伴有不同程度的全身症状。迁延不愈的患者最后可发展为炎症相关恶性肿瘤[1-2]。目前普遍认为，UC的发病是环境、患者自身的易感基因、精神心理因素及肠道微生物等多种因素共同叠加的结果[3-4]，但针对UC及其恶性并发症的致病相关基因尚未得到充分研究。本研究旨在通过分析比较UC样本与正常样本，以及UC样本与UC伴瘤变样本的基因表达谱数据，构建关于差异表达基因编码蛋白质相互作用（proteinprotein interaction，PPI）网络，对UC及其恶性并发症的致病相关基因进行系统的生物信息学分析，从而为进一步探究UC及其恶性并发症的具体分子机制提供新的研究思路。

1 材料与方法

1.1 基因表达数据集获取及数据预处理

从美国国立生物技术信息中心（National Center for Biotechnology Information，NCBI）的基因表达数据库（Gene Expression Omnibus，GEO）（https://www.ncbi.nlm.nih.gov/gds） 中 下 载 编 号 为 GSE13367、GSE9452、GSE36807、GSE37283的基因表达数据集及样本临床信息，其中GSE13367数据集包含16例UC样本及20例对照样本，GSE9452数据集包含8例UC样本及18例对照样本，GSE36807数据集包含15例UC样本及7例对照样本，GSE37283数据集包含5例健康人群对照样本、4例静止期UC样本以及11例UC伴瘤变样本。GSE13367、GSE9452、GSE36807基因表达数据集都是基于GPL570注释平台中Affymetrix（昂飞）Human Genome U133 Plus 2.0 Array [HG-U133_Plus_2]的人类全基因转录本生物信息，GSE37283基因表达数据集是基于GPL13158注释平台中Affymetrix（昂飞）HT HG-U133+ PM Array Plate的人类全基因转录本生物信息，故以GSE13367、GSE9452、GSE36807数据集作为UC组，GSE37283为UC伴腺瘤组（UCN组）。首先依据UC组的K均值聚类算法结果对该组的3个数据集样本进行无监督聚类分析，剔除临床表型与聚类分析结果不符的样本。

1.2 差异表达基因筛选

基因表达谱数据经注释及聚类分析等预处理之后，利用R语言limma函数包分别对UC组与UCN组进行差异基因分析，筛选出P＜0.05且|log2FC|＞1.5的基因作为显著差异表达基因，并通过基于R语言的VennDiagram函数包对UC组中3个数据集进行单独分析以及联合分析得到差异基因，绘制韦恩图以分析差异基因的交并集，再利用ggplot2包绘制差异基因的火山图。

1.3 基于GSEA的KEGG富集分析

基因集富集分析（Gene Set Enrichment Analysis，GSEA）可用来评估一个预先定义的基因集中的基因在以表型相关度排序的基因列表中的分布趋势，从而判断其对表型的贡献。为进一步探究UC致病相关基因的生物学功能，将从MSigDB数据库（http://software.broadinstitute.org/gsea/msigdb/）中获得的KEGG功能基因集作为参照基因集[5]，采用R语言GSEA函数包分析以上2组差异基因数据 在 KEGG（Kyoto Encyclopedia of Genes and Genomes）信号通路中的富集情况。通过定义|NES|＞1，NOM P-value＜0.05，FDR Q-value＜0.25的通路下的基因集合是有意义的，筛选得到2组差异基因在KEGG信号通路中的富集情况及其所包含的具体基因集。

1.4 PPI网络构建及模块分析

STRING数据库（https://string-db.org/）主要通过实验数据、文献报道摘要和来自其他数据库的综合数据实现预测PPI功能。在本研究中，根据显著差异基因符合P＜0.05且|log2FC|＞1.5的筛选条件，使用STRING数据库分析显著差异表达基因的PPI网络，并利用Cytoscape 3.6.1软件的CytoHubba插件[6]提供的12种算法对PPI网络中的各个节点基因进行评分分析以识别网络中的核心基因，从而可以较为准确地找出在UC及其恶性并发症病程中起核心作用的致病相关基因。

1.5 差异基因与免疫细胞浸润的相关性分析

TIMER数据库（https://cistrome.shinyapps.io/timer/）是系统分析不同类型癌症组织中的免疫细胞浸润情况的综合数据库。通过此数据库可以查询特定基因与不同肿瘤组织中6种免疫细胞（B细胞、CD4+T细胞、CD8+T细胞、中性粒细胞、巨噬细胞和树突状细胞）的浸润关系、基因表达量与肿瘤患者预后关系以及基因在多种肿瘤和正常组织中的差异表达情况。

1.6 基于LINCS L1000数据库的药物筛选分析

LINCS L1000数据库（https://clue.io/）是由博德研究所（Broad Institute）开发的基于小分子药物、基因表现与疾病相互关联的药物筛选数据库。目前该数据库已收集7 000余个经小分子药物处理后的人类细胞基因表达谱数据，广泛应用于探索药物作用的机制以及研发新的治疗药物。通过利用该数据库查询特异性针对PPI网络核心基因的小分子药物，以发现具有潜在治疗效果的药物。

1.7 统计学分析

采用Graphpad prism 7.0统计学软件对组间的基因表达量数据进行单因素t检验。通过Pearson相关系数分析评估基因表达量之间的相关性。P＜0.05表示差异具有统计学意义。

2 结果

2.1 差异表达基因的筛选

剔除聚类分析结果中与临床表型不符的样本后（图1），通过limma函数包分析UC组与UCN组的差异表达基因。以P＜0.05且|log2FC|＞1.5为标准，针对23例UC样本与33例正常样本的差异基因分析共发现显著差异基因86个，其中43个为上调基因，43个为下调基因；而针对4例静止期UC样本与11例UC伴瘤变样本的差异基因分析共得到253个显著差异基因，其中160个为上调基因，93个为下调基因。从UC组的差异基因韦恩图中可以发现该组3个数据集联合分析得到的差异基因均处在与各数据集单独分析结果的交集中，这表明多数据集联合分析可以有效提高分析结果的可信度（图2A）。利用ggplot2包绘制2组差异基因中上调基因和下调基因的火山图，能够直观地显示2组差异基因的表达情况（图2B、2C）。

图1 剔除与临床表型不符的样本前后的聚类分析结果Fig 1 Results of cluster analysis before and after excluding the samples inconsistent with clinical phenotype

图2 UC组差异基因分析得到的韦恩图及火山图及UCN组差异基因分析得到的火山图Fig 2 Venn diagram and volcano plot obtained from differential genetic analysis of UC group and volcanic map obtained by differential gene analysis of UCN group

2.2 基于GSEA的KEGG富集分析结果

通过GSEA软件对2组基因全转录本进行KEGG富集分析，以 |NES|＞1，NOMP-value＜0.05 且 FDRQ-value＜0.25为筛选标准，发现UC组差异基因的KEGG富集分析结果中显著上调的通路为IL-17信号通路、Toll样受体信号通路、NOD样受体信号通路等多条信号通路，显著下调的通路为柠檬酸循环（TCA循环）、矿物质吸收、丁酸代谢、氧化磷酸化、脂肪酸降解；UCN组的KEGG富集分析结果中显著上调的通路为IL-17信号通路、TNF信号通路、Toll样受体信号通路等，显著下调的通路为TCA循环、类固醇激素生物合成、矿物质吸收、胆汁分泌、近端小管碳酸氢盐重吸收，富集分析得到的具体通路如图3A、3B所示。此外，在2组KEGG富集分析的结果中都发现了炎症性肠病（inflammatory bowel disease，IBD）信号通路和IL-17信号通路，利用R语言enrichplot函数包绘制的这2个信号通路的GSEA富集分析曲线图如图3D所示。基于基因表达量的变化对IBD信号通路及IL-17信号通路进行KEGG信号通路注释，从IL-17信号通路的KEGG通路注释图中可以发现CXCL8（C-X-C motif chemokine ligand 8）参与了对中性粒细胞的招募作用（图3C、4B）。

图3 2组差异基因KEGG富集分析及GSEA分析得到的部分结果Fig 3 Partial results of KEGG enrichment analysis and GSEA analysis of two groups of differential genes

2.3 差异表达基因的PPI网络模块分析

基于STRING数据库的PPI分析得到了UC组显著差异基因的PPI网络（图4A）。基于Cytoscape的CytoHubba插件提供的12种拓扑分析算法对PPI网络各个节点基因进行评分排序，得到网络中的核心节点基因（表1）。评分分析发现CXCL8在10种算法评分排序中均位列第1。在UC组及UCN组的差异基因分析结果中均发现了上调基因CXCL8（图2B、2C），这表明CXCL8可能是在UC及其恶性并发症的发生发展过程中起核心作用的基因。

2.4 TIMER数据库分析以及基因相关性分析

TIMER数据库的分析显示结直肠腺癌组织中CXCL8主要对中性粒细胞的浸润起到趋化作用（图4C），相比于在正常组织中的表达量，CXCL8在多种肿瘤组织中呈高表达状态（图5A）。基于对患者的预后数据分析发现CXCL8呈高表达状态的结直肠腺癌患者预后较差（图5B）。基于Pearson相关性分析发现CXCL8的表达量与中性粒细胞的2种表面抗原以及程序性死亡配体1（programmed deathligand 1，PD-L1）的表达量呈正相关（图4D～4G、5C、5D）。通过对UCN组的分析，发现CD206及PD-L1的表达量随着UC病程的发展逐渐上调（图5E、5F）。

2.5 小分子治疗药物筛选结果

基于LINCS L1000数据库的查询结果，共筛选出3种小分子药物，分别为紫铆因（butein）、左旋西替利嗪（levocetirizine）、伪麻黄素（pseudoephedrine）（图5G）。

图4 CXCL8参与对中性粒细胞招募作用的多种分析结果Fig 4 Multiple analysis results of CXCL8 involved in the recruitment of neutrophils

表1 利用CytoHubba提供的12种算法对UC组PPI网络中的节点基因评分分析得到的排名前10位的节点基因Tab 1 Top 10 node genes obtained by using 12 algorithms of CytoHubba for scoring analysis of the node genes in PPI networks of UC group

图5 UC患者体内CXCL8与PD-L1表达量的相关性分析以及3种小分子药物对CXCL8抑制效果的分析Fig 5 Correlation analysis of the expression of CXCL8 and PD-L1 in ulcerative colitis patients and analysis of the inhibitory effect of three small molecule drugs on CXCL8

3 讨论

本研究通过分析正常人群与UC患者以及UC患者与UC伴瘤变患者的2组差异基因，发现显著上调的差异基因CXCL8，而后对于通过基于Cytoscape的网络核心基因分析插件CytoHubba深入挖掘筛选到的显著差异基因，CytoHubba的评分分析亦发现CXCL8在10种核心基因算法评分结果中均排名第1。基于以上分析结果，我们推测CXCL8可能在UC及其恶性并发症的病程发展中发挥了重要作用。

CXCL8，又称嗜中性粒细胞因子，属于CXC趋化因子家族，主要由外周血中的淋巴细胞、单核细胞以及中性粒细胞等多种细胞产生，是炎症反应的重要介质，常与局部组织的炎症细胞浸润相一致[7]。研究[8]发现慢性炎症组织中的多种免疫细胞可持续性释放CXCL8，CXCL8及其配体CXCR1/2主要负责中性粒细胞的激活过程和诱导中性粒细胞向炎症部位迁移和着陆。本研究中2组差异基因的KEGG功能富集分析的结果也显示了在IL-17信号通路中CXCL8等趋化因子对中性粒细胞有招募作用。CXCL8在肿瘤微环境中对中性粒细胞的招募功能已得到van den Steen等[9]的证实，该研究结果发现CXCL8可以刺激中性粒细胞释放基质金属蛋白酶9（matrix metalloproteinase-9，MMP-9），而MMP-9可以强化CXCL8的趋化能力，使得更多的中性粒细胞聚集在肿瘤组织周围。

慢性炎症长期以来被认为是促发癌症，尤其是结肠癌的一个重要风险因子，超过20%的UC患者会在确诊后的数十年内转化成结肠癌患者[10]。而结肠炎相关性癌变进程缓慢，干预治疗的效果非常差且死亡率高[11]。但长久以来慢性炎症如何抵御免疫反应并最终导致癌症发生与发展的确切机制尚不明确。作为人体固有免疫应答系统的重要组成成分，中性粒细胞可通过多种方式在结肠炎相关性结肠癌的形成过程中发挥促瘤形成作用，例如促瘤型中性粒细胞（N2型）分泌的趋化因子CXCL8与肿瘤细胞表面的趋化因子受体结合后，既能招募更多的中性粒细胞进入肿瘤微环境，又能促进肿瘤血管生成[12-14]。而近年来针对中性粒细胞的深入研究[15-16]发现，肿瘤微环境中N2型中性粒细胞还可通过PD-1/PD-L1信号转导途径抑制T细胞的正常功能，从而在多种癌症中实现促瘤作用。

PD-L1是程序性死亡受体1（programmed death-1，PD-1）的主要配体，研究[17-18]显示两者在结合后可通过“耗竭”肿瘤效应T细胞的方式起到对免疫应答的负性调控作用，并使肿瘤细胞发生免疫逃逸现象。Nakazawa等[19]通过研究UC患者的手术切除标本，发现患者结肠上皮细胞PD-L1表达水平明显高于健康对照者。Lin等[20]通过研究胃癌肿瘤微环境发现，CXCL8能够通过诱导肿瘤相关巨噬细胞表达PD-L1从而抑制CD8+T细胞的功能活化，导致肿瘤细胞逃过免疫检查。Wang等[15]发现胃癌肿瘤微环境中激活的中性粒细胞同样可以通过PD-1/PD-L1通路抑制T细胞的功能，从而抑制其对肿瘤细胞的杀伤作用。以上这些研究结果表明肿瘤微环境对肿瘤相关免疫细胞的“驯化”作用有助于实现PD-L1在这些细胞中的表达上调。

本研究中，我们通过差异基因分析及CytoHubba核心基因算法分析的结果发现显著差异基因PPI网络的核心基因为CXCL8；通过差异基因的KEGG富集分析发现显著上调的CXCL8基因参与了IL17信号通路对中性粒细胞的招募作用；通过基因表达量之间的相关性分析则发现CXCL8与中性粒细胞表面特异性抗原以及PD-L1的表达量呈正相关。以上分析结果在一定程度上验证了TIMER数据库及前述文献报道中CXCL8可通过招募中性粒细胞的方式以实现对PD-L1表达量调控这一结论。同时，有研究[21]发现调控PD-L1表达量的主要信号通路为JAK-STAT信号通路。此外，还有研究[22]发现N2型中性粒细胞表面特异性抗原CD206的表达量也会随着病程的发展而逐渐升高。我们由此推测肠道炎症组织中的中性粒细胞可能会在慢性炎症组织中逐渐被“驯化”而演变发展为激活/抑制的“双重”表型，并在CXCL8的作用下通过JAK-STAT信号通路上调PD-L1的表达量以实现负性调控T细胞功能，从而促使肠道炎症恶性并发症的发生。

此外，通过LINCS L1000数据库查询特异性针对CXCL8的小分子药物，得到紫铆因、左旋西替利嗪、伪麻黄素3种小分子药物。查阅文献[23]发现紫铆因具有抗氧化、抗炎症、诱导肿瘤细胞凋亡等多种药理作用。紫铆因可以显著降低脂多糖（lipopolysaccharide，LPS）刺激小胶质细胞后细胞内多种炎症因子mRNA的表达上调，这说明紫铆因可以抑制LPS诱导的小胶质细胞的活化，从而明显减少炎症因子的产生[24]。而左旋西替利嗪、伪麻黄素作为在临床上应用已久的药物也都可以通过不同的方式达到抑制炎症的效果，其中第3代抗过敏药物的代表左旋西替利嗪可以直接发挥抑制炎症细胞聚集和浸润的功能，而伪麻黄碱则可通过发挥拟肾上腺素作用以减少中性粒细胞在呼吸道细胞内的聚集、黏附。以上3种小分子药物的药理分析表明，抑制肠道黏膜组织内中性粒细胞的长期异常聚集可能是阻止UC及其恶性并发症发生和发展的关键因素。

总而言之，本研究利用生物信息学方法建立的筛选分析流程为研究UC提供了一个新的思路：基于多种数据库对UC的致病相关基因进行全面分析，并对表达量发生变化的基因进行KEGG信号通路分析，再通过PPI网络分析对处于网络中心节点的基因进行深入挖掘，找出与疾病相关的核心基因，从而为研究提供新的突破点。通过综合利用以上方法，本研究发现了一些可能在UC发展过程中发挥核心作用的基因（如CXCL8）。这为UC及恶性并发症的研究提供了新的思路，为理解UC发生过程中的复杂分子机制提供了新的视角，也为进一步推动临床药物研发提供了依据。然而，考虑到本研究尚有一定的局限性，未来还需要进一步开展更为深入的实验来验证本研究的结果。