心肌梗死潜在生物标志物的生物信息学研究

钱 晗, 曹梦菲, 吕书梅, 袁 伟

(江苏大学附属医院 心血管内科, 江苏 镇江, 212001)

在人口老龄化和代谢危险因素的双重压力下，中国居民的心血管疾病风险持续增加，心肌梗死患者的病死率亦逐渐上升[1]。心肌梗死是指心肌缺血性坏死，通常是在冠状动脉病变的基础上发生冠状动脉血供减少或中断，心肌细胞因严重且持续的缺血缺氧而发生坏死[2]。心肌细胞的死亡会引发一系列免疫炎症反应，这些免疫炎症反应介导的受损心肌修复在梗死后心室重构、心功能和患者结局等方面起着关键作用[3]。急性心肌梗死的诊断主要依据临床表现、血清生物学标志物和心电图，其中临床常用的生物学标志物包括肌红蛋白、肌酸激酶-MB、心肌肌钙蛋白I和心肌肌钙蛋白T[4-5]。然而，这些标志物水平在心力衰竭、肾衰竭和甲状腺疾病患者中也同样升高，故结合多组学平台探寻更多特异性的生物学标志物非常必要[6-7]。外周血单个核细胞(PBMCs)包括T细胞、B细胞、单核细胞、树突状细胞和自然杀伤细胞，在免疫应答、免疫监测和免疫治疗中发挥着重要作用[8]。本研究基于PBMCs分析心肌梗死和梗死恢复过程中潜在的关键基因和生物学改变，旨在为心肌梗死的诊断、评估、管理和预后分析提供新的见解。

1 材料和方法

1.1 数据集获取

从GEO数据库中下载关于心肌梗死PBMCs的转录组测序芯片数据集GSE59867[9-10], 芯片注释平台为GPL6244, 芯片信息为Affymetrix Human Gene 1.0 ST Array。该数据集收集了111例ST段抬高型心肌梗死(STEMI)患者入院、出院、心肌梗死后1个月、心肌梗死后6个月时的外周血样本数据。对照组为46例无心肌梗死病史的稳定型冠心病患者。基于芯片注释平台对数据集进行基因名转换，用于后续共表达网络的构建。

1.2 加权基因共表达网络分析(WGCNA)

通过WGCNA R包筛选与心肌梗死有关的易感基因模块[11]。WGCNA是一种常用的生物信息学方法，可以识别具有相似表达模式的基因模块，分析基因模块与样本表型之间的关系，绘制基因模块中的调控网络，识别关键模块和调控基因[12-13]。首先，选取数据集中方差变异程度最高的前3 000个基因构建共表达网络; 其次，利用R函数pickSoftThreshold计算软阈值1～20; 然后，对这些基因进行层次聚类和动态树切割来划分模块，并合并相似的模块(阈值为0.25); 最后，计算模块与疾病之间的Spearman相关系数，进而确定易感基因模块。为了探讨所选模块的潜在生物学功能，通过ClusterProfiler R包对模块中的基因进行基因本体论(GO)分析和京都基因与基因组百科全书(KEGG)通路分析[14]。GO分析包括生物过程(BP)、分子功能(MF)和细胞成分(CC)。通过Ggplot2 R包进行结果展示。

1.3 基因差异分析和转录因子预测

使用Limma R包进行基因差异分析[15]。差异表达基因(DEGs)的筛选标准为矫正后P<0.05, |log2(倍数变化)|>0.5。基于ClusterProfiler R包对DEGs进行GO分析和KEGG通路分析,P<0.05表示显著富集。对上述易感基因模块和DEGs进行交集分析，通过ChEA3数据库对交集基因进行转录因子预测。ChEA3是一种转录因子富集分析工具，整合了多个数据库，如ENCODE、GTEx、ARCHS4等数据库[16]。由于ChEA3数据库并未对预测结果可靠度进行严格划分，参考相关研究[16-17]结论，本研究对平均排名分≤30分的共同转录因子进行结果展示。基于STRING数据库构建交集基因的蛋白相互作用网络，并将交互评分>0.4分的数据导入Cytoscape 3.8.2软件进行可视化处理[18-19]。

1.4 关键基因的确定和受试者工作特征(ROC)曲线分析

LASSO回归分析是常用的筛选变量的压缩估计方法，通过构造惩罚函数获得较为精炼的模型，使得其压缩一些回归系数，同时设定一些回归系数为0，因此保留了子集收缩的优点，是一种处理具有复共线性数据的有偏估计，常被用于协变量筛选[20]。基于glmnet R包对交集基因中存在蛋白相互作用的基因进行LASSO回归分析，进一步筛选关键基因[21]。采用10倍交叉验证法选择惩罚项(λ), 并选择最小二项偏差在1个标准误差内的最简单模型的λ值(lambda.1se)。通过pROC R包对关键基因进行ROC曲线分析，计算曲线下面积(AUC), 评估关键基因对心肌梗死的诊断价值[22]。

1.5 基因集富集分析(GSEA)和外部数据集验证

根据每个关键基因的表达量将数据集划分为高表达组和低表达组，通过GSEA探讨高表达组和低表达组之间潜在的生物学功能变化，将特征基因集作为预先定义的基因集。作为一种常见的生物信息学方法, GSEA可评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的影响[23]。从GEO数据库中选取另一个数据集GSE123342, 将其用于验证关键基因的表达。GSE123342数据集包括急性心肌梗死(65例)、梗死后30 d(64例)、梗死后1年(37例)、稳定型冠心病(22例)样本数据以及4例技术重复样本数据。

2 结 果

2.1 易感基因模块的确定

根据无标度网络原理，选取10作为共表达网络的软阈值，见图1A。通过层次聚类和动态树切割，将前3 000个基因划分为15个基因模块，见图1B。剔除聚类失败的灰色模块，棕色模块的总体表达与心肌梗死的相关性最高(r=0.4), 且与出院至心肌梗死后6个月的演变过程呈负相关(r=-0.33), 见图1C。棕色模块中的共表达基因伴随着心肌梗死的发生呈现高表达趋势，而在心肌梗死恢复过程中则呈现低表达趋势。由此提示，棕色模块可能是易感基因模块，其中包含243个基因，即ABHD5、ACER3、ACP3、ACSL1、ADAM17、ADAM9、ADM、AHR、ALCAM、ALDH2、ALPK1、ANKRD50、ANO5、AQP9、ARHGAP24、ARHGAP29、ARRDC4、ASPH、ATP6V0A1、ATP6V1A、B3GNT5、BACH1、BLVRB、BST1、C3AR1、C5、C9orf72、CALCRL、CAPG、CAPZA2、CARD6、CASP1、CASP5、CCDC88A、CCR2、CD14、CD163、CD1D、CD33、CD36、CD63、CLEC1A、CLEC4A、CLEC4D、CLEC4E、CLEC7A、CLMN、CNTLN、CPED1、CPM、CPNE8、CPVL、CR1、CR1L、CREG1、CRISPLD2、CSTA、CTNNA1、CTSH、CYBRD1、CYP1B1、CYP1B1-AS1、DACH1、DDIAS、DOCK4、DOCK5、DRAM1、DSC2、DSE、DUSP6、EDNRB、ENTPD1、ERLIN1、EVI5、F5、FAM114A1、FAM13A、FAM151B、FAM20A、FAM198B、FAR2、FBN2、FBP1、FCGR2A、FGD4、FGD6、FLT3、FLVCR2、FMO5、FPR1、FPR2、FRRS1、FUCA1、FUCA2、GAPT、GAS2L3、GCA、GIMAP8、GLT1D1、GPAT3、GPR141、HAL、HGF、HMGB2、HNMT、HORMAD1、HP、HPSE、IDH1、IFNGR1、IL13RA1、IL15、IL18、IMPA2、IRAK3、JAK2、KCNJ15、KCNJ2、KIF13A、KLHL8、KYNU、LHFPL2、LILRA5、LIN7A、LIPN、LMNB1、LPCAT2、LRMDA、LRRK2、LTA4H、LY86、LYVE1、MAP3K20、MCEMP1、ME1、METTL7A、METTL7B、MFSD1、MGST1、MGST2、MILR1、MNDA、MOSPD2、MS4A4A、MSRB1、MTARC1、MYOF、NAAA、NAIP、NETO2、NLN、NLRC4、NPL、OLFML2B、P2RX7、P2RY13、PGD、PIP4P2、PLA2G4A、PLBD1、PLD1、PLIN2、PLSCR1、PPARG、PPT1、PRRG4、PSTPIP2、PYGL、QPCT、RAB39A、RALB、RBP7、RGL1、RNASE2、RNASE6、RNF141、RNF217、RRAGD、S100A12、S100A9、S100Z、SAT2、SEMA3C、SEPTIN10、SERPINB1、SERPINB2、SESTD1、SGMS2、SH3PXD2B、SHTN1、SIPA1L2、SIRPD、SLC15A2、SLC1A3、SLC22A15、SLC22A4、SLC26A8、SLC31A1、SLC36A4、SLC7A7、SLC8A1、SLITRK4、SMAD1、SMPDL3A、SNX10、SOD2、SORT1、SPATA6、ST3GAL6、ST6GALNAC3、STEAP4、SULT1B1、TASL、TBC1D12、TCN2、TDRD9、TFEC、TGFBI、TLR1、TLR2、TLR4、TLR5、TLR6、TLR7、TLR8、TM6SF1、TMEM144、TMEM167A、TMTC2、TNFAIP6、TNFSF13B、TPST1、TSPO、UBE2D1、UGGT2、VNN2、VNN3、WASF1、WDFY3、WLS、ZC3H12C、ZFYVE16、ZNF438。

棕色模块GO分析最显著的结果分别是骨髓白细胞激活(BP)、囊泡(CC)和模式识别受体活性(MF), 见图1D; KEGG分析结果则主要为感染、免疫炎症等相关生物学过程的改变，包括Toll样受体信号通路等，见图1E。

2.2 基因差异分析

进一步对数据集进行基因差异分析，根据筛选标准共得到142个DEGs(图2A), 其中上调的DEGs有77个，下调的DEGs有65个，见表1。BP分析结果主要包括免疫反应、免疫系统过程和免疫效应过程, CC分析结果包括分泌颗粒、分泌小泡和细胞质小泡部分等, MF分析结果主要包括信号受体活性、分子传感器活性和碳水化合物结合，见图2B。KEGG分析提示，这些DEGs主要涉及免疫、炎症和感染的相关生物学改变，见图2C。DEGs与上述棕色模块存在35个交集基因，见图2D。基于ChEA3数据库对35个交集基因进行共同转录因子预测，其中平均排名分≤30分的共同转录因子分别是CREB5、MTF1、NFE4、SPI1、ZNF467、NFE2、TFEC、MXD1、NR1H3、BORCS8MEF2B、NFIL3和CEBPE, 见表2。构建蛋白互作网络，其中有18个基因存在蛋白质相互作用，见图2E。

表1 差异基因与交集基因

表2 共同转录因子预测结果(平均排名分≤30分)

2.3 关键基因的确定

对18个存在蛋白相互作用的基因进一步行LASSO回归分析，这些基因不同惩罚参数值所对应的系数见图3A。本研究选取最小二项偏差在1个标准误差内的λ值(lambda.1se, 6个基因)，该λ值提供更精简的模型，见图3B。ROC曲线分析显示,CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B基因表达水平对心肌梗死的发生均具有良好的诊断价值,AUC分别为0.831、0.798、0.775、0763、0.866和0.829, 见图3C。与对照组稳定型冠心病患者相比，这些关键基因的表达水平在入院时心肌梗死患者中显著上升，并在出院时、心肌梗死后1个月、心肌梗死后6个月逐步下降，见图3D。

2.4 GSEA结果和验证

GSEA结果提示，这些关键基因涉及的生物学改变主要与糖脂代谢、活性氧、免疫炎症等有关，见图4。这些基因在外部数据集中大多也存在差异表达。基于外部验证数据集GSE123342, 本研究同样发现，与稳定型冠心病相比,HP、FAM198B、CD163、FAM20A、MCEMP1表达水平在急性心肌梗死发生时显著上升，而RNASE2仅表现出上升趋势，见图5。

3 讨 论

本研究基于PBMCs转录组学的变化探讨心肌梗死及梗死恢复过程中可能存在的生物学改变和潜在标志物, WGCNA和DEGs分析结果显示，心肌梗死过程伴随着免疫炎症紊乱。基于LASSO分析，本研究从上述2种分析方法得到的共同基因中鉴定出6个关键基因，即CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B。ROC曲线分析结果提示，这些关键基因对心肌梗死的发生均具有较高的诊断价值。本研究还发现，在梗死恢复过程中，这些基因的表达呈现明显下降趋势。GSEA分析结果表明，这些关键基因涉及的生物学改变主要与糖脂代谢、活性氧、免疫炎症等有关。

CD163分子是Ⅰ型膜蛋白，蛋白表达限于单核细胞/巨噬细胞系, CD163抗原特异性释放机制可能在炎症调节过程中起重要作用。在动脉粥样硬化过程中, CD163+巨噬细胞能促进血管生成和增加血管通透性，并伴随炎症反应；此外，破裂的冠状动脉斑块中CD163的表达增加了心肌梗死和冠心病发生风险[24]。RNASE2编码的蛋白质是非分泌型核糖核酸酶，属于胰核糖核酸酶家族[25]。FAM20A是激酶编码基因家族的成员，本身不具有激酶活性，通过与家族成员FAM20C形成复合物，增强FAM20C的激酶活性，从而使分泌通路内的蛋白磷酸化[26]。MCEMP1又称C19ORF59, 可编码表达于肥大细胞、巨噬细胞等的跨膜蛋白[27]。研究[27]表明,MCEMP1基因可能是卒中诊断和预后评估的新生物标志物。生物信息学研究[28-29]提示,RNASE2、FAM20A和MCEMP1可能是心肌梗死的关键基因，然而这些基因在心肌梗死及梗死恢复过程中的具体作用机制有待进一步研究。HP基因编码一种前蛋白，经处理后产生α链和β链，结合为四聚体产生触珠蛋白，其基因型与急性心肌梗死的发病风险密切相关，并且能决定心肌梗死面积[30-31]。FAM198B(GASK1B)是目前未知功能的新基因，可能编码定位在高尔基体上的膜结合糖蛋白，被认为参与癌症的转移与进展[32]。本研究分析结果显示,FAM198B可能是心肌梗死和反映心肌梗死恢复的关键基因。

本研究尚存在一定局限性: 首先，本研究基于GEO数据库中下载的心肌梗死患者PBMCs的测序数据集进行分析，还需要进一步开展分子生物学实验来验证这些基因在心肌梗死中的表达及其可能的生物学机制; 其次，本研究仅分析了现有的样本数据，未来还需基于更大的样本规模和更详细的样本类型进一步深入研究; 最后，由于缺乏临床资料，本研究未分析这些关键基因的表达与心肌梗死患者年龄、性别、射血分数等临床指标的相关性。

综上所述，本研究基于PBMCs确定了6个可能与心肌梗死密切相关的关键基因，即CD163、RNASE2、HP、FAM20A、MCEMP1和FAM198B, 为心肌梗死的诊断、评估、管理和预后分析提供了新的思路。