张飞飞,袁 红,谢映春,韩润川

1.四川省医学科学院·四川省人民医院输血科,四川成都 610072;2.陕西省中医医院检验科,陕西西安 710003

骨髓增生异常综合征(MDS)是一种异质性的髓系克隆性疾病,此类疾病患者的外周血细胞减少,且转化为急性髓细胞性白血病(AML)的风险增加,此类疾病好发于老年男性[1-2]。随着测序技术的发展,已有研究证明MDS患者的基因位点突变(如DNA甲基化、染色质修饰、RNA剪接、转录、信号转导等)与正常造血的破坏,MDS的形成,MDS转化为AML等过程密切相关[1],并已有部分突变基因被推荐为MDS的鉴别诊断和危险度分层的检测项目,如TP53、TET2、DNMT3A等[2-4]。得益于高通量测序的发展,研究者可从各种生物数据库挖掘疾病相关的基因,分析其基因表达情况,为深入研究提供一定的思路。本研究通过基因表达综合数据库(GEO)及生物信息学相关软件,挖掘MDS的疾病特征基因,旨在为探索MDS的发病机制、诊断、治疗及预后判断提供基础。

1 资料与方法

1.1资料搜集 利用国家生物技术信息中心(NCBI)中GEO数据库(https://www.ncbi.nlm.nih.gov/geo/),以“myelodysplastic syndrome”、“homo”为关键词搜索目标芯片。纳入标准:(1)所选数据集为全基因组的mRNA芯片数据;(2)数据集中有疾病组和健康对照组,均不少于10个样本。经筛选后,采用以GPL570为平台的数据集GSE4619、GSE19429、GSE58831进行数据挖掘分析,随机将3组独立数据按2∶1的比例分为训练集和测试集,即GSE4619与GSE19429合并后的校正数据作为训练集,GSE58831作为测试集。在训练集中包含28个正常样本(健康对照组)、238个MDS样本(疾病组);测试集中包含11个正常样本、159个MDS样本。

1.2方法

1.2.1MDS的差异表达基因(DEGs)筛选 应用R语言的Limma包对训练集中的疾病组及健康对照组进行DEGs的筛选(校正后P<0.05、|logFC|>1),并绘制DEGs的火山图。

1.2.2DEGs的基因本体(GO)功能富集及京都基因与基因组百科全书(KEGG)调控通路富集 使用DAVID2021(https://david.ncifcrf.gov/)在线分析工具对1.2筛选出来的DEGs进行GO功能富集和KEGG信号通路分析,并通过ggplot2包进行可视化。

1.2.3MDS的加权基因共表达网络(WGCNA)分析 应用R语言的WGCNA包对训练集的疾病组与健康对照组的表达数据进行共表达网络分析。该方法用来描述不同样本之间的基因关联模式,是对训练集中表达数据的另外一种分析方法。该方法首先通过对基因相关系数取n次幂的方式计算任意两个基因之间的相关系数,并通过合适的power值构建最佳的无尺度网络,再将邻接网络转化为拓扑重叠(TOM)来计算基因之间的关系;另外,再基于TOM值的相异度对基因构建层次聚类树,筛选出连接度高的基因并定义为模块;再者,通过模块与疾病状态的相关系数及P值绘制模块与疾病的相关性热图;基于基因与模块之间的相关性(MM值)、基因的重要性(GS值)得到每个模块中的基因重要性图形和每个模块的散点图;最后,根据基因重要性的过滤条件(geneSigFilter=0.5)和基因与模块相关性的过滤条件(moduleSigFilter=0.8)筛选出每个模块的核心基因(即Hub基因)[5-6]。

1.2.4Hub基因与DEGs取交集 利用R语言的venn包对DEGs和模块Hub基因取交集,得到交集差异基因。

1.2.5交集差异基因的Lasso回归分析 该方法是一种压缩估计,对交集差异基因进行进一步筛选。它通过构造一个惩罚函数得到一个较为精炼的模型,使得它压缩一些回归系数,即强制系数绝对值之和小于某个固定值;同时设定一些回归系数为零。用glmnet包对疾病的特征基因数目进行筛选。

1.2.6训练集和测试集中疾病特征基因的表达量 利用limma包与ggpubr包分别将训练集、测试集的表达数据与Lasso回归分析筛选出来的特征基因进行分析,以MDS患者为疾病组,健康者为健康对照组,对每个基因进行循环并统计分析,最后绘制差异箱线图(P<0.05)。

1.2.7单样本基因富集(ssGSEA)分析 该方法主要用来量化免疫浸润。运用R的扩展包(GSEABase包、GSVA包等)及免疫基因集文件对训练集中每个样本的数据进行免疫基因打分;利用vioplot包对免疫细胞的打分进行统计分析并可视化(P<0.05);最后对特征基因进行免疫相关性统计分析并绘制热图(P<0.05)。

1.3统计学处理 采用R4.1.2进行统计分析。

2 结 果

2.1MDS相关DEGs的筛选及分析 本研究将GSE4619与GSE19429合并后的校正数据作为训练集,绘制疾病组与健康对照组DEGs的火山图,见图1。从中获得88个DEGs,其中上调基因11个,下调基因77个。

注:|logFC|>1,校正后P<0.05;红色表示高表达,绿色表示低表达。

2.2MDS相关DEGs的GO和KEGG富集分析 使用DAVID数据库和ggplot2包对上述DEGs进一步分析,GO富集分析结果显示,DEGs在生物学过程层面主要富集于免疫应答、信号通路传导等;在细胞组分层面,DEGs主要在细胞质膜外侧面执行功能;在分子功能层面主要是结合DNA及蛋白质。KEGG富集分析结果显示,DEGs主要富集于原发性免疫缺陷、造血细胞谱系、B细胞受体信号传导、癌症中的转录失调、Hippo信号通路、中性粒细胞细胞外陷阱形成。

2.3MDS的WGCNA分析 为了构建最佳的无尺度网络和有效的基因连接,最佳软阈值设置为4,见图2。随后将相似基因进行合并,共得到8个模块并分析其与临床性状的相关性。结果显示,黑色模块的基因重要性最强(P<0.01),并且与MDS呈负相关(cor=-0.51,P<0.01),与HD呈正相关(cor=0.51,P<0.01),黑色模块中一共有99种基因,根据geneSigFilter=0.5,moduleSigFilter=0.8两个条件筛选出来5个Hub基因。该模块的MM值与GS值的散点图见图3,右上角的基因即为Hub基因,分别是AKAP12、ARPP21、IGHV5-78、MME、NPY。

图2 软阈值的确认

图3 黑色模块基因模块隶属度与基因显著性相关的散点图

2.4Venn图绘制及关键基因的Lasso回归分析 Venn包取Hub基因与DEGs的交集,得出Hub基因包含于DEGs;Lasso回归分析对上述5个Hub基因进行分析,交叉验证误差最小的点所对应的值为疾病特征基因的数目,最终得到4个疾病特征基因,分别是AKAP12、ARPP21、MME、NPY。见图4。

图4 Hub基因与DEGs的Venn图及核心基因的Lasso回归分析

2.5训练集和测试集中特征基因的表达量 使用limma包与ggpubr包对训练集及测试集中疾病特征基因AKAP12、ARPP21、MME、NPY的表达进行循环和统计并可视化,在训练集中,疾病组上述4种特征基因的表达量下降(P<0.05),见图5;在测试集中也得出了相同的结果,见图6。上述结果显示,与健康对照组相比,这4种基因的表达在疾病组中明显下调(P<0.05)。

注:***P<0.01;HD为健康对照组,MDS为疾病组。

注:***P<0.01。

2.6ssGSEA分析 接下来对训练集中每个样本的数据进行免疫基因打分,进一步统计分析并对其进行小提琴图的绘制。在MDS中活化的B细胞、成熟B细胞、记忆B细胞、嗜酸性粒细胞的免疫浸润明显下降(P<0.01),见图7。而在上述4种特征基因中,免疫细胞相关性分析结果显示,特征基因与活化的B细胞、成熟的B细胞、记忆B细胞、嗜酸性粒细胞、Ⅱ型辅助T细胞、Treg细胞、活化的CD4+T细胞呈正相关;与效应记忆型CD8+T细胞、中央记忆T细胞(CD4+与CD8+)、成熟的树突状细胞,Ⅰ型辅助T细胞呈负相关,见图8。

图7 疾病组与健康对照组的免疫细胞差异分析

注:*P<0.05;**P<0.01;***P<0.001。

3 讨 论

MDS是一种常见的髓系异质性疾病,此类疾病的发生通常与AML的发生有密切联系。此类疾病一般发生于40岁之后,常见的病因有化疗或放疗、烟草或有毒溶剂或农药[7]。基因突变、基因表达失控、表观遗传学改变及免疫失调在此类疾病进程中起着重要作用[8-9]。在先天性免疫应答方面,MDS的造血干细胞中Toll样受体信号通路异常激活,细胞因子及炎症小体异常表达[10];在适应性免疫应答方面,各类T细胞和Treg细胞在疾病不同的阶段呈现不同的反应[11]。另外,有一部分基因的表达异常与免疫失调存在一定联系。笔者运用生物信息学方法挖掘到疾病的4个特征基因(AKAP12、ARPP21、MME、NPY),再对其进行免疫浸润分析,或许对此类疾病的发病机制和预后分析有一定价值。

AKAP12即A激酶锚定蛋白12,是蛋白激酶A与蛋白激酶C的调节剂,也是一种肿瘤生长调节的负调节因子[12]。其甲基化导致的表达下调在肺癌、结直肠癌、肝癌、食管癌等多种癌症中存在,与抑制癌细胞的血管生成、增殖、趋化和侵袭有密切关系。之前的研究显示,该分子不仅在青少年粒单核细胞白血病中呈现DNA高甲基化,而且在髓系恶性肿瘤(如AML、慢性髓性白血病、MDS)中也显著下调[13-14]。Gravin/AKAP12是AKAP的家族成员,其也被证实在急性白血病中低表达,已有研究发现,Gravin基因的低表达与AML的高复发风险之间存在一定联系[15]。而BOULTWOOD等[16]的研究表明,10例MDS患者中有9例患者的Gravin基因呈低表达,综上所述,本研究数据分析发现的AKAP12在MDS患者总体表达量降低与以往研究结果一致。

MME属于基质金属蛋白酶(MMP)家族(又称CD10),称为膜金属内肽酶,MME/CD10编码一种常见的急性淋巴细胞白血病抗原,是诊断儿童急性淋巴细胞白血病的关键细胞表面标志物[17]。同样,该基因在炎症调节、实体肿瘤的进展中有着重要意义。PIETRZAK等[18]的研究表明,MMP2、MMP16在AML中低表达,而CHUNG等[19]及ZHAO等[20]分别通过流式细胞实验和高通量的芯片分析技术证实MDS的骨髓细胞中的MME基因低表达,以上研究结果与本研究结果相符合。

ARPP21表达cAMP调节的磷酸化蛋白21,该蛋白富集于尾状核和小脑皮层,主要与智力发育障碍等疾病有关,另外它还可以作为microRNA的调节蛋白影响靶基因的表达[21];NPY基因编码一种在中枢神经系统中广泛表达的神经肽,神经肽Y受体是Gi/o蛋白偶联受体家族,可能是各种代谢和心血管疾病的危险因素。有研究表明,淋巴母细胞来源的NPY蛋白水平在儿童B细胞前体的急性白血病中增高,且与疾病的临床分型、预后的有利结果有关[22]。NILIUS-ELILIWI等[23]使用优于传统诊断方式的“光学基因组图谱”(一种广泛的基因组检测方式)检测到急性髓系白血病的DDX3X:MLLT10基因融合及ARPP21缺失。但ARPP21及NPY与MDS疾病的具体研究鲜见文献报道。在本研究中,ARPP21及NPY在MDS患者中低表达,且在训练集和测试集中结果一致。

另外,本研究对训练集的基因进行了免疫打分,发现MDS组中的活化、成熟、记忆的B细胞数量低于健康对照组。而对4个关键基因的免疫浸润分析中,笔者发现其与活化、成熟、记忆B细胞呈正相关(P<0.001),提示在MDS疾病中B细胞介导的体液免疫反应处于抑制的状态,这与之前的相关报道相符合[24]。而在细胞免疫应答中,之前的报道显示低危组的免疫系统处于激活、促炎状态;而高危组则与之相反。具体来说,低危组MDS患者的Th1、Ts(细胞毒性T细胞)数目增高,Th1/Th2比例高,Treg细胞减少;而高危组的Th1、Ts数量逐渐减少,Th1/Th2比例降低,Treg细胞增多。这有利于肿瘤细胞的免疫逃逸,以及疾病向白血病快速转化[11,25]。在本研究中,对关键基因的免疫浸润分析结果显示,AKAP12、ARPP21、MME、NPY与Th2、Treg细胞、活化的CD4+T细胞呈正相关,与Th1、Ts(CD8+效应T细胞)及中央型记忆T细胞呈负相关。提示这些关键基因或许对疾病预后分组有一定价值。

最后,MDS中免疫系统的激活与否直接影响到细胞因子的水平及炎症信号通路[10]。在差异基因的GO富集分析中,除了富集于免疫应答,还富集于信号通路的传导。WANG等[26]总结了MDS中的信号通路失调,其中NF-κB信号通路可以影响造血干细胞的功能。而B细胞受体信号传导的失调也许与B细胞活化抑制相关[27],深入探索关键基因与B细胞活化及免疫细胞变化的具体作用机制,或将为疾病的治疗与预后带来一定思考。

本研究主要运用WGCNA及Lasso回归的方式进行分析运算,筛选了MDS的疾病特征基因并对其与免疫细胞的关系进行了分析。本研究结果不仅证实了AKAP12及MME的水平在MDS中下降,更是发现了ARPP21及NPY在MDS中的水平变化,后续还需要收集临床标本进行多方面的验证。同时,在分析特征基因与免疫细胞之间的关系时,本研究尚未区分4种关键特征基因分别在高危组和低危组MDS中的变化,也未深入了解免疫细胞变化的机制,应该在后续研究中进一步探索。本研究为阐述MDS的发病机制、MDS与免疫细胞之间的关系、疾病的诊断及预后提供了重要的参考价值。