基于GEO数据库的狼疮肾炎生物信息学分析及差异表达基因筛选*

陶 丽， 钱诗睿， 苏 华△

华中科技大学同济医学院附属协和医院 1肾内科 2心血管外科，武汉 430022

系统性红斑狼疮(systemic lupus erythematosus，SLE)是一种慢性自身免疫性疾病，可累及多器官、多组织，大约50%的患者会累及肾脏，发展成为狼疮肾炎(lupus nephritis，LN)[1]。LN发病机制复杂，临床表现多样，主要为血尿、蛋白尿、肾功能损伤、急进性肾功能衰竭等[2]。LN患者预后不同，尽管抗炎和免疫抑制疗法有一定的效果，但仍有约10%的LN患者会发展为终末期肾脏病(end stage renal disease，ESRD)，是SLE总体发病率和死亡率的主要危险因素[3]。目前，由于缺乏对分子机制的了解，针对SLE的特异性靶向治疗的发展缓慢。阐明LN的发病机制对于促进靶向治疗研究，进一步改善SLE患者预后，降低病死率具有重要意义。

在基因组水平监测LN的生物学变化是一种很有意义的研究方法。近年来，生物信息学分析已经应用于多种疾病中，可以在极短时间内处理大量的转录组数据并提供相关疾病的有价值信息。分析肾脏组织的基因表达，评估哪些基因对LN更有意义，有助于寻找有效的生物标志物。此外，功能富集分析研究可加深我们对LN的理解。因此，本研究通过生物信息学分析来确定与LN发生发展相关的差异表达基因(differentially expressed genes，DEGs)，为进一步探究LN的致病分子机制提供基础。

1 资料与方法

1.1 研究资料的获取

在美国国立生物技术信息中心(National Center for Biotechnology Information，NCBI)的高通量基因表达(Gene Expression Omnibus，GEO)数据库(http：//www.ncbi.nlm.nih.gov/geo/)中检索得到LN人群研究数据集GSE127797和GSE32591。GSE127797是基于GPL24299构建的芯片，包含LN肾组织(54例)基因表达谱。GSE32591是基于GPL14663构建的芯片，包括正常肾组织(29例)和LN肾组织(64例)基因表达谱。

1.2 DEGs的数据处理

1.2.1 数据预处理 通过GEOquery包(2.56.0版本)下载数据集GSE127797和GSE32591的原始数据，并加载至RStudio(4.0.2版本)软件中。读取各自GPL的探针转换表格，将数据集的ENTREZ ID转换为SYMBOL ID。过滤SYMBOL ID对应多个ENTREZ ID的基因并取表达量的平均值；将数据集的肾小球数据和肾小管数据分开，分别用limma包(3.44.3版本)进行中位数法标准化，以消除非实验误差，使各样本间具有可比性。

1.2.2 质量评估 分别合并2个数据集的肾小球数据和肾小管数据，用R语言的limma包(3.44.3版本)进行批次校正，降低2个数据集因不同实验室、实验批次等原因造成的非实验性误差，使2个数据集的数据具有可比性，并再次用中位数法进行标准化；分别对肾小球和肾小管数据使用ggfortify包(0.4.10版本)进行主成分分析(principal components analysis，PCA)和使用cluster包(2.1.0版本)进行hclust聚类分析。

1.2.3 DEGs选取 用R语言中的limma包(3.44.3版本)进行差异分析，以|log2FoldChange|>1、P值<0.05为截断值挑选DEGs，使用ggplot2包(3.3.2版本)做火山图进行可视化；以|log2FoldChange|>1.5、P值<0.05为截断值挑选DEGs，使用pheatmap包(1.0.12版本)做热图进行可视化。根据log2FoldChange值排序，分别取上调的前250个基因和下调的前150个基因作为肾小球数据和肾小管数据的DEGs，取肾小球数据和肾小管数据交集后的共同DEGs。使用VennDiagram包(1.6.20版本)绘制韦恩图。

1.3 功能富集分析

使用DAVID数据库(https：//david.ncifcrf.gov/tools.jsp)对筛选出的共同DEGs进行基因本体(Gene Ontology，GO)富集分析和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)通路富集分析。用R软件的clusterProfiler包(3.16.1版本)和org.Hs.eg.db包(3.11.4版本)处理数据。对共同DEGs进行GO功能富集分析，主要富集在生物学过程(biological process，BP)、细胞组成(cellular component，CC)、分子功能(molecular function，MF)这3个模块中。利用R软件中的GOplot包(1.0.2版本)使GO富集结果可视化。对共同DEGs进行KEGG通路富集分析，使用R语言ggplot2包(3.3.2版本)绘制气泡图。以P<0.05为差异具有统计学意义。

1.4 建立蛋白质互作网络和挑选关键DEGs

将DEGs的基因名上传到STRING网站(https：//string-db.org/)，设置minimum required interaction score=0.4，然后将计算结果导入Cytoscape(3.8.0版本)中，构建LN的DEGs蛋白质相互作用(protein-protein interaction，PPI)网络。使用Cytoscape软件中的MCODE插件，设置degree cutoff=2，node score cutoff=0.2，k-core=2，max.depth=100，筛选出PPI网络中生物联系最为紧密的子网络。使用Cytoscape软件中的Cytohubba插件，采用Degree算法筛选出PPI网络中连结度最高的前20个关键DEGs，颜色越红表示得分越高。最后，得分前10且存在于子网络中的基因即为最终筛选出的与LN生物学行为密切相关的关键DEGs。

2 结果

2.1 DEGs的选取

基于R语言对数据集GSE127797和GSE32591进行数据的标准化及批次校正处理后(图1)，各样本中位数基本在一条水平线上，说明样本间归一化程度好。之后分别对肾小球和肾小管数据进行质量评估，绘制出PCA和hclust聚类分析图(图2)。PCA图肾小球和肾小管的实验组和对照组数据呈离散分布，综合hclust聚类分析图结果，表明实验组与对照组数据具有很好的区分度。根据|log2FoldChange|和P值对肾小管和肾小球数据进行初步筛选DEGs后的火山图和热图见图3A～3D。肾小球组共筛选出400个(上调250个+下调150个)DEGs，肾小管组共筛选出400个(上调250个+下调150个)DEGs。用韦恩图筛选出这2组数据的共同DEGs，共114个，包含64个上调基因和50个下调基因(图3E)。分别挑取共同DEGs数据集中上调和下调的前10位基因展示，见表1。

A：GSE127797肾小球数据标准化；B：GSE32591肾小球数据标准化；C：GSE127797肾小管数据标准化；D：GSE32591肾小管数据标准化；E：GSE127797和GSE32591肾小球数据合并，进行批次校正及标准化；F：GSE127797和GSE32591肾小管数据合并，进行批次校正及标准化图1 数据预处理Fig.1 Data preprocessing

A：肾小球数据PCA图；B：肾小管数据PCA图；C：肾小球数据hclust图；D：肾小管数据hclust图图2 数据质量评估Fig.2 Data quality assessment

A：肾小球数据火山图；B：肾小管数据火山图；C：肾小球数据热图；D：肾小管数据热图；E：共同DEGs的韦恩图图3 选取DEGsFig.3 DEGs selection

表1 前10个显著上调或下调DEGsTable 1 The top 10 significantly upregulated or downregulated DEGs

2.2 GO富集分析和KEGG信号通路分析结果

通过DAVID数据库对共同DEGs进行GO和KEGG分析(图4)。GO富集分析以人源基因为背景，对共同DEGs进行生物学功能注释。结果显示，BP模块中，共同DEGs主要参与病毒防御反应、Ⅰ型干扰素(type Ⅰ interferon，IFN-Ⅰ)信号通路、病毒基因组复制的调控等生物通路；CC模块中，共同DEGs参与的细胞组分包括胶原蛋白三聚物、血小板α颗粒、胞质囊腔、细胞器外膜等；MF模块中，共同DEGs主要与血小板衍生生长因子结合、细胞外基质结构成分、生长因子结合相关。利用KEGG通路数据库进行信号通路的富集，寻找共同DEGs所富集的信号通路。结果表明，共同DEGs与甲型流感病毒感染、麻疹病毒感染、丙型肝炎病毒感染、百日咳杆菌感染、补体和凝血级联、蛋白质消化吸收、金黄色葡萄球菌感染等通路有关。

A：GO功能富集分析；B：KEGG通路富集分析图4 共同DEGs的GO和KEGG分析Fig.4 GO and KEGG analysis of common DEGs

2.3 PPI网络的构建及关键DEGs的筛选

将上述筛选得到的114个共同DEGs导入STRING数据库并利用Cytoscape软件构建共同DEGs之间的蛋白质相互作用网络，得到PPI网络图(图5A)。PPI网络图中每个节点代表由1个蛋白质编码基因位点产生的所有蛋白质，节点内为该蛋白质的三维结构。两节点间的连线代表蛋白质的相互联系，不同颜色代表不同的生物学意义。使用Cytoscape软件中的MCODE插件筛选出PPI网络图中生物联系最为紧密的子网络(图5B)。同时使用Cytoscape软件中的Cytohubba插件对PPI网络图中每个节点的邻接数进行计算，筛选出得分前20的关键DEGs(图5C)。

A：PPI网络图；B：PPI网络图中生物联系最为紧密的子网络图；C：得分前20的关键DEGs；D：前10位的关键DEGs蛋白互作网图5 建立PPI网络图和挑选关键DEGsFig.5 Construction of PPI network map and selection of key DEGs

最后，得分前10且存在于子网络中的关键DEGs(图5D和表2)：ISG15、MX1、OAS1、MX2、IFIH1、GBP1、IFIT3、OAS2、IFIT1、IFI44，即为最终筛选出的与LN生物学行为密切相关的关键DEGs，颜色越红表示得分越高。

表2 10个关键DEGs及其功能Table 2 The 10 key DEGs and their functions

3 讨论

SLE的确切病因尚未完全阐明，目前认为基因、环境因素刺激(如感染、药物、紫外线、饮食等)及表观遗传修饰异常等均与SLE发病有关[4]。SLE呈现广泛的临床和免疫学表现，其中LN是导致患者致残和死亡的最常见原因。LN的发病机制涉及多种致病途径，包括异常的细胞凋亡、自身抗体产生，免疫复合物沉积和补体激活[5]；并且不同的病理机制可相互作用，最终导致肾脏受损。

本文通过对GEO数据库中正常肾组织和LN肾组织的基因芯片数据集进行差异分析与GO富集分析发现，LN肾小球与肾小管共同DEGs主要参与病毒防御反应及IFN-Ⅰ信号通路等；KEGG分析结果表明，共同DEGs与甲型流感病毒感染、麻疹病毒感染等通路有关。

病毒感染和自身免疫性疾病的因果联系已经在大量临床研究中得到了验证：在5种常见的人类病毒(甲型流感病毒、博尔纳病病毒、麻疹病毒、腮腺炎病毒和风疹病毒)和人类蛋白质组之间出现了大量的肽段重叠，这种序列相似性或许可以解释病毒感染或免疫应答过程中的自身免疫交叉反应[6]。多项研究证实感染性因素可诱导及促进SLE的进展，SLE与病毒感染存在如下可能的关联[7]：①SLE可以发生在慢性病毒感染期间；②SLE的发病可能与病毒感染同时发生；③特定的病毒感染可以触发SLE。如，Epstein-Barr(EB)病毒、细小病毒B19和人内源性逆转录病毒参与SLE的发病过程；除此之外，甲型流感、麻疹、丙肝病毒等也可能与SLE的发病有关。韩国最近的一项研究表明，季节性流感爆发与SLE发病之间存在显著相关性[8]。

干扰素(interferon，IFN)是重要的免疫系统介质，可通过调控树突状细胞、淋巴细胞、自然杀伤细胞以及单核巨噬细胞启动免疫反应与放大组织损伤。IFN-Ⅰ包括IFNα，IFNβ，IFNτ，IFNω和IFNκ，它们在天然免疫和病毒防御中发挥重要作用，亦是SLE发生和发展的核心因素。Ⅱ型IFN以IFNγ为代表，是适应性免疫反应效应机制的关键因子。SLE患者体内表达上调的细胞游离核酸(DNA/RNA)，尤其是双链DNA(double-stranded DNA，dsDNA)，是SLE和LN发病机制中的关键物质；另外，细胞游离DNA/RNA是IFN-Ⅰ最有效的诱导因子，故IFN-Ⅰ是LN发病机制中的重要环节。除循环免疫细胞外，肾脏固有细胞亦是IFN-Ⅰ的主要来源，如DNA/RNA免疫复合物诱导足细胞产生IFN-β[9]。综上所述，GO分析和KEGG分析提示的结果与既往研究报道相符合，并且进一步加深了我们对LN发病机制的理解。

通过构建PPI网络图并进行算法分析和筛选，最终获得10个最重要的DEGs：ISG15、MX1、OAS1、MX2、IFIH1、GBP1、IFIT3、OAS2、IFIT1、IFI44，这10个DEGs均受IFN的调控。

ISG15属于干扰素诱导基因(interferon-inducible gene，IFIG)，Carrillo-Vázquez等[10]研究发现，相比于健康对照组，SLE患者中性粒细胞胞外诱捕网(neutrophil extracellular traps，NETs)中ISG15的表达升高，且SLE患者NETs刺激外周血单核细胞(peripheral blood mononuclear cell，PBMC)产生IFNγ的能力明显增强。此外，SLE患者血液中ISG15 mRNA水平更高，且与治疗前SLE的活动相关；ISG15的表达水平亦与SLE患者体内淋巴细胞的减少有关，提示ISGl5可能参与淋巴细胞的凋亡[11]。Han等[12]研究发现SLE患者的EB病毒感染率较高，EB病毒的潜伏膜蛋白1(latent membrane protein 1，LMP1)的表达较高，并且ISG15与LMP1的表达呈正相关，提示EB病毒的LMP1可能通过激活IFN-Ⅰ途径参与SLE的发生和发展。使用JAK抑制剂托法替尼(tofacitinib，TOFA)通过JAK-STAT途径控制IFN的信号传导可降低NZB/NZW F1小鼠血清中抗dsDNA抗体水平，减少蛋白尿并改善肾炎；且给予TOFA和地塞米松(dexamethasone，DEXA)处理后，SLE易感小鼠CD4+T细胞中ISG15的表达降低[13]。靶向IFN信号传导或许可用于开发新的SLE特异性治疗策略。还有研究发现，姥鲛烷诱导的狼疮小鼠ISG15的表达上调；姥鲛烷诱导的miR155缺陷型狼疮小鼠的血清自身抗体水平下降、肾脏损伤程度减轻且ISGl5表达减少[14]。上述研究均提示ISG15与LN的发病相关。

MX1是一种IFIG。Shimizu等[15]使用酶联免疫吸附测定(enzyme-linked immunosorbent assay，ELISA)发现，与IgA肾病(IgA nephropathy，IgAN)患者、ANCA相关血管炎(ANCA-associated vasculitis，AAV)患者和健康对照者相比，SLE患者外周血中的MX1蛋白浓度显著升高。使用免疫组织化学方法发现，相比于IgAN和AAV肾标本，LN的肾小球和肾小管组织中MX1阳性区域显著增多。与未经免疫抑制剂治疗的患者相比，经免疫抑制剂治疗的LN患者肾组织MX1蛋白水平较低。人类MX2是IFN诱导的GTPase超家族的成员。MX2基因编码的蛋白具有细胞核和细胞质形式，其中MX2核蛋白以颗粒状分布于核膜下的异染色质区，可抑制许多RNA和DNA病毒的初期复制。MX2最接近的家族成员是人类MX1(63%氨基酸序列相同)[16]，因此MX2亦可能与LN的发病相关。

OAS1和OAS2受IFN调控，属于2′-5′-寡腺苷酸合成酶家族成员且均位于12号染色体上。OAS1和OAS2是OAS不同亚型，是机体对病毒感染的固有免疫反应中的必需蛋白。OAS家族蛋白的抗病毒活性主要归因于合成2′-5′-寡腺苷酸的能力，并激活潜在的RNase L导致蛋白质合成受到抑制[17]。Ye等[18]研究发现狼疮活动患者中OAS1和OAS2的mRNA表达高于感染患者和正常对照组。Landolt-Marticorena等[19]的研究显示SLE患者外周血中OAS1和B细胞活化因子(B cell activation factor，BAFF)表达升高，且两者之间存在相关性。与对照组比较，SLE患者外周血中的B细胞异常增生，提示OAS1的高表达可能与B细胞的异常增生有关。

IFIH1是一种IFIG，可以编码黑色素瘤分化相关基因5(melanoma differentiation associated gene 5，MDA5)。MDA5参与病毒感染后干扰素应答调节的过程。Su等[20]发现SLE患者的疾病活动性与PBMC中MDA5的水平呈负相关。Munroe等[21]对SLE患者血液中IFIH1与炎症介质、自身抗体的相关性进行了评估，发现IFIH1基因与炎性介质白介素-6(interleukin-6，IL-6)、干扰素诱导蛋白10(interferon-induced protein 10，IP-10)及自身抗体的表达密切相关，表明IFIH1可能通过调节炎症反应参与SLE发病。此外，IFIH1 G821S错义突变的小鼠可出现狼疮样症状，研究者在这些小鼠血清中检测到抗核抗体和抗dsDNA抗体的表达，同时在肾脏组织中观察到免疫球蛋白和补体的沉积；肾脏中包括IFNβ、IL-6和趋化因子配体10(CXC chemokine ligand 10，CXCL10)在内的炎性细胞因子和趋化因子显著上调，IFNβ和IL-6则在全身器官中表达上调；错义突变的IFIH1可能通过激活IFN-Ⅰ途径触发机体的自身免疫反应[22]。上述研究表明，IFIH1可能参与LN的肾脏损害。

GBP1是一种干扰素刺激基因(interferon-stimulated gene，ISG)。IFN刺激，尤其是IFNγ刺激，会促进该基因的过度表达[23]。在许多类型的细胞，如内皮细胞和单核细胞中，GBP1的表达受到IFNγ的强烈刺激，并且GBP1在炎症环境中抑制细胞增殖。在免疫应答中，GBP1的活性对于细胞内病原体感染的自噬体的成熟和细胞对病原体相关分子模式的反应至关重要[24]。基于GBP1已知的生物学功能，GBP1是否参与LN的致病过程值得我们进一步探索。

IFIT1、IFIT3都是IFIG，属于干扰素诱导的四肽重复蛋白家族。Landolt-Marticorena等[25]发现，与对照组比较，SLE患者外周血的IFIT1表达显著上升，并与疾病的高度活动性相关。Hu等[26]使用MRL/lpr狼疮小鼠模型研究发现IFIT1的表达与MRL/lpr小鼠肾组织中F-actin、Nephrin和Podocin3种足细胞蛋白的表达呈负相关；随着IFIT1表达的增加，MRL/lpr狼疮小鼠肾组织中足细胞大量丢失；上调IFIT1的表达促进足细胞损伤，加重LN肾损伤。Wang等[27]研究发现，与健康对照组比较，IFIT3在SLE患者PBMC中表达明显升高，且与cGAS/STING信号通路的活性呈正相关；提示IFIT3可作为一种新的治疗靶点，用于阻断SLE患者通过cGAS/STING信号通路产生IFN-Ⅰ和其他促炎细胞因子。

IFI44是一种IFIG。有研究表明在SLE患者的PBMC中IFI44的表达显著增高，并且IFNα优先于IFNγ诱导IFI44的表达，这提示IFI44参与IFN-Ⅰ信号通路介导的SLE的致病过程[28]。Shen等[29]最近的研究发现，与健康对照相比，LN患者血清中IFI44显著上调。此外，活动性LN患者血清中的IFI44显著高于非活动性LN患者。

综上所述，本研究采用生物信息学分析方法挖掘了与LN相关的DEGs，并经GO富集分析和KEGG通路分析显示病毒防御反应、IFN-Ⅰ信号通路、病毒基因组复制的调控等生物学过程可能参与LN的发病。此外，通过构建PPI网络并利用算法分析获得10个关键DEGs，这10个DEGs都受到IFN的调控；其中ISG15、MX1、OAS1、IFIH1、IFIT3、OAS2、IFIT1、IFI44等基因在LN中有少量研究，但它们参与LN的具体致病机制尚未阐明；MX2、GBP1基因与LN的具体关系尚未见明确报道。本研究为进一步探究LN致病相关分子机制、发掘潜在治疗靶点提供了新的理论依据与方向。本研究的不足之处在于，我们未能通过实验进一步验证这些关键DEGs在LN发生发展过程中的表达变化及具体作用。因此，未来我们将通过后续的实验研究来验证这些关键DEGs的功能。