m6A RNA甲基化調(diào)節(jié)因子對肝癌影響的統(tǒng)計分析

喻文霞,陳鵬輝,周 霞,鐘 琦,李慧敏

(云南民族大學(xué) 數(shù)學(xué)與計算機科學(xué)學(xué)院,云南 昆明 650500)

基因表觀遺傳學(xué)修飾是腫瘤發(fā)病機制研究中的熱點. 研究表明,包括 DNA甲基化和組蛋白修飾在內(nèi)的表觀遺傳學(xué)異常在腫瘤發(fā)生和發(fā)展中起著重要作用[1～3]. RNA甲基化主要發(fā)生在轉(zhuǎn)錄后水平調(diào)控基因表達(dá),被認(rèn)為是類似于DNA甲基化和組蛋白修飾的另一種表觀遺傳調(diào)控[3,4]. 已知RNA存在超過100種修飾,其中m6A甲基化是最廣泛的 RNA甲基化修飾方式之一,貫穿于RNA生命周期的所有階段,通過影響RNA代謝發(fā)揮生物學(xué)功能[1,3,5]. 研究發(fā)現(xiàn),m6A RNA甲基化與腫瘤的增殖、分化、致瘤、侵襲性和轉(zhuǎn)移性相關(guān),在惡性腫瘤中具有癌基因和癌旁基因的作用[2,4]. 因此,研究m6A甲基化修飾,可以為患者治療提供依據(jù). 因為RNA 甲基化主要在轉(zhuǎn)錄后水平調(diào)控基因表達(dá),轉(zhuǎn)錄本m6A 的修飾受甲基化轉(zhuǎn)移酶、去甲基化酶和甲基化閱讀蛋白的調(diào)控,所以研究m6A甲基化修飾因子(即m6A RNA甲基化調(diào)節(jié)因子)在疾病樣本中的表達(dá)水平具有重要意義[2～10].

肝癌(liver cancer)是一種消化系統(tǒng)惡性疾病. 我國每年患肝癌的人在惡性腫瘤中排行第4,每年因肝癌死亡的人數(shù)約占世界肝癌死亡人數(shù)的40%. 目前我國對肝癌的治療主要采取早期診斷和切除的方法,診斷方法主要是檢查肝癌腫瘤標(biāo)志物. 近期研究發(fā)現(xiàn),m6A RNA甲基化調(diào)節(jié)因子YTHDF2、METTL3和METTL14表達(dá)水平的上調(diào)或下調(diào)與肝癌進(jìn)展、轉(zhuǎn)移及肝癌患者的預(yù)后有關(guān)[2,9,11],因此研究m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平為研究肝癌發(fā)生的表觀遺傳改變提供了一個新的維度.

通過系統(tǒng)分析21個m6A RNA 甲基化調(diào)節(jié)因子在465個肝和肝內(nèi)膽管樣本數(shù)據(jù)中的表達(dá)情況,通過與正常樣本、腫瘤樣本以及臨床性狀的相關(guān)性檢驗分析,發(fā)現(xiàn)絕大部分m6A RNA甲基化調(diào)節(jié)因子在肝組織正常組和腫瘤組的表達(dá)水平具有顯著差異; 利用21個m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平將407個腫瘤樣本分為兩組,并對其分組結(jié)果做分型生存分析和分型臨床相關(guān)性檢驗,結(jié)果提示m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平和肝和肝膽管疾病的發(fā)展具有顯著的相關(guān)性.

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)

肝和肝內(nèi)膽管表達(dá)水平數(shù)據(jù)和臨床性狀數(shù)據(jù)均來源于癌癥數(shù)據(jù)共享系統(tǒng)GDC(Genomic Data Commons,https://portal.gdc.cancer.gov/)[3]. 該系統(tǒng)整合了包括TCGA在內(nèi)的多個癌癥數(shù)據(jù)庫中的信息,提供了癌癥數(shù)據(jù)的統(tǒng)一存儲管理和展示,數(shù)據(jù)可靠且易于處理. 共提取了與m6A甲基化相關(guān)的21個因子的表達(dá)水平數(shù)據(jù),其中包括58個正常樣本和407個腫瘤樣本. 提取了418個肝和肝內(nèi)膽管臨床樣本數(shù)據(jù),剔除信息缺失和分期不明確的數(shù)據(jù),得到404個樣本數(shù)據(jù). 由于因子表達(dá)水平數(shù)據(jù)及臨床樣本數(shù)據(jù)均比較大,只列出了一部分?jǐn)?shù)據(jù),格式見表1和表2.

表1 21個m6A RNA甲基化調(diào)節(jié)因子的部分表達(dá)水平數(shù)據(jù)

1.2 m6A RNA甲基化調(diào)節(jié)因子

21個與m6A甲基化修飾相關(guān)的因子分別為: 甲基化轉(zhuǎn)移酶(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甲基化調(diào)節(jié)因子在腫瘤組和正常組的表達(dá)水平差異用指標(biāo)logFC衡量,logFC計算方法見式(1):

(1)

式(1)中,若logFC>0,表示因子表達(dá)水平在腫瘤組較正常組上調(diào); 若logFC<0,表示基因表達(dá)水平在腫瘤組較正常組下調(diào); 若logFC=0,表示因子表達(dá)水平在兩組樣本中相等.

為更好地說明這些因子表達(dá)水平在兩組樣本中差異的統(tǒng)計顯著性,同時引入Wilcoxon秩和檢驗. Wilcoxon秩和檢驗是一種非參數(shù)檢驗,常用于比較2組數(shù)據(jù)的差異性,檢驗時既考慮了2組數(shù)據(jù)差異的大小,又考慮了差異的方向. Wilcoxon秩和檢驗的統(tǒng)計量由式(2)給出:

(2)

其中,n1和n2分別表示正常組和腫瘤組樣本容量,Wx表示混合正常組和腫瘤組樣本后,正常組平均表達(dá)水平的秩和(即正常組平均表達(dá)水平在混合組的排序之和),τj表示第j個結(jié)值(即2組樣本中表達(dá)平均表達(dá)水平相等)的個數(shù). 由于n1=58,n2=407,均大于8,W近似服從N(0,1),因子表達(dá)水平差異大小可由P值度量.

取顯著性水平為0.05,對21個m6A RNA甲基化調(diào)節(jié)因子在正常組和腫瘤組的肝和肝內(nèi)膽管表達(dá)水平的差異性進(jìn)行檢驗. 當(dāng)|W|>1.64(P<0.05) 時,認(rèn)為m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平在正常組和腫瘤組具有顯著差異.

1.4 K-means聚類與腫瘤樣本分型

目前,通常按照優(yōu)勢成份分型原則,即以腫瘤主要組織學(xué)類型(>50%的組織結(jié)構(gòu))對惡性腫瘤進(jìn)行分型診斷[12]. 然而,我們的分析主要基于m6A RNA甲基化調(diào)節(jié)因子在腫瘤樣本中的表達(dá)水平,優(yōu)勢成份分型原則對該類型數(shù)據(jù)并不適用. 因此,利用K-means聚類方法對腫瘤進(jìn)行分型,選取的分型標(biāo)準(zhǔn)為: ①分組中沒有樣本特別小的組; ②組內(nèi)相關(guān)性高,組間相關(guān)性低; ③共識聚類累積分布函數(shù)CDF增速平緩.

歐式距離作為相似性的評價指標(biāo)， 歐氏距離為:

(3)

式(3)中,K值從2到9依次選擇并進(jìn)行比較;n=21,表示21個m6A RNA甲基化調(diào)節(jié)因子; 聚類在腫瘤組的407個樣本數(shù)據(jù)上進(jìn)行.

1.5 生存分析

生存分析主要應(yīng)用于分析臨床數(shù)據(jù)上患者的生存時間,包括臨床特征對患者的影響,患者個體的幾年內(nèi)生存率,高低風(fēng)險組之間的生存率差異等. 生存分析的方法有多種,我們主要通過對分型結(jié)果繪制Kaplan-meier生存曲線,即以生存時間為橫軸,以生存率為縱軸,繪制連續(xù)型階梯型曲線進(jìn)行生存分析.

利用時序檢驗(log-rank test)比較各組Kaplan-meier生存曲線是否具有差異. 對2組生存曲線λ1(t)和λ2(t),提出假設(shè):

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

利用式(4)的χ2統(tǒng)計量進(jìn)行檢驗:

(4)

其中,di,j表示第i組第j個時刻的實際死亡數(shù),Ti,j表示第i組第j個時刻的理論死亡數(shù),Vi,j表示第i組第j個時刻的實際死亡數(shù)的方差,Ti,j和Vi,j的具體計算方法由文獻(xiàn)[13]給出[13].

2 結(jié)果分析

2.1 m6A RNA甲基化調(diào)節(jié)因子的表達(dá)分析

計算logFC值,并對21個m6A RNA甲基化調(diào)節(jié)因子在正常組和腫瘤組的表達(dá)水平進(jìn)行Wilcoxon檢驗, 結(jié)果見表3.

表3 21個m6A RNA甲基化調(diào)節(jié)因子的表達(dá)量差異

由表3可知,除因子ZC3H13在正常組和腫瘤組的表達(dá)水平不具有顯著差異(P=0.806>0.05)外,其余20個因子在正常組和腫瘤組的表達(dá)水平均具有顯著差異(P<0.05),并且與正常組相比,其因子表達(dá)水平在腫瘤組中顯著上調(diào)(logFC>0).

為了更直觀地顯示21個m6A RNA甲基化調(diào)節(jié)因子在正常組與腫瘤組的表達(dá)差異情況,引用R語言中的pheatmap包繪制它們的表達(dá)水平熱圖(圖1). 圖中橫坐標(biāo)為m6A甲基化相關(guān)的465個樣本(Type:N為藍(lán)色，表示58個正常樣本；T為紅色，表示407個腫瘤樣本),縱坐標(biāo)標(biāo)明了21個甲基化調(diào)節(jié)因子名稱，圖中綠色代表低表達(dá),白色代表中表達(dá),,紅色代表高表達(dá);“A***”表示P<0.001,“A**”表示P<0.01,“A*”表示P<0.05(A表示因子名). 從圖1可以看出,在正常組中,因子ZC3H13的表達(dá)水平情況為白色,其余因子的表達(dá)水平情況大部分顯示為綠色,由此可知,除ZC3H13因子是中表達(dá)外,其余19個因子在正常樣本里基本都是低表達(dá); 而在腫瘤樣本中,除因子ZC3H13和METTL14之外,其余19個因子在腫瘤樣本中顯著都是高表達(dá)(P<0.01). 由此也提示肝和肝內(nèi)膽管的癌變或許與m6A RNA甲基化調(diào)節(jié)因子表達(dá)水平上調(diào)有關(guān).

圖1 21個m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平熱圖

檢驗因子表達(dá)水平的相關(guān)性,可以為之后進(jìn)行風(fēng)險生存分析,建立模型預(yù)測生存時間及疾病的治療提供依據(jù). 為了更好地理解21個m6A RNA甲基化調(diào)節(jié)因子之間的相互作用,計算它們表達(dá)水平之間的Pearson相關(guān)系數(shù),結(jié)果見圖2.

圖2 m6A RNA甲基化調(diào)節(jié)因子之間的相關(guān)性圖

從圖2可以看出,HNRNPC和HNRNPA2B1表達(dá)水平之間的相關(guān)系數(shù)高達(dá)0.8,表明這2個因子之間高度正相關(guān),當(dāng)因子HNRNPC在肝癌中高表達(dá)時,HNRNPA2B1也高表達(dá); YTHDC1與METTL16、HNRNPA2B1、YTHDC2、HNRNPC、RBM15B、YTHDF1、METTL3、WTAP等因子的表達(dá)水平之間的相關(guān)系數(shù)均大于0.5,由此可知YTHDC1與上述7個甲基化因子呈正相關(guān); YTHDF1與RBM15B、METTL3、HNRNPC、HNRNPA2B1、WTAP等因子的表達(dá)水平之間的相關(guān)系數(shù)均大于0.5,因此YTHDF1與上述5個甲基化因子正相關(guān); ZC3H13和FTO與ELF3、IGF2BP3、IGF2BP2、IGF2BP1、ALKBH5等因子間相關(guān)性較弱; 其余少部分因子之間也不相關(guān). 由上結(jié)果分析表明,大部分的m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平在肝癌中相互正相關(guān).

2.2 腫瘤樣本分型結(jié)果

基于21個m6A RNA調(diào)節(jié)因子的表達(dá)水平,對肝癌腫瘤樣本進(jìn)行分型,結(jié)果見圖3. 根據(jù)分型標(biāo)準(zhǔn),當(dāng)k=4時,CDF增速平緩,曲線下的面積變化開始變小,因此可以將腫瘤樣本最多分為4組,即k=2,3,4(圖3a,圖3b); 但由于k=3,k=4時, 從圖3d可以看出,組1和組3、組2都具有較高的相關(guān)性,且在圖3(e)中,k=4,存在樣本較小的組(組3),組與組之間的相關(guān)性相關(guān)性較高,因此不滿足分型標(biāo)準(zhǔn); 當(dāng)k=2時滿足所有分型標(biāo)準(zhǔn),可以將腫瘤樣本分為2組(3c); 圖3f是利用主成分分析的方法對k=2這個分型結(jié)果進(jìn)行檢驗所得,由圖中可以清楚的看到紅色和藍(lán)色的點重復(fù)部分較少,分型結(jié)果較好,因此選取k=2為最適當(dāng)?shù)慕Y(jié)果. 同時,結(jié)合基因的表達(dá)水平,21個m6A RNA甲基化調(diào)節(jié)因子中絕大部分在cluster1中的表達(dá)水平低于在cluster2中的表達(dá)水平,分型具有明顯意義.

圖3 m6A甲基化腫瘤分型結(jié)果圖

2.3 分型與臨床結(jié)果的相關(guān)性檢驗

對407個肝癌腫瘤樣本的分型結(jié)果和404個只含有生存時間和生存狀態(tài)的臨床數(shù)據(jù)做分型生存分析,結(jié)果見圖4.

在圖4中,橫坐標(biāo)為生存時間,縱坐標(biāo)為生存率. 由該圖可以看出,cluster1和cluster2之間的生存率具有較大的差異(p=2.428×10-5). cluster1的患者生存率顯著高于cluster2,且生存時間曲線長于cluster2,即cluster1的總生存周期較cluster2更長,說明肝癌患者生存率和生存時間跟分型有較高的相關(guān)性. 5年生存率是癌癥患者較關(guān)心的一個問題,從圖中也可以看出2組之間5年的生存率存在很大的差異. 由于m6A RNA調(diào)節(jié)因子在cluster1中的表達(dá)水平低于在cluster2中的表達(dá)水平,因此可以得出m6A RNA調(diào)節(jié)因子的表達(dá)水平與肝癌患者生存時間相關(guān). 即因子表達(dá)水平越高,患者生存時間越短,5年生存率越低,反映了m6A RNA甲基化調(diào)節(jié)因子表達(dá)水平與腫瘤發(fā)展程度有關(guān).

基于卡方檢驗(顯著性水平α=0.05)對分型結(jié)果和臨床性狀(生存狀態(tài)、性別、年齡、分級、分期、原發(fā)腫瘤(T分期)、遠(yuǎn)處轉(zhuǎn)移(M分期)、淋巴結(jié)(N分期))進(jìn)行相關(guān)性檢驗,其結(jié)果如圖5所示.圖中橫坐標(biāo)為腫瘤分型樣本cluster； 縱坐標(biāo)為各臨床性狀和21個甲基化因子的名稱；因子名稱所對應(yīng)的熱圖中綠色代表低表達(dá),白色代表中表達(dá),,紅色代表高表達(dá);“B***”表示P<0.001,“B*”表示P<0.01,“B*”表示P<0.05(B表示臨床性狀).

由圖5可以看到臨床性狀分級(grade)、臨床性狀患者的性別(gender)、分期(stage)、N分期(N)、T分期(T)和生存狀態(tài)(fustat)與腫瘤的分型均具有顯著相關(guān)性(P<0.05); 患者年齡(age)、M分期(M)與腫瘤的分型相關(guān)性不顯著(P>0.05); 從因子的表達(dá)水平上看,21個m6A RNA甲基化調(diào)節(jié)因子中絕大部分因子在cluster2中高表達(dá),而在cluster1中基本都是低表達(dá). 因此也可以得出,甲基化調(diào)節(jié)因子的表達(dá)水平與臨床性狀相關(guān),也與肝癌的嚴(yán)重程度有關(guān).

3 結(jié)語

癌癥一直是威脅人類生命健康的最嚴(yán)重的疾病之一,同時也是生物學(xué)家研究的重點. 基因表觀遺傳學(xué)修飾是腫瘤發(fā)病機制研究中的熱點. 研究表明,包括DNA甲基化和組蛋白修飾在內(nèi)的表觀遺傳學(xué)異常在腫瘤發(fā)生和發(fā)展中起著重要作用. 基于肝和肝內(nèi)膽管疾病中的轉(zhuǎn)錄組數(shù)據(jù)和臨床數(shù)據(jù),通過研究21個m6A RNA甲基化調(diào)節(jié)因子在肝和肝內(nèi)膽管樣本中的表達(dá)水平,發(fā)現(xiàn)大部分m6A RNA甲基化調(diào)節(jié)因子在癌癥樣本和正常樣本中的表達(dá)情況具有顯著差異: 大部分甲基化調(diào)節(jié)因子的表達(dá)水平在肝癌樣本中相對于正常樣本顯著上調(diào),這表明m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平與肝癌發(fā)生相關(guān). 此外,還對m6A RNA甲基化調(diào)節(jié)因子在407個肝和肝內(nèi)膽管癌癥樣本中的表達(dá)水平分型,并進(jìn)行分型生存分析和分型臨床相關(guān)性檢驗,發(fā)現(xiàn)腫瘤的分型結(jié)果與患者生存狀況以及患者大部分臨床性狀具有高度相關(guān)性. 我們的分析結(jié)果提示,m6A RNA甲基化調(diào)節(jié)因子的表達(dá)水平或許為研究肝癌發(fā)生的表觀遺傳改變提供了一個新的維度,研究m6A RNA甲基化調(diào)節(jié)因子表達(dá)水平的差異性和臨床性狀的相關(guān)性,可為今后檢測m6A RNA甲基化調(diào)節(jié)因子在肝和肝內(nèi)膽管癌癥中的作用提供基礎(chǔ),也為肝癌患者后期的治療和生存分析提供了理論依據(jù).