基于基因芯片整合分析筛选慢性牙周炎相关基因和转录因子

曾小丽 李生娇 单铮男 殷俊豪 姜吉蕊 郑章龙 李家

1.上海牙组织修复与再生工程技术研究中心，同济大学附属口腔医院·同济大学口腔医学院口腔颌面外科，上海200072；2.上海牙组织修复与再生工程技术研究中心，同济大学附属口腔医院·同济大学口腔医学院修复科，上海200072

牙周炎是一种常见的口腔科疾病，慢性牙周炎的持续进展将导致牙齿的松动、脱落乃至口腔咬合功能的完全丧失[1-3]，作为被多种因素混杂影响的慢性炎症性疾病，牙周炎的发生、发展过程中有大量基因及基因产物参与其中，因此，寻找牙周炎发病过程中的关键基因及相关因子是目前研究的一大热点。

随着高通量技术的发展，通过芯片或测序等技术研究疾病与对照组间的基因表达谱的差异成为筛选致病基因的有力工具。基因表达量的变化往往预示着生理过程的改变或病理过程的发生，与此同时，基因编码的产物所涉及的分子功能、生物过程、信号通路则可能与疾病的发生发展相关。由于慢性牙周炎的发生发展是多种基因及其产物的共同作用效果，单个基因或蛋白的表达变化往往不足以完全反映疾病的发展过程。此外，由于样本数目小以及芯片检测平台的差异，过去不同研究者获得结果并不完全一致，甚至可能互相矛盾。因此，本实验基于大样本量、相同检测平台等检索要求，从基因表达综合（gene expression omnibus，GEO） 数据库中筛选并获得符合标准的两个慢性牙周炎相关数据集GSE10334 和GSE16134 并深入分析其中差异表达基因的功能、相关通路以及核心模块、基因，为未来的研究指明方向[4-5]。

1 材料和方法

1.1 基因芯片数据获取

在GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/） 中以 “牙周炎（periodontitis）” 为关键词检索相关芯片，并筛选基于同一检测平台且样本量大于100 的基因表达数据集，获取了2 组符合标准的数据集：GSE10334 和GSE16134。两组数据集均采用美国Affymetrix 公司提供的Human Genome U133 Plus 2.0 Array 基因芯片进行检测，平台号为GPL570。两组芯片数据的样本采集于中重度牙周炎患者炎性区域牙龈组织（试验组） 和健康牙龈组织（对照组），两组使用相同的采样方法，并具有一致的样本纳入标准和临床分型指标[4-5]。其中，GSE10334 的芯片数据来源于90 位患者的247 例样本（183例病变，64例健康）；GSE16134的芯片数据来源于120 位患者的310 例样本（241 例病变，69例健康）[4-5]。

1.2 筛选差异表达基因

使用GEO2R 在线分析工具（http://www.ncbi.nlm.nih.gov/geo/geo2r/） 进行试验组与对照组之间差异表达基因（differentially expressed gene，DEG）的筛选，筛选阈值设定为调整P值（adjustedPvalue，adj.P） ＜0.05，差异倍数（fold change，FC） ＞2。使用R 语言的ggplot2 软件包绘制火山图和热图。

1.3 进行功能富集分析

对两组数据集的差异基因取交集后进行后续分析，以增加结果可信度。使用g: Profiler （https://biit.cs.ut.ee/gprofiler/gost） 分别对上、下调基因进行功能富集分析，功能富集分析包括对基因的基因本体论（gene ontology，GO） 分析和通路分析。GO 分析涵盖了分子功能、生物过程和细胞组分3个下级类目。通路分析中的生物通路信息来源于京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）、Reactome和WIKI共3个数据库。adj.P＜0.05 视为差异具有统计学意义。

1.4 建立蛋白互作网络

为了探索差异表达基因之间的相互作用关系，笔者通过在线工具STRING 数据库（https://stringdb.org/） 建立蛋白相互作用（protein-protein interaction，PPI） 网络，然后使用Cytoscape （www.cytoscape.org/） 软件将该网络进行可视化。

1.5 核心模块与核心基因确立

使用Cytoscape 软件中的MCODE 插件并基于默认阈值（K-core=2，Node Score Cutoff=0.2，Degree Cutoff=2，Maximum Depth=100） 筛选PPI 网络中的重要功能模块。使用cytoHubba 插件进行PPI 网络中关键节点的筛选。cytoHubba 插件使用12 种拓扑分析策略进行重要节点的预测和评估，在6 种以上算法中得分均位于前10 的节点被视作核心基因。为进一步了解PPI网络中重要子网络的功能，笔者对MCODE核心模块中的基因进行了功能富集分析，并绘制通络分析气泡图。

1.6 预测转录因子

使用Cytoscape 软件中的iRegulon 插件进行转录因子的预测。该插件使用标准化富集分数（normalized enrichment score，NES） 来评估预测结果的可信度，NES 值越大，说明可信度越高，笔者将NES＞4.5的转录因子将用于调控网络的建立。

2 结果

2.1 筛选出的差异表达基因

基于adj.P＜0.05 且FC 绝对值＞2 的筛选标准，在GSE10334 中筛选出221 个DEG，其中150 个上调，71 个下调（图1A）；在GSE16134 中筛选出297 个DEG，其中212 个上调， 85 个下调（图1C）；两组DEG 在试验组与对照组的表达量变化及其基因、样本聚类结果见图1B、D。

图1 试验组与对照组间差异基因的筛选Fig 1 Identification of differentially expressed genes between test group and control group

2.2 差异表达基因的功能富集分析结果

在两组数据集中，196 个DEG 表达趋势一致，包括139 个上调基因和57 个下调基因（图2A、D）。GO 分析结果显示：上调基因共富集于10 项分子功能、99 项生物过程以及20 项细胞组分，下调基因共富集于1项分子功能、11项生物过程以及7项细胞组分。图2B、E展示了上、下调DEG adj.P前10 的GO 条目。在分子功能方面，上调基因参与趋化因子受体结合、免疫球蛋白结合，而下调基因参与皮肤表皮结构形成。在生物过程方面，上调基因涉及免疫应答、刺激应答、白细胞活化等，而下调基因涉及皮肤及表皮的发育和角化。在细胞组分方面，上调基因主要位于胞外间隙、细胞表面、胞外体、质膜外侧，而下调基因主要位于超分子纤维、超分子聚合物、聚合细胞骨架纤维、透明角化颗粒等处。

通路分析结果显示：上调基因富集于26 条信号通路，其中包括12 条KEGG 通路、12 条Reactome 通路和2 条WIKI 通路；而下调基因富集于2条Reactome 通路。图2C、F 分别展示了上调基因adj.P前15 的通路条目和下调基因所富集的2 条通路条目。上调基因参与细胞因子-细胞因子受体相互作用、白细胞经内皮细胞迁移、趋化因子信号通路以及人补体系统、造血干细胞分化等免疫反应相关通路，而下调基因涉及角质化包膜形成、角化等2条通路。

图2 差异基因功能富集分析结果Fig 2 Functional enrichment analysis of differentially expressed genes

2.3 PPI网络建立及核心模块、核心基因筛选结果

由交集基因构成的PPI 网络包含132 个节点和519 个节点作用关系（图3）。图4A、B 显示：PPI网络的第一模块涉及病毒蛋白与细胞因子及细胞因子受体的相互作用、细胞因子-细胞因子-受体相互作用和趋化因子信号通路等通路。图4C、D 显示：第二模块参与B细胞受体信号通路。图4E、F显示：第三模块与白细胞介素介导的信号通路、JAK-STAT信号通路等通路有关。

图3 构建PPI网络并筛选核心基因Fig 3 Construction of PPI network and identification of hub genes

图4 PPI网络核心模块筛选与通路分析结果Fig 4 Identification and pathway analysis of key modules in PPI network

cytoHubba 插件的分析结果显示：满足筛选标准的基因为CXC 基序趋化因子配体（C-X-C motif chemokine ligand，CXCL） 8、CXCL1、CXC 基序趋化因子受体（C-X-C motif chemokine receptor，CXCR） 4、选择素（selectin，SEL） L， CD19 分子（CD19 molecule，CD19） 和IKAROS 家族锌指（IKAROS family zinc finger，IKZF） 1。在PPI 网络中，CXCL1、CXCL8、CXCR4 和SELL 位于第一模块，IKZF1位于第三模块。

2.4 转录因子预测结果与调控网络的建立

依据iRegulon预测结果，NES＞4.5的转录因子为干扰素调节因子（interferon regulatory factor，IRF） 4、组蛋白乙酰转移酶（E1A binding protein，EP） P300、易洛魁同源盒蛋白（iroquois homeobox，IRX） 4 和TATA 序列结合蛋白（TATA-box binding protein，TBP），分别调控60、18、24 和36个DEG（图5），其中IRF4属于上调DEG且位于第二核心模块，参与B细胞受体信号通路（图4 C、D）。

图5 转录因子及其调控网络Fig 5 Transcription factors and their regulatory networks

3 讨论

牙周炎不仅是成年人失牙的主要原因，也是多种疾病如糖尿病、动脉粥样硬化、慢性阻塞性肺疾病等系统性疾病的潜在风险因素[3,6]。鉴定与牙周炎进展相关的相关分子及信号通路，有利于了解牙周炎的发病机制，为牙周炎的早期诊治及伴发疾病的预防提供理论基础。本研究基于GEO数据库公开的牙周炎基因表达谱，通过系统性生物信息学分析，为牙周炎的机制研究提供理论依据和研究方向。由于样本量小会在一定程度上降低试验结果的可信度，因此本研究筛取大样本数据集GSE10334 和GSE16134 （分别含有247 例和310 例样本） 进行生物信息学分析。联合分析两组芯片数据集，在两组中均差异表达的基因共196个，其中139 个为上调基因，57 个为下调基因。功能富集分析结果显示，上调基因参与免疫系统，病毒蛋白与细胞因子及细胞因子受体的相互作用，细胞因子-细胞因子受体相互作用，白细胞跨内膜迁移和趋化因子受体结合趋化因子等免疫反应相关途径，由此可见免疫应答相关途径在牙周炎发病过程中发挥重要作用。而下调基因涉及角质化包膜形成、角化等途径，说明上皮屏障形成能力降低，提示牙周炎的发生与牙龈上皮完整性丧失相关。由差异基因构建的PPI 网络中，MCODE 得分前三的模块被视作核心模块，分别富集于趋化因子信号通路、B细胞受体信号通路以及白细胞介素介导的信号通路等途径，由此推测，这些分子事件在牙周炎发病过程中活跃。CXCL8、CXCL1、CXCR4、SELL、CD19和IKZF1为核心基因，位于PPI网络的中心位置，因此最有可能成为牙周炎的治疗靶点。

CXCL8 和CXCL1 属于趋化因子CXC 亚家族，通过募集白细胞来增强或延长炎症反应。多项研究[7-9]发现，相较于健康牙龈，炎症牙龈组织切片中CXCL8 和CXCL1 的mRNA 水平和蛋白水平均显著升高。这些结果提示：炎性环境下牙周细胞可能通过分泌更多的趋化因子以对抗牙周病原体的感染，而这种过度的免疫应答正是牙周组织破坏的主要原因。此外Souto 等[1]研究发现，CXCL8水平与牙周炎患者的探诊出血、探诊深度和临床附着丧失等临床检测指标正相关，其结果提示CXCL8 表达量与牙周炎严重程度相关，有望成为牙周炎进展标志物。

CXCR4 是CXC 趋化因子受体家族成员，Hajishengallis 等[10]发现，牙龈卟啉单胞菌可通过结合巨噬细胞表面CXCR4，逃避宿主的免疫清除。口服CXCR4 的拮抗剂AMD3100，可抑制牙龈卟啉单胞菌的定殖，进而预防或减轻小鼠实验性牙周炎及其引起的牙槽骨吸收[11]。因此推测，CXCR4是牙龈卟啉单胞菌在牙周炎致病作用的一个关键位点。

SELL 编码淋巴细胞归巢受体L-选择素，该因子属于细胞表面黏附分子，调控白细胞在血管内的黏附与迁移。L-选择素主要表达于中性粒细胞，在炎性疾病中，L-选择素结合其相应配体，调控中性粒细胞活化与聚集，并使其从血管中迁移至炎症区域。中性粒细胞是固有免疫中的重要组成细胞，在免疫防御中起先驱作用，但近年来的研究[12]发现，过度活化的中性粒细胞可通过产生大量活性氧和基质降解酶造成牙周组织的破坏。因此，SELL 可能通过促进中性粒细胞迁移和活化参与牙周炎病理发生。

CD19 是特异性表达于B 细胞表面的免疫球蛋白基因超家族成员[10-11,13-14]，为B细胞抗原受体复合物的共受体，可降低激活下游信号通路和触发B细胞对抗原的反应的阈值[14]。因此，CD19 在试验组的表达量升高，提示在炎性牙龈中B 细胞数目和免疫应答活性增加。有研究[15]表明，B 细胞是牙周炎骨缺损区分泌RANKL的主要细胞之一。使用抗B 细胞药物如利妥昔单抗则可以缓解牙周炎的多种症状[2]。因此，以CD19 为靶点拮抗B 细胞的药物研发，可能成为药物治疗牙周炎的一个研究方向。

IKZF1 编码淋巴系特异性转录因子IKAROS，IKAROS 是造血细胞分化的关键调节因子，调控T淋巴细胞和B 淋巴细胞抗原受体介导的增殖反应以及相关效应细胞的功能[16]，T 细胞与B 细胞参与牙周炎结缔组织区域的炎症浸润及骨丧失[2,15]。

IRF4 为iRegulon 预测结果中可信度最高的转录因子， 调控包括核心基因CXCR4、 CD19 和IKZF1在内的60个DEG的转录，IRF4 mRNA与蛋白质在炎性牙龈中的表达已在过去的文献[7]报道中得到支持。虽然IRF4 在牙周炎发病过程中的具体作用尚不明确，但IRF4 参与系统性红斑狼疮、类风湿性关节炎等自身免疫性疾病的发生[17]，提示其与机体免疫亢进相关。

在本研究中，为在多样本高通量数据中快速高效找到显著上调或下调的基因，笔者使用了严格的阈值筛选方法，但值得注意的是，这种方法可能会遗漏部分具有重要生物学意义但表达量变化并不显著的基因。例如，白细胞介素-1β 和肿瘤坏死因子-α 作为牙周炎疾病中重要的炎性因子和骨吸收介质，可激活下游核因子κB 信号通路参与牙周炎相关骨丢失[18]。然而由于不满足FC＞2 的阈值标准，这些牙周炎相关基因在本实验中不作为“差异表达基因” 纳入后续分析，其在牙周炎发病过程中的作用也未得到体现。因此本研究尚存在一定局限性，相关分析结果需要通过后续的分子生物学实验和动物实验进一步验证。

综上所述，本研究基于大样本芯片数据集研究，通过生物信息学分析筛选、鉴定与牙周炎发病相关的易感基因、转录因子及相关信号通路。共发现196 个牙周炎相关差异基因，包括6 个核心基因（CXCL8、CXCL1、CXCR4、SELL、CD19和IKZF1） 和1 个关键转录因子（IRF4），这些分子通过细胞因子-细胞因子受体相互作用、趋化因子受体结合趋化因子等多条免疫相关途径参与牙周炎的病理发生。本研究首次报道了IKZF1、IRF4与中重度牙周炎患者牙龈组织炎性病变相关的研究。病毒蛋白与细胞因子及细胞因子受体的相互作用、白细胞跨内膜迁移、B细胞受体信号通路等途径是本研究新发现的与牙周炎发病相关的信号通路。针对以上分子及通路展开深入研究，将有助于牙周炎发病机制的阐明和治疗靶标的确立。

利益冲突声明：作者声明本文无利益冲突。