m6A RNA甲基化调节因子对肝癌影响的统计分析

喻文霞,陈鹏辉,周 霞,钟 琦,李慧敏

(云南民族大学 数学与计算机科学学院,云南 昆明 650500)

基因表观遗传学修饰是肿瘤发病机制研究中的热点. 研究表明,包括 DNA甲基化和组蛋白修饰在内的表观遗传学异常在肿瘤发生和发展中起着重要作用[1～3]. RNA甲基化主要发生在转录后水平调控基因表达,被认为是类似于DNA甲基化和组蛋白修饰的另一种表观遗传调控[3,4]. 已知RNA存在超过100种修饰,其中m6A甲基化是最广泛的 RNA甲基化修饰方式之一,贯穿于RNA生命周期的所有阶段,通过影响RNA代谢发挥生物学功能[1,3,5]. 研究发现,m6A RNA甲基化与肿瘤的增殖、分化、致瘤、侵袭性和转移性相关,在恶性肿瘤中具有癌基因和癌旁基因的作用[2,4]. 因此,研究m6A甲基化修饰,可以为患者治疗提供依据. 因为RNA 甲基化主要在转录后水平调控基因表达,转录本m6A 的修饰受甲基化转移酶、去甲基化酶和甲基化阅读蛋白的调控,所以研究m6A甲基化修饰因子(即m6A RNA甲基化调节因子)在疾病样本中的表达水平具有重要意义[2～10].

肝癌(liver cancer)是一种消化系统恶性疾病. 我国每年患肝癌的人在恶性肿瘤中排行第4,每年因肝癌死亡的人数约占世界肝癌死亡人数的40%. 目前我国对肝癌的治疗主要采取早期诊断和切除的方法,诊断方法主要是检查肝癌肿瘤标志物. 近期研究发现,m6A RNA甲基化调节因子YTHDF2、METTL3和METTL14表达水平的上调或下调与肝癌进展、转移及肝癌患者的预后有关[2,9,11],因此研究m6A RNA甲基化调节因子的表达水平为研究肝癌发生的表观遗传改变提供了一个新的维度.

通过系统分析21个m6A RNA 甲基化调节因子在465个肝和肝内胆管样本数据中的表达情况,通过与正常样本、肿瘤样本以及临床性状的相关性检验分析,发现绝大部分m6A RNA甲基化调节因子在肝组织正常组和肿瘤组的表达水平具有显著差异; 利用21个m6A RNA甲基化调节因子的表达水平将407个肿瘤样本分为两组,并对其分组结果做分型生存分析和分型临床相关性检验,结果提示m6A RNA甲基化调节因子的表达水平和肝和肝胆管疾病的发展具有显著的相关性.

1 数据和方法

1.1 数据

肝和肝内胆管表达水平数据和临床性状数据均来源于癌症数据共享系统GDC(Genomic Data Commons,https://portal.gdc.cancer.gov/)[3]. 该系统整合了包括TCGA在内的多个癌症数据库中的信息,提供了癌症数据的统一存储管理和展示,数据可靠且易于处理. 共提取了与m6A甲基化相关的21个因子的表达水平数据,其中包括58个正常样本和407个肿瘤样本. 提取了418个肝和肝内胆管临床样本数据,剔除信息缺失和分期不明确的数据,得到404个样本数据. 由于因子表达水平数据及临床样本数据均比较大,只列出了一部分数据,格式见表1和表2.

表1 21个m6A RNA甲基化调节因子的部分表达水平数据

1.2 m6A RNA甲基化调节因子

21个与m6A甲基化修饰相关的因子分别为: 甲基化转移酶(Writers): METTL3、METTL14、METTL16、WTAP、KIAA1492(Vir)、RBM15、RBM15B和ZC3H13; 去甲基化酶(Erasers): FTO和ALKBH5; 甲基化阅读蛋白(Readers): WTHDF2、YTHDF3、YTHDC1、YTHDC2、YTHDF2、HNRNPA2B1、elF3、IGF2BP1、IGF2BP2、IGF2BP3和HNRNPC[3].

1.3 logFC值和Wilcoxon检验

21个m6A RNA甲基化调节因子在肿瘤组和正常组的表达水平差异用指标logFC衡量,logFC计算方法见式(1):

(1)

式(1)中,若logFC>0,表示因子表达水平在肿瘤组较正常组上调; 若logFC<0,表示基因表达水平在肿瘤组较正常组下调; 若logFC=0,表示因子表达水平在两组样本中相等.

为更好地说明这些因子表达水平在两组样本中差异的统计显著性,同时引入Wilcoxon秩和检验. Wilcoxon秩和检验是一种非参数检验,常用于比较2组数据的差异性,检验时既考虑了2组数据差异的大小,又考虑了差异的方向. Wilcoxon秩和检验的统计量由式(2)给出:

(2)

其中,n1和n2分别表示正常组和肿瘤组样本容量,Wx表示混合正常组和肿瘤组样本后,正常组平均表达水平的秩和(即正常组平均表达水平在混合组的排序之和),τj表示第j个结值(即2组样本中表达平均表达水平相等)的个数. 由于n1=58,n2=407,均大于8,W近似服从N(0,1),因子表达水平差异大小可由P值度量.

取显著性水平为0.05,对21个m6A RNA甲基化调节因子在正常组和肿瘤组的肝和肝内胆管表达水平的差异性进行检验. 当|W|>1.64(P<0.05) 时,认为m6A RNA甲基化调节因子的表达水平在正常组和肿瘤组具有显著差异.

1.4 K-means聚类与肿瘤样本分型

目前,通常按照优势成份分型原则,即以肿瘤主要组织学类型(>50%的组织结构)对恶性肿瘤进行分型诊断[12]. 然而,我们的分析主要基于m6A RNA甲基化调节因子在肿瘤样本中的表达水平,优势成份分型原则对该类型数据并不适用. 因此,利用K-means聚类方法对肿瘤进行分型,选取的分型标准为: ①分组中没有样本特别小的组; ②组内相关性高,组间相关性低; ③共识聚类累积分布函数CDF增速平缓.

欧式距离作为相似性的评价指标， 欧氏距离为:

(3)

式(3)中,K值从2到9依次选择并进行比较;n=21,表示21个m6A RNA甲基化调节因子; 聚类在肿瘤组的407个样本数据上进行.

1.5 生存分析

生存分析主要应用于分析临床数据上患者的生存时间,包括临床特征对患者的影响,患者个体的几年内生存率,高低风险组之间的生存率差异等. 生存分析的方法有多种,我们主要通过对分型结果绘制Kaplan-meier生存曲线,即以生存时间为横轴,以生存率为纵轴,绘制连续型阶梯型曲线进行生存分析.

利用时序检验(log-rank test)比较各组Kaplan-meier生存曲线是否具有差异. 对2组生存曲线λ1(t)和λ2(t),提出假设:

H0:λ1(t)=λ2(t)，H1:λ1(t)≠λ2(t).

利用式(4)的χ2统计量进行检验:

(4)

其中,di,j表示第i组第j个时刻的实际死亡数,Ti,j表示第i组第j个时刻的理论死亡数,Vi,j表示第i组第j个时刻的实际死亡数的方差,Ti,j和Vi,j的具体计算方法由文献[13]给出[13].

2 结果分析

2.1 m6A RNA甲基化调节因子的表达分析

计算logFC值,并对21个m6A RNA甲基化调节因子在正常组和肿瘤组的表达水平进行Wilcoxon检验, 结果见表3.

表3 21个m6A RNA甲基化调节因子的表达量差异

由表3可知,除因子ZC3H13在正常组和肿瘤组的表达水平不具有显著差异(P=0.806>0.05)外,其余20个因子在正常组和肿瘤组的表达水平均具有显著差异(P<0.05),并且与正常组相比,其因子表达水平在肿瘤组中显著上调(logFC>0).

为了更直观地显示21个m6A RNA甲基化调节因子在正常组与肿瘤组的表达差异情况,引用R语言中的pheatmap包绘制它们的表达水平热图(图1). 图中横坐标为m6A甲基化相关的465个样本(Type:N为蓝色，表示58个正常样本；T为红色，表示407个肿瘤样本),纵坐标标明了21个甲基化调节因子名称，图中绿色代表低表达,白色代表中表达,,红色代表高表达;“A***”表示P<0.001,“A**”表示P<0.01,“A*”表示P<0.05(A表示因子名). 从图1可以看出,在正常组中,因子ZC3H13的表达水平情况为白色,其余因子的表达水平情况大部分显示为绿色,由此可知,除ZC3H13因子是中表达外,其余19个因子在正常样本里基本都是低表达; 而在肿瘤样本中,除因子ZC3H13和METTL14之外,其余19个因子在肿瘤样本中显著都是高表达(P<0.01). 由此也提示肝和肝内胆管的癌变或许与m6A RNA甲基化调节因子表达水平上调有关.

图1 21个m6A RNA甲基化调节因子的表达水平热图

检验因子表达水平的相关性,可以为之后进行风险生存分析,建立模型预测生存时间及疾病的治疗提供依据. 为了更好地理解21个m6A RNA甲基化调节因子之间的相互作用,计算它们表达水平之间的Pearson相关系数,结果见图2.

图2 m6A RNA甲基化调节因子之间的相关性图

从图2可以看出,HNRNPC和HNRNPA2B1表达水平之间的相关系数高达0.8,表明这2个因子之间高度正相关,当因子HNRNPC在肝癌中高表达时,HNRNPA2B1也高表达; YTHDC1与METTL16、HNRNPA2B1、YTHDC2、HNRNPC、RBM15B、YTHDF1、METTL3、WTAP等因子的表达水平之间的相关系数均大于0.5,由此可知YTHDC1与上述7个甲基化因子呈正相关; YTHDF1与RBM15B、METTL3、HNRNPC、HNRNPA2B1、WTAP等因子的表达水平之间的相关系数均大于0.5,因此YTHDF1与上述5个甲基化因子正相关; ZC3H13和FTO与ELF3、IGF2BP3、IGF2BP2、IGF2BP1、ALKBH5等因子间相关性较弱; 其余少部分因子之间也不相关. 由上结果分析表明,大部分的m6A RNA甲基化调节因子的表达水平在肝癌中相互正相关.

2.2 肿瘤样本分型结果

基于21个m6A RNA调节因子的表达水平,对肝癌肿瘤样本进行分型,结果见图3. 根据分型标准,当k=4时,CDF增速平缓,曲线下的面积变化开始变小,因此可以将肿瘤样本最多分为4组,即k=2,3,4(图3a,图3b); 但由于k=3,k=4时, 从图3d可以看出,组1和组3、组2都具有较高的相关性,且在图3(e)中,k=4,存在样本较小的组(组3),组与组之间的相关性相关性较高,因此不满足分型标准; 当k=2时满足所有分型标准,可以将肿瘤样本分为2组(3c); 图3f是利用主成分分析的方法对k=2这个分型结果进行检验所得,由图中可以清楚的看到红色和蓝色的点重复部分较少,分型结果较好,因此选取k=2为最适当的结果. 同时,结合基因的表达水平,21个m6A RNA甲基化调节因子中绝大部分在cluster1中的表达水平低于在cluster2中的表达水平,分型具有明显意义.

图3 m6A甲基化肿瘤分型结果图

2.3 分型与临床结果的相关性检验

对407个肝癌肿瘤样本的分型结果和404个只含有生存时间和生存状态的临床数据做分型生存分析,结果见图4.

在图4中,横坐标为生存时间,纵坐标为生存率. 由该图可以看出,cluster1和cluster2之间的生存率具有较大的差异(p=2.428×10-5). cluster1的患者生存率显著高于cluster2,且生存时间曲线长于cluster2,即cluster1的总生存周期较cluster2更长,说明肝癌患者生存率和生存时间跟分型有较高的相关性. 5年生存率是癌症患者较关心的一个问题,从图中也可以看出2组之间5年的生存率存在很大的差异. 由于m6A RNA调节因子在cluster1中的表达水平低于在cluster2中的表达水平,因此可以得出m6A RNA调节因子的表达水平与肝癌患者生存时间相关. 即因子表达水平越高,患者生存时间越短,5年生存率越低,反映了m6A RNA甲基化调节因子表达水平与肿瘤发展程度有关.

基于卡方检验(显著性水平α=0.05)对分型结果和临床性状(生存状态、性别、年龄、分级、分期、原发肿瘤(T分期)、远处转移(M分期)、淋巴结(N分期))进行相关性检验,其结果如图5所示.图中横坐标为肿瘤分型样本cluster； 纵坐标为各临床性状和21个甲基化因子的名称；因子名称所对应的热图中绿色代表低表达,白色代表中表达,,红色代表高表达;“B***”表示P<0.001,“B*”表示P<0.01,“B*”表示P<0.05(B表示临床性状).

由图5可以看到临床性状分级(grade)、临床性状患者的性别(gender)、分期(stage)、N分期(N)、T分期(T)和生存状态(fustat)与肿瘤的分型均具有显著相关性(P<0.05); 患者年龄(age)、M分期(M)与肿瘤的分型相关性不显著(P>0.05); 从因子的表达水平上看,21个m6A RNA甲基化调节因子中绝大部分因子在cluster2中高表达,而在cluster1中基本都是低表达. 因此也可以得出,甲基化调节因子的表达水平与临床性状相关,也与肝癌的严重程度有关.

3 结语

癌症一直是威胁人类生命健康的最严重的疾病之一,同时也是生物学家研究的重点. 基因表观遗传学修饰是肿瘤发病机制研究中的热点. 研究表明,包括DNA甲基化和组蛋白修饰在内的表观遗传学异常在肿瘤发生和发展中起着重要作用. 基于肝和肝内胆管疾病中的转录组数据和临床数据,通过研究21个m6A RNA甲基化调节因子在肝和肝内胆管样本中的表达水平,发现大部分m6A RNA甲基化调节因子在癌症样本和正常样本中的表达情况具有显著差异: 大部分甲基化调节因子的表达水平在肝癌样本中相对于正常样本显著上调,这表明m6A RNA甲基化调节因子的表达水平与肝癌发生相关. 此外,还对m6A RNA甲基化调节因子在407个肝和肝内胆管癌症样本中的表达水平分型,并进行分型生存分析和分型临床相关性检验,发现肿瘤的分型结果与患者生存状况以及患者大部分临床性状具有高度相关性. 我们的分析结果提示,m6A RNA甲基化调节因子的表达水平或许为研究肝癌发生的表观遗传改变提供了一个新的维度,研究m6A RNA甲基化调节因子表达水平的差异性和临床性状的相关性,可为今后检测m6A RNA甲基化调节因子在肝和肝内胆管癌症中的作用提供基础,也为肝癌患者后期的治疗和生存分析提供了理论依据.