COPD 差异表达基因的生物信息学分析及在LUSC样本中的表达

2022-12-17刘艳玲童明琼王晓玥

生物信息学 2022年4期

刘艳玲，刘静，童明琼，范娜，王晓玥，孙婉

(德州学院医药与护理学院，山东德州 253023)

慢性阻塞性肺病(Chronic Obstructive Pulmonary Disease，COPD)是一种典型的具有进行性、不可逆的气流限制和反复呼吸道症状的慢性呼吸系统疾病，致残率和病死率很高，全球40岁以上发病率已高达9-10%[1-2]。COPD并非单一疾病，它不仅损伤肺组织，同时还会损伤包括心脏、骨骼、免疫系统、消化系统和中枢神经系统在内的肺外器官和组织，COPD目前已成为全球第三大死因[3-4]。

COPD患者肺部的慢性炎症反应可导致肺实质组织受损，并破坏正常的修复功能和防御机制，导致小气道纤维化，最终产生气体陷闭和进行性气流受限，成为导致气道阻力增加的主要因素[5-6]。由此可见，COPD与小气道病变息息相关，尽早干预轻中度COPD患者的小气道病变，对于延缓病程具有重要意义。研究表明，被诊断为肺癌的患者中40%-70%患有COPD[7]，而COPD患者的肺癌发病率为16.7‰，且最长见的病理类型是肺鳞状细胞癌(Lung Squamous Cell Carcinoma，LUSC)[8]。对于吸烟者来说，COPD患者患肺癌的几率是非COPD患者的4～6倍[9]，并且轻中度COPD患者患肺癌的几率较重度COPD更高，这可能与轻中度COPD患者存在炎症和氧化应激增高有关[10]。目前，COPD与肺癌共存机制尚不清楚，如果能够探寻COPD合并肺癌的预测因子，研究两者共同的发病机制，发现新的治疗靶点，可为COPD的预防、治疗和肺癌的早期诊断、降低死亡率提供理论依据。

目前，关于COPD的生物信息学研究还较少，王小乐等通过对COPD患者肺泡巨噬细胞的芯片数据集进行生物信息学分析，共获得43个差异表达基因，且与免疫和炎症的发生和发展有关[11]。谷雨等对COPD 患者肺组织的芯片数据进行了生物信息学分析，找到了FGG, FGA, IL-6, SERPINE1和SPP1 5个Hub基因，根据调控途径的进一步分析发现差异基因表达可能与代谢及细胞生物功能异常有关[12]。

本研究采用生物信息学方法，从GEO(Gene Expression Omnibus)数据库中筛选3个COPD芯片数据集，挖掘COPD患者小气道上皮细胞(Small airway epithelium，SAEC)的差异表达基因(Differentially Expressed Genes，DEGs)以及潜在的生物标记物，并通过基因本体(Gene Ontology，GO)和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)富集分析预测DEGs的功能及参与的代谢途径。通过对DEGs进行COPD的蛋白质网络互作分析，构建PPI(Protein-Protein Interaction)网络，并使用Cytoscape软件中的CytoHubba和MCODE筛选子模块和Hub基因。最后将10个Hub基因通过TCGA和GTEx数据库验证其在LUSC样本中的差异表达情况及差异基因间的相关性，为确定COPD 的分子标记物及COPD与LUSC共存的DEGs奠定一定的基础。

1 材料与方法

1.1 材料

使用“COPD”作为关键词，在GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)中搜索并选择了3个芯片数据集：GSE11784、GSE11906和GSE20257，均来自于GPL570平台(Affymetrix Human Genome U133 Plus 2.0 Array，Affymetrix，Inc)。由于芯片数据都出自威尔·康奈尔医学院遗传医学系实验室，其中的部分样本在不同芯片间存在重复，去除重复样本后，最终GSE11784纳入53例正常对照和14例COPD患者，GSE11906纳入46例正常对照和15例COPD患者，GSE20257纳入22例正常对照和9例COPD患者(见表1)，3个芯片共纳入159例样本，其中包含121例正常对照和38例COPD患者。

表1 GEO数据库COPD芯片数据集Table 1 COPD microarray datasets from GEO database

1.2 方法

1.2.1 微阵列数据处理

分别对3个数据集探针的矩阵文件取log2，然后使用Perl 5.0软件(http://www.perl.org)对数据进行注释、合并，再通过R软件(4.0.3版本)的sva和limma包对合并后的数据进行批次校正。

1.2.2 DEGs分析

使用R软件中的Limma包对正常样本和COPD样本进行比较，在校正后的P值(AdjP)<0.05的条件下鉴定DEGs。

1.2.3 Go和KEGG 富集分析

使用R软件中的ClusterProfiter和org.Hs.eg.db包进行GO和KEGG富集分析，以P<0.05定义差异具有统计学意义。

1.2.4 PPI网络分析

使用在线网站STRING(http://string-db.org)分析DEGs的相互作用关系，选择最小交互作用得分>0.4的蛋白质对构建PPI网络。然后分别使用Cytoscape中的MCODE(degreecut≥2，nodescorecut≥0.2，K-core≥2，maxdepth=100)和CytoHubba插件筛选重要的子模块和Hub基因。

1.2.5 Hub基因在TCGA和GTEx数据库中的显著性验证

使用在线网站GEPIA(Gene Expression Profiling Interactive Analysis) (http://gepia.cancer-pku.cn/index.html)，选择TCGA和GTEx数据库中的LUSC样本(486例)和正常样本(338例)对Hub基因进行差异性验证(|Log2FC|>2，p-value<0.05，JitterSize=0.4)。

1.2.6 Hub基因的相关性分析

使用在线网站GEPIA网站，选择TCGA数据库中的LUSC样本(486例)和正常样本(50例)对Hub基因进行差异分析。

2 结果与分析

2.1 COPD患者的DEGs

通过对GSE11784、GSE11906和GSE20257进行RNA水平上正常样本和COPD样本分析，共获得76个 DEGs ，包括 52个上调基因和24个下调基因(见表2)，DEGs的火山图和热图(见图1和图2)。

图1 DEGs的火山图Fig.1 Volcano plot for DEGs注：图中红点代表显著上调基因，绿点代表显著下调基因，黑点代表非显著表达基因，所有差异表达基因符合矫正后p值<0.05.

图2 COPD和正常样本DEGs热图Fig.2 Heat map for DEGs between COPD and normal samples注：图中从红色到绿色，基因表达水平逐渐下降，所有差异表达基因符合矫正后p值<0.05；N代表正常样本，T代表肿瘤样本.

表2 三个GEO数据集中的76个DEGsTable 2 A total of 76 DEGs in three GEO datasets

2.2 GO富集分析

通过R软件的ClusterProfiter和org.Hs.eg.db包进行DEGs的GO 富集分析(见表3)，DEGs主要富集在受体-配体活性、信号受体激活剂活性、四吡咯结合、内肽酶活性、血红素结合和G蛋白偶联受体结合等功能上。

表3 COPD差异表达基因的GO富集分析Table 3 GO enrichment analysis of DEGs in COPD samples

2.3 KEGG代谢通路分析

通过R软件的ClusterProfiter包进行DEGs的KEGG代谢通路分析,结果(见如图3和表4)，DEGs主要富集在细胞色素P450对外源物质的代谢、化学致癌、花生四烯酸代谢和甲状腺激素合成4条代谢途径中。

表4 COPD差异表达基因的KEGG代谢通路富集分析Table 4 KEGG pathway analysis of DEGs in COPD samples

图3 COPD差异表达基因的KEGG代谢通路富集分析Fig.3 KEGG enrichment analysis of DEGs in COPD注：X轴代表差异基因中与该通路相关的基因数与整个差异基因总数的比值，Y轴代表KEGG富集通路.

2.4 COPD差异表达基因蛋白质网络互作分析

使用在线网站STRING (http://string-db.org) 对52个上调基因和24个下调基因进行COPD的蛋白质互作分析，构建PPI网络(见图4)，每个基因代表一个节点，每2个蛋白质之间的连线代表一个关系对，总共有36个节点和51个关系对。

图4 COPD差异基因的PPI网络Fig.4 PPI network of DEGs in COPD samples注：图中圆圈代表基因，连线代表基因间蛋白质的相互作用.

使用Cytoscape软件中的MCODE从PPI网络中筛选得到2个功能模块(见图5)，进一步使用CytoHubba来确定Hub基因(见图6)，评分前十位的基因分别是SPP1、EGF、CCL2、ALDH3A1、BMP4、SPRR3、KRT6A、BPIFB2、CYP1A1和SPRR1B，且都属于上调基因。与图5比较发现，这10个基因有7个存在于两个功能模块中，因此将这10个基因作为Hub基因做进一步分析对于COPD患病机理及分子

图5 用MCODE插件从PPI网络中筛选获得两个关键功能模块(A和B)Fig.5 PPI network of two key modules obtained by MCODE in cytoscape(A and B)

图6 通过PPI网络筛选10个Hub基因Fig.6 Ten hub genes screened by PPI network注：图中颜色越深，基因评分越高.

机制具有一定的研究意义。

2.5 Hub基因在TCGA和GTEx数据库中的差异性验证

TCGA (The Cancer Genome Atlas)是一个涵盖33种癌症及对应正常样本的数据库，GTEx (Genotype Tissue Expression)数据库仅包含正常人的样本，通常与TCGA数据库结合使用，以弥补TCGA正常样本的不足。本研究中，为探寻COPD与LUSC共同的生物标记物，将GEO数据库中获得的COPD的10个Hub基因在TCGA和GTEx数据库的LUSC样本中进行差异性验证。基因SPP1、ALDH3A1、SPRR3、KRT6A和SPRR1B均具有显著性差异(见图7)，表明这5个基因是COPD和LUSC共同的DEGs。进一步分析发现，SPP1和ALDH3A1、SPRR3和KRT6A呈显著正相关(见图8)，其它基因间无显著性关系。

图7 Hub基因SPP1(a)、ALDH3A1(b)、SPRR3(c)、KRT6A(d)、SPRR1B(e)在TCGA和GTEx数据库中LUSC样本中的差异性验证Fig.7 Significance verification of SPP1(a), ALDH3A1(b), SPRR3(c), KRT6A(d), and SPRR1B(e) in LUSC samples in TCGA and GTEx databases

图8 Hub基因相关性分析Fig.8 Correlation analysis of hub genes注：A:SPP1和ALDH3A1相关性分析，B:KRT6A和SPRR3相关性分析.

3 讨论

COPD是一种常见的以持续气流受限为特征的疾病，气流受限进行性发展，与气道和肺脏对有毒颗粒或气体的慢性炎性反应增强有关[13]。小气道管腔狭窄，空气阻力将成倍增加，周围气道阻力增高是小气道病变的重要病理生理学特征，故当小气道有炎症或痰液阻塞，很容易造成闭合、萎陷，COPD病变也多从小气道开始[14]。

本研究以小气道上皮细胞为研究对象，筛选3个COPD的GEO芯片数据集进行生物信息学分析，获得了52个上调基因和24个下调基因。GO分析发现DEGs主要富集在受体-配体活性、信号受体激活剂活性、四吡咯结合、内肽酶活性、血红素结合和G蛋白偶联受体结合等功能上，KEGG代谢通路主要集中在细胞色素P450对外源物质的代谢、化学致癌、花生四烯酸代谢和甲状腺激素合成4条途径上。进一步分析发现参与GO和KEGG代谢通路的DEGs均为上调基因，这些功能和代谢通路均与COPD发病机理有密切关系。程序性细胞死亡受体1(PD-1)通过与其配体PD-L1结合，负调控T淋巴细胞活化及效应功能，参与免疫耐受的维持，从而促进了COPD的发生和发展过程[15]。乙酰胆碱受体M3是引起气道平滑肌收缩的主要受体，激活后发生构象改变，并进一步激活磷脂酶C，催化二磷酸磷脂酰肌醇水解为三磷酸肌醇，促进钙离子通道开放，引起气道平滑肌收缩，气流阻力增大[16]。花生四烯酸是在炎性反应中由细胞膜磷脂在磷脂酶A2催化下产生的，并经过酶促反应生成多种代谢物，它们可以通过多种途径刺激呼吸道中的伤害性感受器，导致炎性细胞聚集而发生局部效应[17]。由此可见筛选得到的差异表达基因在COPD的发病机制中占有重要地位，可为后续研究奠定良好的基础。

进一步取交互作用得分>0.4的蛋白质构建PPI网络，发现共有36个基因参与到网络构建中，为获得关键基因，使用Cytoscape软件中的MCODE和CytoHubba插件对PPI网络图进行分析，分别获得了两个重要的子模块和10个Hub基因，比较发现有7个Hub基因位于子模块中，其中ALDH3A1和CYP1A1为细胞色素P450对外源物质的代谢和化学致癌途径中的重要基因(见表4)，因此可初步认为筛选得到的Hub基因对于进一步研究COPD患病机理及分子机制具有一定的研究意义。为寻求COPD与LUSC共存的DEGs，将10个Hub基因通过TCGA和GTEx数据库验证其在LUSC中的差异表达情况，从而确定SPP1、ALDH3A1、SPRR3、KRT6A和SPRR1B为COPD与LUSC共存的DEGs。进一步分析其相关性发现仅SPP1和ALDH3A1、SPRR3和KRT6A间呈显著正相关(见图8)，目前尚未有该基因对间作用机制的相关报道，是否存在共同的调节因子还有待进一步研究。

早期研究并没有发现SPP1与COPD之间的联系[18]，血液分析结果显示癌症分期和SPP1表达水平相互矛盾[19]，而最新报道表明SPP1在COPD和肺癌患者中均高表达，并且与肺癌较低的生存率密切相关[20]，这与本研究结果一致，可以初步认为COPD中SPP1的高表达与患肺癌风险增加有关，SPP1可能是一个治疗靶点，用于延缓COPD患者肺癌的发展，提高其生存时间。

ALDH3A1属于乙醛脱氢酶家族3中的成员A1，在人体的胃、食道和肺中含量较高。通过前列腺肿瘤异种移植的动物实验发现ALDH3A1可以促进前列腺肿瘤细胞肺转移的发生[21]。ALDH3A1基因表达还与大肝癌早期复发转移关系密切，并且该基因表达越高，复发转移时间越长[22]。目前关于ALDH3A1对COPD和LUSCD之间的关系还鲜有报道，通过KEGG代谢通路分析我们可以看到ALDH3A1参与了细胞色素P450对外源物质的代谢和化学致癌两个途径(见表4)。细胞色素P450是人体内主要的多环芳烃代谢酶，参与内外源性物质的代谢，对于化学致癌物来说，可把无活性的前致癌物激活转变为电子化合物，与 DNA或蛋白质结合，导致癌变。由于ALDH3A1在COPD和LUSC中均上调，可以认为高表达的ALDH3A1加速了细胞色素P450对外源物质尤其是化学致癌物质的代谢，从而促进了COPD和LUSC的发生和发展。

SPRRs为一类富含脯氨酸蛋白的基因家族，该家族共有4个成员：SPRR1、SPRR2、SPRR3和SPRR4。它们蛋白产物的中央片段均由串联重复的氨基酸(XKXPEPXX(X为任意氨基酸))序列组成：SPRR1、SPRR3和SPRR4为8个，SPRR2为9个。本研究涉及到的Hub基因为SPRR1B和SPRR3。研究表明这两个基因参与多种癌症的发生和发展。SPRR1B在口腔鳞癌组织中过表达，从而影响细胞增长[23]，它还参与黑色素瘤的转移[24]，并且是气道上皮细胞鳞状细胞分化的早期生物标志物，其表达下调伴随细胞的恶性转化[25]。通过meta和Kaplan-Meier分析，发现SPRR1B在肺腺癌组织中表达增强，敲除该基因可抑制肺腺癌细胞增殖、迁移和侵袭；此外，通过火山图发现肺癌组织中SPRR3的表达水平也明显增高[26]。同样，SPRR1B和SPRR3在肺癌H1299细胞系中均过表达[27]。本研究与前人结果一致，SPRR1B和SPRR3在COPD和LUSC中均表达上调，提示这两个基因在COPD和LUSC发生机制的探讨中具有潜在的研究价值，建议SPRR家族基因与COPD和肺癌的关系应进一步深入研究。