活动性肺结核患者与潜伏感染者差异表达基因的生物信息学分析

杨秉芬 翟斐 安红娟 曹志红 刘艳华 王若 程小星

结核病作为世界性传染病已严重威胁人类健康，世界卫生组织(WHO)估计，2017年约有1000万新增结核病患者，死亡人数高达130万人，在致死性传染性疾病中排名第一[1]。流行病学研究表明，全球约有17亿人感染结核分枝杆菌，但其中绝大部分人处于潜伏感染状态而并不发病，约有百分之五到十的感染者最终会发展为活动性结核，在结核病的发病和发展过程中，人体的机体免疫功能发挥了关键作用[2]，但是目前并不是十分清楚具体的分子机制。本文通过转录组测序，获取活动性肺结核患者和潜伏感染者PBMCs差异表达基因，并对这些基因进行生物信息学分析，期望找到活动性肺结核患者与潜伏感染者免疫功能差异的重要分子，为结核病的发病机制研究寻找新的突破点。

资料与方法

一、病例收集

患者样本均来自解放军总医院第八医学中心结核病研究所住院治疗的患者，通过影像学检查与痰菌培养确诊为活动性肺结核患者，并排除其他免疫性疾病或使用免疫抑制剂，共6例，男性2例，女性4例，平均年龄(36.0±4.34)岁。潜伏感染者样本选自总医院第八医学中心体检中心的体检人员，无发热、咳嗽和咳痰等症状，影像学检查肺部无病灶，通过干扰素释放实验，筛选出潜伏感染者，共6例，男性2例，女性4例，平均年龄(30.0±2.59)岁。

二、分离外周血单个核细胞(PBMCs)

采集EDTA抗凝外周血2 mL，按照1∶1的比例使用1640基础培养基 稀释，然后使用等比Ficoll(GE Biosciences, USA)，1000 g，18 min密度梯度离心，中间细胞层为单个核细胞，收集后使用1640基础培养基洗两遍后备用。

三、干扰素释放实验

将PBMCs悬液细胞浓度调整为2.5×106/mL，预包被γ-干扰素抗体的ELISPOT板每孔加入100μL细胞悬液，用结核分枝杆菌特异混合抗原多肽A和B进行刺激，设置阴性对照孔(不加任何抗原)和阳性对照空(PHA)，设置孵箱温度37 ℃、5%CO2培养20 h。用去离子水裂解细胞，扣干水后再用洗涤液洗5遍，拍干后加又入亲和素标记的γ-干扰素抗体，放入培养箱37 ℃孵育60 min。洗去未结合的抗体，然后加入HRP-标记的链霉素，37℃孵育60 min。洗去未结合的链霉素，拍干，加入酶底物显色，室温避光显色，观察斑点形成，用去离子水终止反应，晾干后用免疫斑点计数仪(Cellular Technology Ltd, USA)进行图像和斑点计数。如果阴性对照孔斑点数为0～5点，用测试孔的斑点数减去阴性对照孔斑点数≥6；或者当阴性对照孔斑点数≥6点时，测试孔斑点数≥2倍阴性对照孔斑点数，以上两种情况为阳性反应，可判断为结核菌感染。

四、RNA提取

将细胞用无血清培养基AIM Ⅴ进行重悬，加入到96孔细胞培养板中，按终浓度10 μg/mL加入结核分枝杆菌H37Rv裂解液，37℃、5%CO2培养箱培养过夜。使用TRIzolTMReagent(Invitrogen，USA)进行RNA提取：每个样本加入1 mL TRIzol，0.2 mL氯仿，剧烈振荡15 s，室温放置3 min；12 000 g 4 ℃离心15 min，将上层无色水相分离，并加入异丙醇0.5 mL，上下颠倒混匀，室温内静置10 min；12 000 g 4℃离心10 min，弃上清，使用1 mL无核酸酶水稀释的75%乙醇清洗RNA沉淀；7500 g 4℃离心10 min，弃上清，RNA沉淀置于冰上晾干(5～10 min)，用20μL无核酶水溶解RNA沉淀。使用Agilent 2100对RNA样本进行浓度、纯度和完整性检测，挑选合格的样本，并根据测序样本总量要求混合样本，每个混合样本包含2例女性和1例男性。

五、RNA测序及数据分析

委托华大基因对活动性肺结核患者及潜伏感染者PBMCs进行转录组测序，在得到测序数据raw reads后先去除原始序列中带有的adaptor序列、N的比例大于5%和低质量的reads，获得clean reads，再使用SOAPaligner/SOAP2工具来做参考序列比对。通过统计基因在不同位置的reads数获得reads的分布随机性。基因表达量的计算使用RPKM法(reads per kilobase transcriptome per million mapped reads)，根据基因的RPKM值计算该基因在不同样本间的差异表达倍数。对差异检验的P值作多重假设检验校正，通过控制FDR(False Discovery Rate)来决定p-value的域值。在本次分析中，差异表达基因定义为FDR≤0.001，且倍数差异在2倍以上的基因，共获得差异表达基因98个，其中上调基因67个，下调基因31个。

六、GO及KEGG分析

使用DAVID 6.8(the Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/)对差异表达基因进行GO分析：细胞组分(cell complement)、分子功能(molecular function)和生物进程(biology process)[3]。使用Catespace插件GoClue分析差异表达基因的KEGG通路富集[4]。

七、蛋白相互作用网络分析

使用STRING(Search Tool for the Retrieval of Interacting Genes, http://string db.org/)数据库对差异表达基因进行蛋白相互作用分析，筛选出综合得分(combined scores)大于或等于0.4的相互作用蛋白用于构建蛋白相互作用网络[5]，使用软件Cytoscape 3.7.1(www.cytoscape.org)将相互作用网络绘制出来[6]。用Cytoscape插件MCODE(Molecular Complex Detection, http://apps.cytoscape.org/apps/mcode)对蛋白相互作用网络进行蛋白复合物分析[7]，设定degree cut-off为3，node score cut-off为0.2，k-score为2，maximum depth为100。

结 果

一、差异表达基因的聚类分析

根据华大基因测序结果，对活动性肺结核患者和潜伏感染者差异表达基因使用MeV4.9.0(Multi Experiment Viewer 4.9.0)软件进行聚类分析[8]，结果如图1所示，67个基因在活动性肺结核患者上调显著，31个基因下调显著。

图1 差异表达基因的聚类分析

二、差异表达基因的GO分析和KEGG分析

使用DAVID 6.8对差异表达基因进行GO分析，根据富集分数——-log10(p-value)进行富集条目排序，筛选出富集程度最高的10个条目，结果如图2所示蓝色条形表示“生物进程”富集，在“免疫反应”、“天然免疫反应”和“中性粒细胞趋化”3个词条富集显著；绿色条形表示“细胞组分”富集，在“细胞外空间”、“胞外区”和“外泌体”3个词条富集显著；红色条形表示“分子功能”富集，在“丝氨酸型内肽酶活性”、“糖结合”和“RAGE受体结合”3个词条富集显著。使用Catespace插件GoClue分析差异表达基因的KEGG通路富集如图3所示，在紫色圆点所示“自然杀伤细胞介导的细胞毒作用”、“抗原处理和递呈”，绿色圆点所示的“IL17信号通路”，青色圆点所示的“链球菌感染”4个信号通路显著富集。其中细胞外细胞因子γ-干扰素(interferon-gama, IFNG)，膜蛋白整合素亚单位αM(intergrin subunit alpha M, ITGAM)，胞内蛋白基质金属蛋白酶1(matrix metallopeptidase 1, MMP1)和转录因子FOS连接多个信号通路，是重要的分子，其功能有待进一步研究。

图2 差异表达基因的GO分析 蓝色条形表示“生物进程”富集，绿色条形表示“细胞组分”富集，红色条形表示“分子功能”富集

三、差异表达基因的蛋白相互作用网络分析

使用STRING数据库对差异表达基因进行蛋白相互作用分析，筛选出综合得分大于或等于0.4的相互作用蛋白，用于构建蛋白相互作用网络，使用Cytoscape 3.7.1软件绘制相互作用网络，用NetworkAnalyzer分析网络拓扑学[9]，结果如图4所示：STRING分析共有95个节点(dots)、423对连接(edges)构成网络，平均聚集程度(degree)为8.91。节点的颜色表示基因变化：红色是上调表达基因，其中防御素α1(defensin alpha1, DEFA1)、叶酸受体γ(folate receptor gama, FOLR3)、肽聚糖识别蛋白1(peptidoglycan recognition protein 1, PGLYRP1)和触珠蛋白(haptoglobin, HP)上调倍数(fold change)的大于32，绿色是下调表达基因，其中KIR3DL1、KIR3DS1和KIR2DL4下调较明显，这些基因与NK细胞杀伤作用相关。节点的大小表示节点聚集程度：其中中性粒细胞表达弹性蛋白酶(elastase, neutrophil expresse，ELANE)(degree=32)，ITGAM(degree=31)，基质金属蛋白酶9(matrix metallopeptidase，MMP9)(degree=30)和 精氨酸酶1(arginase 1，ARG1)(degree=26)聚集程度最高。连接线的粗细表示相互作用的综合得分。

图3 差异表达基因的KEGG分析 不同颜色表示不同的KEGG信号通路富集，大的圆点是信号通路，由大到小表示富集程度的变化，最小的圆点是代表性基因

图4 差异表达基因的蛋白相互作用网络分析 圆点的颜色表示基因表达差异，红色表示基因表达上调，绿色表示基因表达下调；圆点的大小表示基因在分子网络中的聚集程度；连接线的粗细表示相互作用的综合得分

图5 差异表达基因的蛋白复合物分析 使用Cytoscape插件MCODE对蛋白相互作用网络进行蛋白复合物分析获得A、B、C三个复合物。

四、差异表达基因的蛋白复合物分析

使用Cytoscape插件MCODE对蛋白相互作用网络进行蛋白复合物分析，设定degree cut-off为3，node score cut-off为0.2，k-score为2，maximum depth为100，结果如图5所示，共有3个蛋白复合物，其中复合物1由上调表达基因PGLYRP1、FOLR3和HP等组成，复合物2上调表达基因CD177、MS4A3和CEACAM1等组成，复合物3由下调表达基因KIR3DL1、KIR3DL2和KIR3DL3等组成，复合物3组成基因与NK细胞介导的杀伤作用相关，说明活动性肺结核患者NK细胞可能比潜伏感染者有显著功能异常，提示我们可以研究NK细胞的功能在结核病发病中的作用。

讨 论

机体的免疫功能在结核病的发病和发展中起了关键的作用[2]，但是具体的分子机制目前并不是十分清楚。生物信息学的发展，特别是一些数据库的开放为科研人员提供了平台，科研人员可以结合自己的关注点进行分析，生成的可视化的结果能更直观地找到新的研究切入点。本文目的旨在通过转录组测序，获取活动性肺结核患者和潜伏感染者PBMCs差异表达基因，并对这些基因进行生物信息学分析，结果显示GO分析生物进程富集的词条是“免疫反应”、“天然免疫反应”和“中性粒细胞趋化”。KEGG分析结果显示富集词条是“自然杀伤细胞介导的细胞毒作用”、“抗原处理和递呈”和“IL17信号通路”等信号通路，其中IFNG、ITGAM、MMP1和FOS基因连接多个信号通路。构建蛋白相互作用网络，筛选到聚集程度高的分子ELANE、ITGAM、MMP9和ARG1。本文测序使用的PBMCs在体外经过结核抗原刺激，活动性结核病患者和潜伏感染者对结核抗原的反应性得到放大，以上的结果提示结核病患者PBMCs在天然免疫如“自然杀伤细胞”和“抗原递呈”与潜伏感染者存在显著差异。筛选到重要的节点蛋白，其中IFN-γ在抗结核免疫具有重要的作用，IFN-γ和TNF-α释放实验可用于鉴别潜伏感染和活动性肺结核，IFN-γ促进结核菌感染的巨噬细胞形成脂滴，加强宿主的保护性免疫反应[10-11]，MMP1和MMP9在结核病的肺损伤中具有重要的作用，可作为结核病诊断的标志物，机制研究发现miR-223可以下调MMP1和MMP9的表达[12-14]。在寄生虫和结核菌共感染模型，ARG1表达上调引起炎症增强和结核病情加重，Ⅰ型干扰素可以通过诱导NOS2和抑制ARG1起到抗结核的作用，又有研究表明抑制肉芽肿中缺氧巨噬细胞ARG1表达可导致肉芽肿病理变化增强和荷菌量增加[15-17]，这些基因的差异与文献报道相符，说明转录组测序和信息学分析结果的可靠。而我们发现新的关键分子ELANE、ITGAM和FOS，其在结核病免疫中的作用知之甚少，为我们下一步研究提供依据。而蛋白复合物分析结果显示由上调基因组成的复合物1和复合物2、由下调基因组成的复合物3，值得注意的是模块3中KIR3DL1、KIR3DL2和KIR3DL3等基因与NK细胞介导的杀伤作用相关，最近多中心研究表明NK细胞亚群在活动性结核病患者和潜伏感染者存在显著差异，并且能预示疾病的进展和抗结核治疗的效果[18]，我们分析获得的功能分子KIR3DL1、KIR3DL2和KIR3DL3为NK细胞在结核病发病中的作用机制研究提供依据。