基于基因测序和GEO数据挖掘探讨结核性溃疡的发病机制

高贝贝，黄子慧，陈思琪，曹丽敏，翁嘉晨，郭丹丹，陈 悦

(南京中医药大学附属南京市中西医结合医院，江苏 南京 210014)

结核病是在世界范围内流行的传染病，是单一传染源导致死亡的主要原因之一，仅次于新型冠状病毒肺炎(COVID -19)。根据世界卫生组织(WHO)[1]估计，我国2020年新发结核病患者数为84.2万，发病率为8.5%，居结核病高负担国家的第二位；同时，受COVID -19的影响，基本结核病服务中断，全球因结核病死亡的患者约增加10万例。近年来，在COVID -19肆虐、结核分枝杆菌变异以及耐药菌种增多等多种因素影响下，结核病发病率上升，治疗难度加大，结核病的防治显得尤为重要。

结核性溃疡(tuberculous ulcer，TU)，又称皮肤结核，是指结核分枝杆菌、牛分枝杆菌以及卡介苗侵犯机体局部组织，引起病灶周围软组织、皮下及皮肤坏死，最终液化破溃形成的创面，是临床上最常见的一类特异感染性疾病。临床主要表现为创周皮色暗红，肉芽苍白水肿，脓液稀薄，形成窦道或慢性溃疡[2]。虽破溃口一般较小，但深层次多有潜行空腔，迁延难愈。目前TU的病因及发病机制不清，早期诊断困难，缺乏敏感和特异性高的诊断标志物，常与其他慢性创面混淆，发病周期长，缠绵难愈，严重影响患者日常生活，给TU的防治带来极大的负担。迄今为止，国内外文献对结核病的研究以肺结核为主，TU属于肺外结核，其发病机制研究较少，以个案报道较为多见。因此，寻找TU病变发生的机制，找到精确的治疗靶点，获得良好预后，从根本上解决TU的问题已成为当务之急。

高通量测序技术迅速发展，生物信息学分析逐步成为科研领域中重要的研究方法。该方法在识别可靠的功能差异表达基因(differentially expressed genes,DEGs)及筛选疾病候选生物标志物方面显示出了巨大的优势。本文通过对基因表达综合数据库及课题组前期基因测序数据集进行筛选，应用多种生物信息技术进行基因分析，采用临床试验进行验证，以期发现与TU高度相关的特异性指标，为TU的靶向研究提供依据。

1 材料与方法

1.1 数据采集 选取课题组前期基因测序的相关数据集[3]，纳入3例TU组织标本和3例TU旁正常组织标本，所有标本均经过病理切片鉴定。研究从公共数据库GEO(https://www.ncbi.nlm.nih.gov/geo/)下载肺外结核(extrapulmonary tuberculosis，EPTB)相关数据集GSE83456，其平台为GPL10558(Illumina HumanHT-12 V4.0 Expression BeadChip)。其中纳入GSE83456数据集47例EPTB血标本和61例健康人体对照(healthy controls，HC)的血标本。GEO芯片样本数据均下载以供进一步分析。

1.2 差异基因的筛选和京都基因与基因组百科全书(KEGG)/基因本体论(GO)富集分析 用R语言(版本：4.1.2)对基因测序及GSE83456的数据列进行处理。采用LIMMA软件包比较试验组和正常对照组之间的差异，识别DEGs，对于具有重复基因符号的探针，使用平均值作为唯一表达值。两组设定P值<0.05和|logfc|>1为筛选条件，筛选获得DEGs。用维恩图将基因测序的DEGs和GSE83456的DEGs相交，筛选与该病相关的靶基因。通过使用R语言中的‘clusterProfiler’包对交集基因进行富集分析，然后使用“ggplot2”包进行可视化。设定P值<0.05为筛选条件，分别筛选出明显富集的GO分析，以及KEGG相关富集通路分析。

1.3 PPI网络的构建 使用检索工具检索相互作用基因(STRING)在线数据库(http://string-db.org；version11.5)构建蛋白质-蛋白质相互作用(protein-protein interaction,PPI)网络。PPI可以揭示疾病发生发展的机制。在此项研究中，使用STRING数据库构建相交基因的PPI网络，得分>0.4的交互作用被认为具有统计学意义。Cytoscape(版本：3.9.0)是一个开源软件平台，用于可视化复杂网络并将其与任何类型的属性数据集成。CytoHubba是一个用于发现PPI网络中枢基因的Cytoscape插件。使用CytoHubba插件中的最大相关准则(MCC)对相交基因进行MCC评分，并将评分按照从高到低的顺序排列，顺序越高，基因的红色越深。

1.4 基因组富集分析(GSEA) GSEA是一种计算方法，用于确定一组先验定义的基因是否在两种生物状态(如表型)之间存在差异。使用GSEA(4.2.1)确定GBP1参与淋巴结结核病潜在的生物学机制。通过进行1 000次基因组排列获得标准化的富集分数(NES)。在当前的研究中探索特征[h.all.v7.5.1.symbols.gmt(Hallmarks)]基因组。设定P值<0.05，错误发现率(FDR)<0.05为筛选条件，以量化统计上的显著富集。

1.5 临床试验验证

1.5.1 试验材料及仪器 包括TRIzolTM LS Reagent(美国Thermo Fisher Scientific公司)，DEPC水，500 mL BR 0.1%(上海源叶生物科技有限公司),PrimScript RT Master Mix(日本TaKaRa公司)，ChamQ SYBR qPCR Master Mix(南京Vazyme公司),引物(上海生工生物公司),高速低温离心机(德国Hettich公司),IMS-25全自动雪花制冰机(常熟雪科公司)，超微量分光光度计(美国Thermo Fisher Scientific公司)，T100梯度PCR仪(购自美国BIORAD公司)，全自动荧光PCR分析仪(瑞士Roche公司)。

1.5.2 纳入及排除标准 纳入标准：(1)符合TU中西医诊断标准的患者；(2)溃疡面积4～25 cm2；(3)患者未接受抗结核及免疫学治疗；(4)患者一般临床资料完整且同意参加本次研究并签署知情同意书。排除标准：(1)合并有造血系统、肝、肾及心脑血管等严重原发性疾病及继发疾病者；(2)免疫功能受损者，如人类免疫缺陷病毒(HIV)感染、白血病、淋巴瘤或全身性恶性疾病等；(3)正在接受免疫抑制治疗的患者，如使用糖皮质激素、烷化剂、抗代谢药物以及放射治疗等；(4)活动性肺结核患者[4]。

1.5.3 研究对象 选择2021年7月—2022年3月在南京中医药大学附属南京市中西医结合医院瘰疬科就诊的60例TU患者为研究对象，所有患者均经过病理诊断确诊，其中男性23例，女性37例，年龄23～72岁，平均年龄(39.2±8.9)岁。同期选取本院健康体检者60例作为对照组，其中男性26例，女性34例，年龄25～67岁，平均年龄(41.7±11.4)岁。受试者自愿参加此项研究并已签署同意书，本研究已获得南京中医药大学附属南京市中西医结合医院伦理科的批准，符合世界医学会赫尔辛基宣言。

1.5.4 血标本采集 采集所有HC及TU患者治疗前和治疗2周后清晨空腹静脉血3 mL，置于离心管中，4℃ 3 000×g，离心10 min，之后转移上清至新离心管中，保存在-80℃冰箱备用。

1.5.5 RNA提取及逆转录聚合酶链反应(qRT-PCR)检测 按照试剂说明书对血清进行提取，并逆转录为cDNA，随后将cDNA产物进行10倍稀释，并使用qPCR试剂盒进行实时荧光定量qRT-PCR分析。每个样品设1个复孔，重复三次试验。qRT-PCR对所有受试者TU及HC血清中GBP1的表达水平进行检测。qRT-PCR反应条件为：95℃预变性30 s，95℃变性10 s，60℃退火与延伸30 s，共计40个循环。采用2-△△Ct对检测数据进行相对定量分析，通过内参基因GAPDH的Ct值对靶基因GBP1的相对表达水平进行归一化处理，并利用GraphPad Prism(8.0.1)对试验结果进行处理。引物序列见表1。

表1 qRT-PCR引物序列

1.6 统计分析 统计分析均由R语言(版本：4.1.2)软件及GraphPad Prism(8.0.1)软件进行。t检验用于检验两组之间的差异，P值为双侧，P≤0.05认为差异具有统计学意义。

2 结果

2.1 EPTB中DEGs的鉴定 从TU患者组织的基因表达序列中筛选出2 490个差异基因，其中1 276个基因上调，1 214个基因下调，绘制相应的火山图(见图1A)和热图(见图1B)。此外，从EPTB患者血液基因表达序列中筛选出79个差异基因，其中69个基因上调，10个下调,绘制相应的火山图(见图1C)和热图(见图1D)。设置P<0.05和|logFC|>1为选择标准。基因测序序列与GSE83456的交叉点分别为：CXCL10、SLAMF8、C1QB、WARS、IDO1、ANKRD22、CARD17、LAP3、GBP1、ETV7、APOL6、GBP3、TRIM22、RTP4、OAS1、SAMD9L、P2RY14、DDX60、XAF1(见表2)，绘制相应的韦恩图(见图1E)和圈图(见图1F)。

注：A为基因测序中DEGs的火山图；B为基因测序中DEGs的热图；C为GSE83456中DEGs的火山图；D为GSE83456中DEGs的热图；E为显示基因测序和GSE83456中差异基因交集的韦恩图；F为显示交集基因在染色体上相应位置的圈图。

2.2 差异基因的富集分析 富集分析用于分析DEGs的生物学功能，对19个交集DEGs进行GO和KEGG途径富集分析(见图2)。GO分析结果表明，DEGs在抗病毒免疫反应、单核细胞趋化性调节、干扰素-β相关信号通路中显著富集，主要集中在生物学功能(BP)上。KEGG途径分析表明，交集GEGs在NOD样受体信号通路、感染性疾病方面显著富集。

表2 基因测序和GSE83456中DEGs的交集候选靶基因

注：A为GO功能富集柱形图；B为GO功能富集圈图；C为KEGG的通路富集分析气泡图。

2.3 PPI网络构建及相应生物标志物的选择 使用String Online数据库对19个交叉基因进行分析，获得19个交叉基因的相互作用数据，构建PPI网络(见图3A)。统计基因的连接节点数并进行可视化(见图3B)，将数据作为节点导出供进一步分析。导出的数据采用Cytoscape软件进行可视化处理，并使用cytoHubba插件分析MCC的HUB基因(见图3C)，其中GBP1的连接节点数及MCC评分均处于核心地位。根据GBP1的表达量，将GSE83456的样本分为两组：GBP1低表达组和GBP1高表达组，然后采用GSEA对两组数据进行分析(见图3D)。再对GSE83456样本中的EPTB组和HC组进行GSEA分析(见图3E)。GSEA结果显示，高表达GBP1组血标本中显著富集的途径与EPTB组基本一致。免疫系统相关途径如白细胞介素(IL)-6/JAK/STAT3和干扰素-γ(INF-γ)反应等与EPTB高度相关。

注：A为19个交集基因的相互作用；B为PPI网络核心基因的连接节点数；C为19个交集基因中HUB基因的相互作用，颜色越红，MCC评分越高；D为GSE83456中GBP1高表达组和低表达组中Hallmarks基因集的GSEA；E为GSE83456中EPTB组和HC组中Hallmarks基因集的GSEA。

2.4 qRT-PCR检测血清中GBP1的mRNA表达水平 qRT-PCR检测TU患者及HC血清中GBP1的表达水平，结果发现TU患者血清内GBP1的表达水平较HC中高，差异具有统计学意义(P=0.007，见图4A)。治疗2周后，TU患者血清内GBP1的表达水平较治疗前明显下降，差异具有统计学意义(P<0.001，见图4B)。

注：A为qRT-PCR方法检测GBP1在TU患者及HC血清中mRNA的相对表达量；B为qRT-PCR方法检测GBP1在TU患者治疗前及治疗2周后血清中mRNA的相对表达量。

3 讨论

目前临床上报道的结核病以肺结核为主，TU属于EPTB，较为少见。加上临床上多药耐受性结核分枝杆菌层出不穷，TU又常易与其他慢性创面混淆导致误诊漏诊，给TU的防治带来极大的负担，探索其潜在的发病机制和治疗指标是当务之急。

本研究从TU组织基因测序及GEO数据库提取EPTB的原始数据着手，使用生物信息学分析分别挖掘TU组织与TU旁正常组织的差异基因，以及EPTB血液与HC血液之间的DEGs。总共筛选出19个交集的DEGs，均为上调基因。对比这些DEGs在组织和血液中的表达水平，发现组织中的差异表达普遍较血液更为明显。据相关文献[5]报道，超过半数的EPTB患者可引起病灶周围软组织、皮下及皮肤的损伤，最终导致TU的形成。TU形成后，结核分枝杆菌在病灶组织中聚集，通过直接蔓延的方式继续侵犯周围软组织，扩大感染范围，而血液是结核分枝杆菌传播的一个途径，并在全身扩散，推测可能是因为结核分枝杆菌在血液中相较于局部组织含量较低，故因感染引发的变化不明显。

GO分析结果显示，DEGs主要富集在抗病毒免疫反应、单核细胞趋化性调节、干扰素-β(IFN-β)等方面。单核细胞的趋化性是指单核细胞移行至炎症区域或其他趋化刺激物处执行功能，趋化因子在单核细胞趋化的过程中行使调节作用[6]。在结核病中，趋化因子可趋化单核细胞向被结核分枝杆菌感染的部位聚集，参与机体的慢性炎症反应及肉芽肿的形成[7]。IFN-β属于Ⅰ型干扰素家族，常在临床上应用于调节免疫、抑制肿瘤生长和抗病毒等方面[8-9]。Ⅰ型干扰素会促进慢性细菌感染，据相关文献[10-11]报道，IFN-β在结核分枝杆菌感染期间会抑制宿主抗菌免疫反应并促进发病机制。

KEGG富集分析显示，DEGs主要富集在NOD样受体信号通路、感染性疾病等信号通路上。NOD样受体(Nucleotide binding oligomerization domain-like receptors，NLR)为模式识别受体，由介导配体蛋白的C端亮氨酸富集重复结构域(Leucine-rich repeat domain，LRRs)、具有ATP酶结构域的中心核酸结合结构域(Nucleotide-binding domain，NBD)和结合下游分子的N端结构域组成[12]。NOD样受体是结核分枝杆菌激活机体促炎反应的微生物传感器。通过其LRR结构域感知细菌，从而激活NF-κB和MAPK等信号通路，启动机体炎症反应。同时，NOD样受体会激活巨噬细胞及其他免疫细胞，将其招募到感染部位参与慢性炎症的形成。NOD1样受体可以促进趋化因子的产生和中性粒细胞的募集，NOD2样受体可以调节人类巨噬细胞中结核分枝杆菌的促炎反应和宿主反应[13-15]。

对交集基因进行蛋白相互作用分析，筛选出关键基因GBP1。GBP1又称鸟苷酸结合蛋白，是一种在INF-γ刺激下大量表达的基因，属于GTPase家族，在巨噬细胞等多种细胞类型中表达，并在炎症组织中上调[16]。GBP1的高表达可增加结核分枝杆菌感染后细胞因子IL-1β的释放，保护结核分枝杆菌在细胞内存活。GBP1为脂多糖(LPS)的传感器，可作为一种模式识别受体与LPS直接相连，从而激活炎性小体，触发机体炎性反应[17]。同时，GBP1具有抗血管生成活性，一方面通过抑制VEGF-A的表达直接抑制血管的生成，另一方面通过诱导血管抑制细胞因子/趋化因子的释放如CXCL9-11的表达间接抑制血管生成[18-19]。而血管新生与创面的愈合息息相关，是愈合过程中的重要机制，更是创面修复愈合的关键环节。

通过GSEA富集分析对GSE83456芯片进行验证，发现高表达GBP1组血标本中显著富集的途径与EPTB组基本一致。随后，通过qRT-PCR验证人体血液中GBP1的表达水平，证实GBP1在TU患者的血清中高表达并在TU患者接受治疗2周后明显下降。据此推测，TU的创面久不愈合，与GBP1促进炎症反应的发生、结核分枝杆菌感染及抑制机体血管的形成密不可分。抑制GBP1可能加速TU创面的愈合，减轻患者痛苦，提高患者的生存质量。

综上所述，本研究揭示GBP1在TU的发病机制中发挥了重要作用，有可能成为TU的药物靶点和诊断标志物。但是，由于慢性创面久不愈合时，都有着类似的炎性过程，以GBP1能否作为区别TU与其他慢性创面的指标，尚待进一步验证。本研究对TU的分子生物学有了新的认识，也为将来研究TU新的生物标志物及作用机制提供了新思路。

利益冲突：所有作者均声明不存在利益冲突。