王益民 苏小艺 苏文革 彭伟 李焱 王怡斐
(1.山东中医药大学第一临床医学院,济南 250014;2.山东中医药大学附属医院,济南 250014)
高血压肾病(hypertensive nephropathy,HN)由慢性高血压引起,是终末期肾病的主要病因之一。高血压是一种容易导致靶器官损伤的疾病,肾脏损伤是其重要的一部分,包括肾小球硬化和萎缩、足细胞丢失以及肾小管间质纤维化等[1-2]。然而,有研究发现即使血压降至目标水平,也只能减缓而不能阻止HN的进展[3]。肾小管间质部分占肾脏总质量的95%,间质炎症和纤维化是导致肾功能下降的重要原因,也是肾硬化症的主要特征[4-5]。大量研究已证实免疫系统在高血压的发病机制中发挥关键作用,淋巴细胞是高血压发生和靶器官损伤的重要参与者[6]。目前尚无反映高血压肾损伤的早期诊断标志物能对HN的发生起到预警作用。基于此,本研究拟以肾小管间质组织为研究对象,进一步探讨HN发病的相关机制,并发掘其潜在的诊断标志物;并利用评估22种免疫细胞相对含量的CIBERSORT算法对HN肾小管间质组织的免疫细胞浸润情况进行分析研究,以揭示免疫细胞浸润模式。
1.1 GEO数据集的获取 从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)下载HN患者和健康人(Normal)肾小管间质组织样本的数据集,GSE编号分别为GSE37455、GSE99325和GSE104954[4,7-8]。其中GSE37455和GSE104954均包括20例HN患者和21例健康人肾小管间质组织样本,GSE99325包括20例HN患者和4例健康人肾小管间质组织样本。所选数据集均具有完整的基因表达谱,同时不存在道德问题、伦理问题和其他利益冲突。本研究中实验组和验证组数据集选择如表1所示。
表1 GEO数据集Tab.1 GEO datasets
1.2 数据处理和差异表达基因(differentially expressed genes,DEGs)的筛选 根据探针注释文件将每个数据集中的探针名称转换为基因名称,应用R软件(版本4.1.2)中“SVA”包的combat功能消除批次效应,使用R软件的“limma”包对实验组中的40个HN样本和24个Normal样本进行差异表达分析。调整后的假阳性率P<0.05和|log2FC|>0.5的样本作为DEGs的阈值。
1.3 功能富集分析 为进一步探索DEGs的生物学意义,采用Metascape[9]在线分析平台(https://metascape.org/gp/index.html)对DEGs进行功能富集分析。将P<0.01、最小计数为3,且富集因子>1.5的项分组到聚类中,相似性>0.3的项由边连接,并呈现为网络图。
1.4 机器学习筛选特征基因 本研究采用3种机器学习算法进一步筛选DEGs。分别是最小绝对收缩和选择算子(LASSO)、支持向量机递归特征消除法(SVM-RFE)和随机森林(RF),在R软件中分别用“glmnet”包、“e1071”包和“randomForest”包实现,然后将3种方法筛选的结果取交集,进而确定HN的特征基因。
1.5 特征基因的诊断价值 分别使用实验组和验证组数据将得到的特征基因进行验证,并绘制受试者工作特征(receiver operating characteristic,ROC)曲线,曲线下面积(area under curve,AUC)值用于确定特征基因在HN和Normal样本中的诊断有效性。
1.6 免疫细胞浸润分析 将基因表达矩阵数据上传至CIBERSORT,筛选出P<0.05的样本,得到免疫细胞浸润矩阵。然后利用R软件中“corrplot”包绘制相关性热图,将22种免疫细胞浸润水平的相关性可视化;“ggplot2”包绘制小提琴图,将免疫细胞浸润的差异可视化。
1.7 特征基因与免疫细胞浸润水平的相关性 为探讨特征基因与免疫细胞浸润水平的相关性,采用R软件进行Spearman等级相关分析,并用“ggplot2”包进行可视化。
2.1 DEGs的结果 本研究对实验组40例HN和24例Normal样本进行了分析,在去除批次效应后,使用“limma”包对数据集进行分析,共获得277个DEGs,其中128个DEGs下调,149个DEGs上调,DEGs的结果显示在热图和火山图中(图1、图2)。
图1 DEGs热图Fig.1 Heat map of DEGs
图2 DEGs火山图Fig.2 Volcano map of DEGs
2.2 功能富集分析结果 在Metascape在线分析平台上进行了DEGs的功能富集分析,图3、图4列出了排在前20的具有显著意义的富集结果,并呈现出关系网络,其中基因本体论(gene ontology,GO)生物过程主要与先天免疫反应、对细胞因子的反应、体液免疫反应、细胞因子产生的正向调节、炎症反应和免疫系统过程的负向调节等有关。反应组基因集(Reactome Gene Sets)主要与免疫系统中的细胞因子信号传导、干扰素信号传导和中性粒细胞脱颗粒等有关。维基通路(WikiPathways)主要与同种异体移植排斥等有关。
图3 DEGs富集项的条形图Fig.3 Bar graph of enriched terms of DEGs
图4 DEGs富集项的网络图Fig.4 Network diagram of enrichment terms of DEGs
2.3 特征基因筛选结果 使用LASSO回归算法确定17个基因作为HN的候选特征基因。使用SVMRFE算法确定40个基因作为HN的候选特征基因。使用RF算法得出最佳决策树数目为136个,在DEGs中挑选出评分最高的30个基因作为候选特征基因。取三种不同算法得到基因的重叠部分。最后获得了3个特征基因:CISH、GADD45A和ZFP36(图5)。
图5 三种机器学习算法筛选结果图Fig.5 Screening results of three machine learning algorithms
2.4 特征基因的诊断价值结果 如图6所示,3个特征基因在实验组和验证组中均表现出良好的诊断价值,在实验组数据集中,CISH的AUC值为0.963(95%CI:0.915~0.993),GADD45A的AUC值为0.953(95%CI:0.893~0.994),ZFP36的AUC值为0.916(95%CI:0.835~0.972)。在验证组数据集中,CISH的AUC值为0.793(95%CI:0.648~0.923),GADD45A的AUC值为0.775(95%CI:0.616~0.895),ZFP36的AUC值为0.755(95%CI:0.595~0.889)。
图6 验证3个特征基因诊断价值的ROC曲线Fig.6 ROC curve for verifying diagnostic value of three characteristic genes
2.5 免疫细胞浸润分析结果 得到的22种免疫细胞浸润水平相关性热图如图7,其中单核细胞和活化肥大细胞存在明显的正相关关系(R=0.41),幼稚CD4+T细胞和活化肥大细胞存在明显的正相关关系(R=0.37)。幼稚B细胞和记忆B细胞存在明显的负相关关系(R=-0.64),活化肥大细胞和静息肥大细胞存在明显的负相关关系(R=-0.62),浆细胞和M2巨噬细胞存在明显的负相关关系(R=-0.45)。免疫细胞浸润差异的小提琴图显示,与Normal组相比,HN组中肾小管间质组织调节性T细胞和M1巨噬细胞浸润较多(图8)。
图7 HN中免疫细胞间的相关性分析热图Fig.7 Analysis of correlation between immune cells in HN by heat map
图8 HN与Normal组间免疫细胞浸润差异的小提琴图Fig.8 Violin diagram of the difference of immune cell infiltration between HN and Normal groups
2.6 特征基因与免疫细胞浸润水平的相关性 相关性分析结果显示(图9),CISH与静息树突状细胞呈正相关(R=0.33,P=0.017),与M1巨噬细胞呈负相关(R=-0.41,P=0.002 6)。GADD45A与浆细胞(R=0.59,P=5.6e-06)、活化CD4+记忆T细胞(R=0.29,P=0.035)和幼稚CD4+T细胞(R=0.28,P=0.046)呈正相关,与M1巨噬细胞(R=-0.38,P=0.005)、M2巨噬细胞(R=-0.30,P=0.029)和调节性T细胞(R=-0.27,P=0.048)呈负相关。ZFP36与静息树突状细胞(R=0.37,P=0.007 1)和活化肥大细胞(R=0.37,P=0.006 1)呈正相关,与M1巨噬细胞(R=-0.37,P=0.006 9)和静息肥大细胞(R=-0.36,P=0.007 6)呈负相关。
图9 3个特征基因与免疫细胞浸润水平的相关性Fig.9 Relationship between three characteristic genes and level of immune cell infiltration
目前HN的治疗仍以控制血压为主,尚无特异性治疗方式。由于人口老龄化和心血管疾病存活率的提高,预计HN的发病率在未来几十年将进一步增加,这也对HN的防治提出了新的挑战。本研究通过分析包含40个HN和24个Normal肾小管间质组织的基因表达谱,得到277个DEGs。利用Metascape在线分析平台进行生物信息学分析,发现DEGs与多种免疫相关过程有关,如先天免疫反应、体液免疫反应、对细胞因子的反应、细胞因子产生的正向调节、免疫系统中的细胞因子信号传导、干扰素信号传导及中性粒细胞脱颗粒等有关。研究发现,先天免疫反应的激活伴随着炎症反应的发展,进而引起纤维化、基质沉积和进行性肾损伤被认为是包括高血压肾损害在内的多种肾脏疾病发病机制的关键因素[10]。体液免疫系统激活与高血压发病之间的关联也得到证实,系统性红斑狼疮(systemic lupus erythematosus,SLE)小鼠体液免疫系统激活和高血压的发病具有明确的因果关系,用小鼠抗CD20抗体以消耗B细胞,可显著减弱自身抗体的产生,并防止SLE小鼠模型高血压的发展,强烈支持体液免疫系统激活可导致高血压的发病[11-12]。多种细胞因子参与HN的发病机制,炎症相关细胞因子如IL-1β、IL-6、单核细胞趋化蛋白-1(MCP-1)、细胞间黏附分子-1(ICAM-1)和TNF-α等参与了包括高血压肾损伤在内的炎症反应[13-14]。纤维化相关细胞因子如TGF-β、结缔组织生长因子(CTGF)等,在HN肾间质纤维化的形成中发挥关键作用[15-17]。由此可见,免疫相关机制在HN的发病中扮演重要角色,为揭示HN的发病机制提供了新方向。
本研究采用三种机器学习算法进一步筛选DEGs。LASSO是一种回归分析算法,其使用正则化来提高预测精度[18]。支持向量机(SVM)是一种有监督的机器学习技术,在分类和预测方面具有优异的性能[19]。RF是通过集成学习的思想将多棵树集成的一种算法,在处理多特征数据方面表现出优异的性能[20]。这三种机器学习算法目前广泛应用于特征基因的筛选[21-22]。本研究基于这三种机器学习算法,最终在DEGs中筛选出3个特征基因,分别是CISH、GADD45A和ZFP36。细胞因子诱导的含SH2蛋白(CISH)属于细胞因子信号传导抑制因子(SOCS)家族成员,SOCS家族包括SOCS1~SOCS7和CISH,最新研究表明,SOCS家族成员可通过介导细胞因子信号传导的负反馈抑制,在先天性和适应性免疫应答中起关键作用[23]。生长停滞和DNA损伤诱导蛋白45α(GADD45A)是生长停滞和DNA损伤诱导蛋白45(GADD45)蛋白家族成员,GADD45A在各种细胞功能中发挥至关重要的作用,包括DNA修复、细胞凋亡和DNA甲基化等[24]。锌指蛋白36(ZFP36)是一种ARE结合蛋白,可促进促炎细胞因子降解,例如TNF-α、粒细胞-巨噬细胞集落刺激因子(GM-CSF)、IL-6和环氧合酶-2(COX-2)等,在限制炎症反应方面具有重要作用[25]。3个特征基因在实验组和验证组均表现出良好的诊断价值,在生理功能上与细胞因子功能的调控、细胞功能的调节以及炎症反应的调节等相关,有望成为HN潜在的诊断标志物和治疗靶点。
免疫细胞浸润的结果显示HN组中肾小管间质组织调节性T细胞(regulatory T cells,Tregs)和M1巨噬细胞浸润较多,提示Tregs和M1巨噬细胞在HN中发挥关键作用。根据目前的研究,Tregs可影响免疫系统的多个环节来预防高血压,以及减轻靶器官损伤[26]。缺乏Tregs会加重血管紧张素Ⅱ(Ang Ⅱ)依赖性高血压,而反复过继转移Tregs则会减弱血压升高[27]。Tregs产生的IL-10可通过减轻血管氧化应激改善高血压微血管内皮功能[28]。Tregs还与补体系统相互作用。有研究表明,特异性靶向Tregs中的补体受体C3aR和C5aR可能是治疗高血压的另一种新方法[29]。M1巨噬细胞是极化巨噬细胞的一种,可产生活性氧和促炎细胞因子加剧炎症,特别是TNF-α和IL-1β,介导了高血压发病中的众多机制[30]。还有研究指出,CD14+M1巨噬细胞通过其强烈表达血管紧张素转换酶促进血压升高,提示它们可能通过RAS系统参与高血压的发生[31]。22种免疫细胞浸润水平的相关性结果显示,多种免疫细胞在HN的发病中具有相关性,单核细胞和活化的肥大细胞存在较强的正相关关系,幼稚B细胞和记忆B细胞存在较强的负相关关系。特征基因与免疫细胞浸润水平的相关性结果显示,3个特征基因与M1巨噬细胞均具有显著的负相关性,CISH和ZFP36与静息树突状细胞均呈正相关,GADD45A与浆细胞、活化CD4+记忆T细胞等呈正相关。目前关于各免疫细胞之间、特征基因与免疫细胞之间相关性的基础研究尚缺乏大规模验证,本研究结果具有一定的参考价值。
本研究也存在一定的局限性,数据来源依赖于GEO数据库、样本数量有限均可能造成分析结果的偏倚,同时还需要进一步的实验验证。总之,本研究利用生物信息学和机器学习的方法对HN的DEGs进行了分析研究,并探讨了HN的免疫细胞浸润机制,对于HN发病机制的探讨、诊断标志物的筛选及治疗均具有参考意义。