慢性阻塞性肺疾病(chronic obstructive pulmonary disease,COPD)是全球四大慢性疾病之一,其具有迁延难愈、反复发作和持续进展的特性,2019年,全球报告了2.123亿例COPD流行病例,其中COPD导致330万人死亡
[1-3]。COPD是一项重要的公共卫生挑战,如何从危险因素入手进行早期预防和治疗,成为目前COPD研究的关键。
COPD的危险因素具有多样性,为个体易感因素与环境因素共同作用的结果,其中遗传因素在个体易感因素中占据重要位置。近年来,测序技术的迅速发展和孟德尔随机化(mendelian randomization,MR)方法的广泛应用,为挖掘COPD的遗传相关靶点提供了机会,研究
[4-5]显示:编码基质金属蛋白酶12(matrix metallopeptidase 12,MMP12)和谷胱甘肽S-转移酶(glutathione s-transferase,GST)的基因多态性可能与肺功能的下降有关。全基因扫描研究
[6]显示:α尼古丁乙酰胆碱受体和刺猬因子相互作用蛋白等与COPD或者肺功能有关。国际COPD遗传学联盟最新的研究
[7]结果显示:有82个与COPD有关的基因位点,不同基因与COPD的不同病理或临床特征关联,从遗传基因的角度支持COPD存在异质性。一项MR分析
[8]显示:抗炎性肺蛋白棒细胞分泌蛋白16在血清中的高表达可增加COPD发病的风险并有助于疾病进展,可能成为COPD潜在的生物标志物。但上述研究缺乏大样本数据进行独立验证,本研究利用基因表达综合数据库(Gene Expression Omnibus,GEO)中多个独立的COPD数据集进行联合分析,整合了表达数量性状位点(eQTL-expression quantitative trait locus,eQTL)数据和全基因组关联研究(Genome-Wide Association Study,GWAS)数据,选用MR方法探讨相关基因与COPD间的因果关系,挖掘COPD潜在的治疗靶点,探讨其作用机制,为COPD患者的临床诊疗提供参考。
1 资料与方法
1.1 数据集来源
使用GEO数据库(
https://www.ncbi.nlm.nih.gov/)获取COPD相关基因表达数据集, 以“chronic obstructive pulmonary disease”为检索词检索后获得GSE130928、GSE11906、GSE13896、GSE64614和GSE212331数据集。为了扩大数据来源,将GSE130928、GSE11906、GSE13896和GSE64614O数据集作为训练集进行合并,使用R软件(版本:4.4.0)对上述数据集进行数据读取和校正,使用sva包(版本:3.52.0)进行归一化处理。利用ggplot2(版本:3.5.1)和ggpubr包(版本:0.6.0)进行主成分分析(principal component analysis,PCA),消除批次效用。将GSE212331数据集作为外部验证集对筛选的基因进行验证。见
表1。
COPD的结局数据来源于GWAS汇总数据集(Integrative Epidemiology Unit,IEU)的遗传关联数据库(
https://gwas.mrcieu.ac.uk/)。所使用的GWAS ID为ebi-a-GCST90018807,涉及13 530例COPD病例和454 945例欧洲血统对照者,共包括24 180 654个单核苷酸多态性(single nucleotide polymorphisms,SNPs)。eQTL暴露数据来源于eQTL基因联盟(
https://www.eqtlgen.org/)
[9],包含31 684份全血样本的eQTL数据,用于后续MR分析。
1.2 差异表达基因(differentially expressed genes,DEGs)筛选和可视化
使用R软件中limma包(版本:3.58.1)对归一化后的数据集中的DEGs进行筛选,设置差异倍数(fold change,FC),筛选条件为|log2FC|>0.585,矫正后的P<0.05,筛选后获得DEGs,并绘制火山图和热图,进行可视化展示。
1.3 eQTL数据筛选
获得的eQTL数据作为暴露因素,使用R软件中TwoSampleMR包(版本:0.6.2)筛选出强相关的SNPs(P<5×10-8)作为工具变量(instrumental variable,IV),为了去除连锁不平衡,设置r2 <0.001,遗传距离为10 000 kb,同时设置F>10排除弱关联或表型方差解释不足的SNP。
1.4 MR分析
本次研究中首先将eQTL数据作为暴露因素、SNPs作为工具变量、GWAS疾病数据集(ebi-a-GCST90018807)作为结局数据进行MR分析。采用R软件中TwoSampleMR”包和“forestplot”包进行数据分析及可视化。本次分析采用逆方差加权法(inverse-variance weighted,IVW)为主分析方法,MR-Egger法和加权中位数方法进行数据分析。为了更好地确定相关基因,设定该基因至少应满足IVW的结果有显著性差异(
P<0.05),同时加权中位数和MR-Egger的结果与IVW法随机效应模型的结果必须方向相同。MR分析后,采用Cochran’s Q检验评估IVW 随机效应模型异质性,若检验的结果显示
P<0.05,则表示存在异质性,将该基因剔除
[8]。采用MR-Egger回归确定潜在的水平多效性,结果中截距项的值接近于0且
P>0.05时,则该MR结果不存在多效性
[9],最后,再通过漏斗图对结果进行评判,漏斗图对称则表示结果无潜在的异质性
[10]。
对“1.2”中得到的DEGs与结果进行交集,获得与疾病相关的交集基因,同时根据其在试验组和对照组的表达量差异,分为上调基因与下调基因,再对得到的交集基因作为暴露因素再次进行MR分析,评估其与疾病的因果关系。
1.5 交集基因的基因本体论(Gene Ontology,GO)富集分析和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)信号通路富集分析
利用R软件中clusterProfiler包(版本:4.10.1)对共表达基因进行GO生物功能富集分析和KEGG信号通路富集分析,设置过滤条件为P<0.05,对富集显著的生物功能或通路以柱状图进行表示。为了明确交集基因的染色体位置,使用OmicCircos(版本:1.32.0)对交集基因的基因组位置进行可视化。
1.6 基因集富集分析(Gene Set Enrichment Analysis,GESA)
为了进一步挖掘上述基因的功能和通路富集情况,本课题组利用从分子特征数据库(
https://www.gsea-msigdb.org/gsea/msigdb)下载所得的c2.cp.kegg.Hs.symbols.gmt基因集为参考,对交集基因进行分析,以
P<0.05为筛选条件确定显著富集的通路,对其中归一化富集评分(normalized enrichment score,NSE)>0和NSE<0的前5个通路进行可视化。
1.7 外部验证组差异分析
使用limma包(版本:3.58.1)读取GSE212331数据集中试验组与对照组的表达数据,以验证共表达基因是否在对照组与试验组之间表现出差异,将上述结果与MR分析结果进行比较,并进行可视化。
2 结 果
2.1 DEGs鉴定结果
对GSE130928、GSE11906、GSE13896和GSE64614数据集中的每个基因的表达值进行校正及合并,并通过PCA分析消除批处效应。经过校正后,数据集中的所有样本都达到了可接受的同质性(
图1)。通过对校正后的合并数据集中对照组和试验组的表达数据进行分析,共获得1 571个DEGs,其中820基因表达上调,751个基因表达下调(
图2A)。按表达量进行排序,前10个表达上调的基因:细胞因子样蛋白1(cytokine-like 1,
CYTL1)、安眠蛋白1α(meprin alpha,
MEP1A)、泛素羧基端酯酶L1(ubiquitin C-terminal hydrolase-L1,
UCHL1)、髓系C型凝集素结构域家族5成员A(myeloid C-type lectin domain family 5 member A,
CLEC5A)、脱氧核糖核酸酶2B(deoxyribonuclease 2 beta,
DNASE2B)、基质金属蛋白酶12(matrix metalloproteinase 12,
MMP12)、趋化因子2(chemokine 2,
CCL2)、几丁质酶1(recombinant chitinase 1,
CHIT1)、细胞色素P450家族成员1B1(cytochrome P450 1B1,
CYP1B1)和细胞色素P450家族成员1A1(cytochrome P450 1A1,
CYP1A1);前10个表达下调的基因:分泌球蛋白家族3A成员1(secretoglobin family 3A member 1,
SCGB3A1)、分泌珠蛋白家族1A成员1(secretoglobin family 3A member 1,
SCGB1A1)、维生素K依赖性蛋白S(vitamin K-dependent protein S,
PROS1)、 角 蛋 白 19(keratin 19,
KRT19)、人类白细胞抗原B(human leukocyte antigen-B,
HLA-
B)、半胱氨酸丰富蛋白1(cysteine rich protein 1,
CRIP1)、视黄酸受体应答蛋白3(retinoic acid receptor response protein 3,RARRES3)、水通道蛋白3(aquaporin 3,
AQP3)、四旋蛋白6(tetraspanin 6,
TSPAN6)和溶质载体家族 22 成员 4(solute carrier family 22 member 4,
SLC22A4)(
图2B)。
2.2 MR分析结果
经过筛选,最终获得25 457个强关联SNPs作为工具变量,所有SNPs的
F值均>10。同时鉴定出了286个COPD相关基因,通过与“2.1”中获得的DEGs取交集得出3个上调基因,包括二酰甘油激酶γ(diacylglycerol kinase gamma,
DGKG)、重肽神经丝蛋白(neurofilament heavy polypeptide,
NEFH)和Fc受体样B(Fc receptor like B,
FCRLB);8个下调基因,包括HCP5(HLA complex P5,
HCP5)、PTGFR(prostaglandin F receptor,
PTGFR)、金属还原酶4(
STEAP4 metalloreductase,
STEAP4)、含普列克底物蛋白家族F成员2(pleckstrin homology domain containing family F member 2,
PLEKHF2)、分化簇3D(CD3d molecule,
CD3D)、转胶蛋白2(transgelin 2,
TAGLN2)、三基序蛋白22(tripartite motif containing 22,
TRIM22)和核糖体蛋白L9(ribosomal protein L9,
RPL9)(
图3)。将上述11个交集基因作为暴露因素分别进行MR分析,评估上述基因与COPD的因果关系,结果显示:
HCP5与
PTGFR截距项的值接近于0且
P>0.05,存在多效性,将其剔除。最终结果显示:在采用IVW方法进行的MR分析中,3个上调的基因与COPD存在显著的正相关关系,
DGKG[比值比(odds ratio,OR)=1.184,95%置信区间(confidence interval,CI):1.019~1.376,
P=0.027];
NEFH (OR=1.078,95%CI:1.011~1.148,
P=0.021);
FCRLB(OR=1.085,95%CI:1.004~1.172,
P=0.04)。6个下调的基因降低COPD发生风险,
STEAP4(OR=0.963,95%CI:0.934~0.993,
P=0.017);
PLEKHF2(OR=0.942,95%CI:0.892~0.995,
P=0.032);
CD3D(OR=0.855, 95%CI: 0.773~0.946,
P=0.002);
TAGLN2(OR=0.955, 95%CI: 0.912~1.000,
P=0.049);
TRIM22 (OR=0.921, 95%CI:0.855~0.992,
P=0.03);
RPL9(OR=0.888,95%CI:0.823~0.958,
P=0.002),见
图4。表明3个上调基因增加COPD发生风险(OR>1),6个下调基因降低COPD发生风险(OR<1)。交集基因在染色体中分布图见
图5。
2.3 GO功能富集分析和KEGG信号通路富集分析结果
GO功能富集分析共得到64条结果,其中生物过程(biological process,BP)43条,主要涉及:铁输入细胞、蛋白质同源三聚化和外周神经系统神经元的发育及分化等;细胞组分(cellular component,CC)11条,主要涉及α-β T细胞受体复合体、突触后细胞骨架和早期内质体膜等;分子功能(molecular function,MF)10条,主要涉及:ATP 依赖性二酰甘油激酶活性和氧化还原酶活性等。KEGG通路富集分析富集到11条通路,主要涉及原发性免疫缺陷症和甘油酯代谢以及辅助性T(helper T,Th)1和Th2细胞分化等(
图6)。
2.4 GESA
MR分析中,上调基因中的
NEFH基因表达量升高增加了COPD发生风险(
P=0.021),下调基因中
CD3D和
RPL9基因表达量升高降低了COPD发生风险(
P=0.002),为了进一步挖掘上述3个基因的功能,进行了GESA,结果显示:
NEFH基因高表达组中前5个通路是细胞因子-细胞因子受体相互作用、成熟期发病的青年糖尿病、神经活性配体受体相互作用、嗅觉传导和味觉传导;
NEFH基因低表达组中前5个通路是亨廷顿氏症、氧化磷酸化、帕金森病、核糖体和鳞状细胞;
CD3D基因高表达组中前5个通路是细胞因子-细胞因子受体相互作用、Janus激酶(Janus kinase,JAK)-信号转导和转录激活因子 (signal transducer and activator of transcription,STAT)信号通路、成熟期发病的青年糖尿病、神经活性配体受体相互作用和嗅觉传导;
CD3D基因低表达组中前5个通路是亨廷顿症、溶酶体、帕金森病、核糖体和鳞状细胞;
RPL9基因高表达组中前5个通路是亨廷顿氏症、氧化磷酸化、帕金森病、核糖体和鳞状细胞;
RPL9基因低表达组中前5个通路是细胞因子-细胞因子受体相互作用、成熟期发病的青年糖尿病、神经活性配体受体相互作用、嗅觉传导和味觉传导(
图7)。上述结果中
NEFH基因为上调基因,
RPL9基因为下调基因,所富集到的通路在高表达组与低表达组中结果相反,证实本研究结果的可靠性。
2.5 外部验证组差异分析
COPD样本组中
DGKG、
NEFH和
FCRLB基因表达水平均高于健康对照组,其中FCRLB表达量升高最明显(
P<0.01)。与健康对照组比较,COPD样本组中的
STEAP4、
PLEKHF2、
CD3D、
TRIM22和
RPL9这5个基因表达量下调,其中
CD3D和
RPL9基因表达量下调最显著(
P<0.05或
P<0.01),以上结果与MR分析结果基本一致,证明了本研究结果的可靠性。见
图8。
3 讨 论
COPD是呼吸系统常见病,具有遗传易感性,目前相关研究一直在致力于COPD早期诊断或疾病前诊断,因此利用MR方法筛选合适的诊断标志物或治疗靶点,是目前较为合适且可靠的方法
[11-13]。COPD的药物治疗主要包括支气管扩张剂、吸入性糖皮质激素、磷酸二酯酶4(phosphodiesterase 4,PDE‐4)抑制剂和祛痰药及抗氧化剂等
[1,14],而相关生物制剂或靶向药物的使用也为COPD治疗提供了新选择
[15],因此,研发相关生物制剂或靶向药物,可能是COPD未来治疗的主要发展方向
[16]。
本研究利用GEO数据库中的多个数据集和GWAS数据库中COPD遗传数据及eQTL数据进行MR分析,鉴定出9个共表达基因,包括3个上调基因:
DGKG、
NEFH和
FCRLB基因,其表达量升高增加了COPD的发生风险;6个下调基因:
STEAP4、
PLEKHF2、
CD3D、
TAGLN2、
TRIM22和
RPL9基因,其表达量减少也与疾病风险的降低有关。上述基因的表达模式为COPD的遗传基础提供了新的见解,可能反映了疾病病理中的关键分子变化,如缺氧诱导的肿瘤血管内皮细胞中DGKG过度表达可通过E盒结合锌指蛋白2(zinc finger E-box-binding protein, ZEB2)/转化生长因子β1(transforming growth factor-β
1,TGF-β
1)轴促进肿瘤血管生成和免疫抑制调节性T淋巴细胞分化,与COPD中的气道重塑与免疫调节机制基本一致
[17-18]。SONG等
[19]利用机器学习算法以
FCRLB等7个核心基因构建肺腺癌焦亡相关的预后风险评分模型,具有较好的准确度和灵敏度,而细胞焦亡在COPD的发生发展中起关键作用
[20]。WANG等
[21]发现:
CD3D基因与儿童哮喘密切相关,并在哮喘患者的肺泡灌洗液细胞中进行了验证。免疫细胞浸润分析
[22]结果显示:
CD3D基因与静息肥大细胞和嗜酸性粒细胞增加呈负相关关系,与Th1、Th2和Th17细胞中的几种细胞标志物高度相关,与目前COPD的免疫炎症反应机制结果一致。RPL9是一种新型晚期糖基化终末产物结合蛋白,可以减少脂多糖和高迁移率族蛋白B1(high mobility group box 1,HMGB1)刺激下的巨噬细胞中TNF-α的表达,从而减轻炎症
[23-24]。另一项针对运动机能损伤的研究
[25]显示:COPD患者股四头肌中微小RNA542-3p/5p(microRNA 542-3p/5p,miR-542-3p/5p)水平升高,抑制核糖体蛋白合成,并与其肺部疾病的严重程度成正比,该结果表明核糖体蛋白在COPD进程中可能发挥着重要作用。目前,关于上述基因在COPD发病机制中的研究较少,但本文作者发现多数蛋白与COPD发病机制密切相关,如干预慢性炎症、减轻细胞焦亡以及免疫调节等方面,但上述机制仍需进一步深入探讨和验证。
本研究中,GO功能富集分析结果显示:交集基因主要涉及铁输入细胞、α-β T细胞受体复合体和氧化还原酶活性等;KEGG信号通路富集分析结果显示:共表达基因主要涉及甘油酯代谢以及Th1和Th2细胞分化等。上述结果与既往研究
[26-29]一致,NAJAFINOBAR等
[26]利用飞行时间二次离子质谱技术(time-of-flight secondary ion mass spectrometry,ToF-SIMS)法对COPD患者肺组织中的铁积累进行量化,结果显示:与对照组比较,COPD患者肺组织切片中铁含量升高12倍。另一项研究
[27]也证实:气道铁含量升高会导致更频繁的COPD急性加重。COPD是一种异质性慢性炎症性疾病,Th1和Th2细胞分化在COPD的发病机制中扮演了重要角色,其产生的多种细胞因子或介质,可直接介导COPD炎症的发生
[28-29]。基于上述研究,本文作者推测:交集基因可能是通过干预Th1和Th2细胞分化影响自身免疫,以及抑制铁离子流入细胞来减轻COPD慢性炎症,影响疾病进程。
为了进一步探讨共表达基因在染色体中的分布,本文作者确定了基因位点,发现共表达基因与多条染色体相关。相关研究
[30]显示:X染色体与COPD表型有关,
TMSB4X基因附近的1个变体rs5979771与肺功能相关并呈现全基因组显著性。本研究结果显示:位于1号和11号染色体上的基因最多,该结果与既往的GWAS分析
[31-32]具有一致性。上述染色体可能与COPD显著相关,对COPD相关SNPs的鉴定也具有一定的意义。
基于MR分析结果,本文作者选择了3个显著关联的
NEFH、
CD3D和
RPL9基因,利用GESA富集分析探讨其潜在影响通路,同时基于外部验证数据集对交集基因进行验证,结果显示:
CD3D和
RPL9基因在COPD样本组中下调;GESA富集分析结果显示:3个基因主要富集于细胞因子-细胞因子受体相互作用、氧化磷酸化、核糖体和鳞状细胞等通路,上述通路均与COPD的发生发展密切相关
[33-35]。YEW-BOOTH等
[33]利用蛋白印迹或免疫组织化学技术对吸烟者、COPD患者和健康人肺组织中JAK-STAT信号通路蛋白进行检测,结果表明:JAK-STAT信号通路在COPD患者中被异常激活,证实了JAK-STAT信号通路在COPD发生发展中的作用。另一项动物实验
[36-37]结果显示:过氧化物酶体增殖物激活受体-γ激动剂可以通过抑制JAK-STAT信号通路明显改善香烟烟雾或脂多糖刺激的巨噬细胞的极化障碍和功能活性,减轻气道炎症,恢复肺功能。研究
[38]显示:受环境因素等影响,COPD患者外源性活性氧和肺部炎症及结构细胞产生的内源性活性氧增多,造成线粒体功能障碍,引起气道上皮细胞损伤,诱导平滑肌细胞增生,导致肺功能下降,与氧化还原及氧化磷酸化反应密切相关。目前,临床中抗氧化剂的应用是COPD临床治疗的关键措施,如N-乙酰半胱氨酸等,可以改善并延缓COPD的急性发作。上述研究表明:
NEFH、
CD3D和
RPL9基因可能通过JAK-STAT信号通路影响COPD进程中的氧化应激及线粒体功能,进一步影响COPD发生发展,该推论可能是COPD的潜在机制,为后续治疗及机制探讨提供了参考。
本研究基于MR和GEO数据库挖掘出在COPD疾病过程中与遗传因素高度相关的9个基因,通过后续的深入分析基因的MF和BP等,明确出NEFH、CD3D和RPL9基因为3个核心基因,较完整地揭示了基因与疾病的相关性以及未来向临床应用转化的可能性,但仍需要后期进一步实验验证。本研究中相关步骤虽然较为完整,但仍有分析技术存在局限性,相关数据集选择存在主观性等不足。但该研究筛选出的基因,可为COPD精准医疗的发展提供新参考,为后续研究和相关药物的研发提供了一定的方向,具有一定的临床意义。