羅斯琦，婁世鋒

(重慶醫(yī)科大學附屬第二醫(yī)院血液內(nèi)科，重慶 400010)

1 研究背景與思路分析

急性早幼粒細白血病(acute promyelocytic leukemia，APL)是一種以骨髓分化在早幼粒細胞階段停滯導致細胞蓄積為特征的急性髓系細胞白血病(acute myelocytic leukemia，AML)，在WHO分型中屬M3型急性髓系白血病(AML-M3)，APL患者t(15；17)(q22；q21)染色體易位產(chǎn)生PML-RARα融合基因[1]，該融合基因是導致APL中骨髓細胞分化停滯的主要機制[2]。全反式維甲酸(all-trans-retinoic acid，ATRA)聯(lián)合三氧化二砷對APL患者的誘導分化治療是目前常用的靶向治療藥物，治愈率高，可以將完全緩解(completeremission，CR)率 提 高 到90%[3]，ATRA在APL細 胞 中 誘 導 分 化 的 機制包括：(1)ATRA將PML-RARa從轉(zhuǎn)錄阻遏物轉(zhuǎn)化為轉(zhuǎn)錄激活物[4]；(2)ATRA通過半胱氨酸蛋白酶或泛素/蛋白酶體系統(tǒng)降解PML-RARa融合蛋白，同時釋放野生型RXR和PML，并誘導PMLNBs重組，PML/RARalpha直接參與賦予APL細胞ATRA敏感性，并 且ATRA誘導的PMLNBs重 組是PML/RARalpha消 失 的 結(jié) 果[5]；(3)RAF-1/MEK/ERK信 號通路被激活并通過調(diào)節(jié)C/EBPb，C/EBPe和PU.1的蛋白質(zhì)水平來調(diào)節(jié)ATRA誘導的APL細胞分化[6]，第四為一個非基因組途徑，ATRA激活PKA從而導致ATRA靶基因的轉(zhuǎn)錄增加和APL的分化[7]。在ATRA治療過程中，可能出現(xiàn)分化綜合征(Differentiationsyndrome，DS)，又稱維甲酸綜合征(Retinoicacidsyndrome，RAS)[8]，同時仍有ATRA治療后復發(fā)的情況發(fā)生，減少APL患者的早期死亡率以及減少復發(fā)率，是APL治療的研究熱點。因此，本研究旨在通過生物信息學工具分析ATRA處理APLNB4細胞(一種從人急性早幼粒細胞白血病(M3)中分離的具有t(15；17)標記的成熟誘導細胞系，NB4)前后基因表達譜差異，并尋找差異基因富集的信號途徑，探討全反式維甲酸治療急性早幼粒細白血病的可能分子作用機制，同時篩選出信號通路中高表達提示預后不良的關(guān)鍵基因，為尋找APL治療的新靶點提供思路。

2 資料和方法

2.1 一般資料

以美國國立生物技術(shù)信息中心(NCBI)的高通量基因表達數(shù)據(jù)庫(GEO)為工具，以“all-trans-retinoic acid”為關(guān)鍵詞，選取并下載下載ATRA分化的野生型和TG2基因敲除NB4細胞的基因表達譜數(shù)據(jù)集(GSE23702)，其來源于德布勒森大學，采用平臺為GPL6244[HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]。選取其中ATRA處理前、48h、72h后的野生型APL NB4細胞基因表達數(shù)據(jù)，每個時間點選取3個樣本。

2.2 方法

2.2.1 基因芯片質(zhì)量驗證

使用GEO2R(https：//www.ncbi.nlm.nih.gov/geo/geo2r) R boxplot對表達式數(shù)據(jù)應用分位數(shù)生成箱線圖，如果若某個芯片的箱線圖出現(xiàn)較大偏差則表明該芯片存在問題。使用R limma (plotDensities)生成表達密度圖，如果圖中曲線相差偏離大，則在差異表達分析之前進行數(shù)據(jù)歸一化。使用limma(qqt) 生成t統(tǒng)計分位線圖評估limma測試結(jié)果的質(zhì)量。

2.2.2 差異基因篩選

在GEO2R使 用Bioconductor項 目 中 的GEOquery和limma R包對分析篩選GEO數(shù)據(jù)集中不同數(shù)據(jù)組的差異表達水平[9]。在GSE23702數(shù)據(jù)集中以未經(jīng)過ATRA處理過的野生型APL NB4細胞基因表達數(shù)據(jù)(GSM594067、GSM594068、GSM594069)為對照組分別與ATRA處理48h(GSM594073、GSM594074、GSM594075)、72h(GSM594079、GSM594080、GSM594081)后的野生型APL NB4細胞基因表達數(shù)據(jù)對比，篩選標準為矯正后P＜0.05并且l log(Fc)I＞2，使用GEO2R在線工具limma(volcanoplot)包生成火山圖，并顯示統(tǒng)計顯著性(-log10P值)與變化幅度(log2倍變化)的對比，使差異基因可視化，最后使用limma(vennDiagram)包生成維恩圖，用于探索和下載多個對比之間重要基因的重疊，可聯(lián)合使用在線工具Draw Venn Diagram(http：//bioinformatics.psb.ugent.be/webtools/Venn/)篩選出在ATRA處理48h和72h后同時上調(diào)或下調(diào)的DEGs用于進一步分析。

2.2.3 基因本體論(gene ontology，GO)富集分析和京都基因與基因組百科全書(Kyotoencyclopedia of genes and genomes，KEGG)通路篩選

用DAVID 6.8在線數(shù)據(jù)庫(https：//david.ncifcrf.gov/tools.jsp)對DEGs進行基因本體論(gene ontology，GO)注釋和京都基因與基因組百科全書(Kyotoencyclopedia of genes and genomes，KEGG)通路分析[10]，聯(lián)合基于KEGG的注釋系統(tǒng)KOBAS(http：//kobas.cbi.pku.edu.cn/kobas3/?t=1)尋找差異基因富集的關(guān)鍵信號通路[11]。

2.2.4 PPI網(wǎng)絡(luò)構(gòu)建和分析

使 用String在 線 數(shù) 據(jù) 庫(https：//string-db.org/)對 上 述DEGs進行分析，構(gòu)建PPI網(wǎng)絡(luò)[12]后下載并使用Cytoscape(Version3.8.2)軟件進行可視化[13]，利用cytoHubba插件分析預測該網(wǎng)絡(luò)中最高相關(guān)度的前10個DEGs作為ATRA誘導細胞分化相關(guān)的關(guān)鍵基因，即hub基因[14]。

2.2.5 統(tǒng)計學分析

運用GEO2R工具的Benjamini & Hochberg法分析GSE23702表達譜芯片數(shù)據(jù)，采用SAM法篩選芯片數(shù)據(jù)中DEGs。在DAVID數(shù)據(jù)庫中，通過Bonferroni 校正法和Bootstrap 法等進行GO分析的檢驗，利用Fisher精確概率檢驗和基因富集分析等方法進行KEGG分析。在Cytoscape中，選取CytoHubba插件中MCC拓撲分析算法進行數(shù)據(jù)分析。

3 結(jié)果

3.1 芯片質(zhì)量驗證

用GEO2R對表達式數(shù)據(jù)生成箱線圖(圖1)，其中以中位數(shù)為中心的值表示數(shù)據(jù)已標準化且可交叉比較；使用R limma(plotDensities)生成表達密度圖(圖2)，圖中9條曲線之間趨勢相對平行，使用limma(qqt) 生成t統(tǒng)計分位線圖(圖3)，其中圈基本沿一條直線走行，說明limma測試結(jié)果的質(zhì)量較好，該芯片質(zhì)量優(yōu)，檢測結(jié)果穩(wěn)定，可為后續(xù)分析提供可靠數(shù)據(jù)。

圖1 箱線圖

圖2 表達密度圖

圖3 t統(tǒng)計分位線圖

3.2 DEGs篩選

用GEO2R使用limma(vennDiagram)對3組數(shù)據(jù)：ATRA處理 前(UNTREATED)、48h(ATRA48H)、72h(ATRA72H)后 的 野 生型APL NB4細胞基因表達數(shù)據(jù)進行綜合分析(設(shè)置P＜0.05)并生成韋恩圖(圖4)，可以看出ATRA處理48h后(ATRA48H)與未處理組的差異基因與處理72h(ATRA72H)后與未處理組的差異基因有相同部分，但ATRA處理48h與ATRA處理72h組無差別，因此我們分別選取ATRA48H、ATRA72H與UNTREATED組進行差異分析，用P＜0.05并且l log(Fc)I＞2進行篩選；與未用ATRA處理相比ATRA處理48h后篩選出上調(diào)基因140個，下調(diào)基因96個，共計236個 DEGs；ATRA處理72h后篩選出上調(diào)基因191個，下調(diào)基因131個，共計322個 DEGs，結(jié)果分別用火山圖(圖5和圖6)進行展示，上調(diào)和下調(diào)最顯著的前10個DEGs具體信息分別整理為表格(表1和表2)，同時兩組數(shù)據(jù)上調(diào)和下調(diào)最顯著的前10個 DEGs取交集得出韋恩圖(圖7和圖8)。

圖4 ATRA處理前(UNTREATED)、48h(ATRA48H)、72h(ATRA72H)后的野生型APL NB4細胞基因表達交集數(shù)據(jù)韋恩圖

圖5 與未用ATRA處理相比ATRA處理48h后基因分析火山圖，黑色代表非 DEGs,紅色代表上調(diào) DEGs,藍色代表下調(diào) DEG

圖6 與未用ATRA處理相比ATRA處理72h后基因分析火山圖，黑色代表非 DEGs,紅色代表上調(diào) DEGs,藍色代表下調(diào) DEG

圖7 ATRA分別處理48h和72h后上調(diào)最顯著的前10個 DEGs取交集

圖8 ATRA分別處理48h和72h后下調(diào)最顯著的前10個 DEGs取交集

表1 UNTREATED vs ATRA48H DEGs中上、下調(diào)最顯著的前10個基因具體信息

表2 UNTREATED vs ATRA72H DEGs中上、下調(diào)最顯著的前10個基因具體信息

3.3 DEGs的GO 和KEGG富 集 分 析

DEGs的GO富集分析結(jié)果用柱狀圖(圖9和10)表示如下，結(jié)果表明ATRA48H后上調(diào)DEGs在生物過程(biological process，BP)中主要與一氧化氮的生物合成過程、細胞周期調(diào)控、凋亡過程的負調(diào)控、葡萄糖進入細胞的負調(diào)控、異源代謝過程、RNA聚合酶II轉(zhuǎn)錄預起始復合物裝配、建立平面極性、宿主共生體生長的負調(diào)控、腎臟發(fā)育等過程相關(guān)；在細胞組成(cellular component，CC)層面主要與細胞纖維網(wǎng)、分泌顆粒、細胞的頂端部分等結(jié)構(gòu)相關(guān)；在分子功能(Molecular Function，MF)上主要與苦味受體活性、轉(zhuǎn)錄因子活性，RNA聚合酶Ⅱ核心啟動子序列特異性結(jié)合參與預啟動復合物組裝、碳酸鹽脫水酶活性、味覺受體活性等分子功能相關(guān)。ATRA48H后下調(diào)DEGs在BP中主要與對抗原刺激的急性炎癥反應、激活先天免疫反應、B細胞增殖的負調(diào)控、白三烯代謝過程、凋亡細胞清除、干擾素-γ介導和整合素的信號通路、中性粒細胞趨化途徑、細胞粘附、調(diào)節(jié)細胞形狀、白細胞遷移、細胞防御反應、細胞對干擾素-β的反應途徑、激活先天免疫反應、蛋白質(zhì)結(jié)合的正調(diào)控、經(jīng)由TAP依賴的MHC I類進行抗原加工和呈遞外源肽抗原、細胞表面受體信號通路、炎癥反應的正調(diào)控、血管生成的正調(diào)控、Wnt信號通路的正調(diào)控、蛋白質(zhì)瓜氨酸化等過程相關(guān)；在CC層面主要與細胞外泌體、胞漿、細胞膜與質(zhì)膜、蛋白酶體核心復合體、精原蛋白酶體復合物等結(jié)構(gòu)相關(guān)；在MF上主要與細胞粘附分子、鈣離子、整聯(lián)蛋白、花生四烯酸結(jié)合，蛋白質(zhì)精氨酸脫亞氨酶、蘇氨酸型內(nèi)肽酶活性等分子功能相關(guān)。用同樣的方法分析ATRA72H后上調(diào)和下調(diào)DEGs與ATRA48H差別不大，選取主要富集通路用柱狀圖(圖11和12)顯示。結(jié)果用KEGG可視化分析(篩選條件為校正后P＜0.05)結(jié)果表明,ATRA48H上調(diào)DEGs主要與流體剪切應力與動脈粥樣硬化條信號及代謝通路(hsa05418)密切相關(guān)，與下調(diào)DEGs相關(guān)的通路共篩選出15條，具體信息顯示在表3中，其中包括趨化因子信號通路(hsa04062)、細胞粘附分子(hsa04514)、NOD樣受體信號轉(zhuǎn)導通路(hsa04621)、白細胞跨內(nèi)皮遷移(hsa04670)、碳水化合物消化和吸收(hsa04973)、PI3K-Akt信號通路(hsa04151)、造血細胞譜系(hsa04640)、蛋白酶體(hsa03050)等。用同樣的方法分析ATRA72H后上調(diào)和下調(diào)DEGs與ATRA48H差別不大，結(jié)果顯示在表4中。

表3 UNTREATED vs ATRA48H 上、下調(diào)DEGs通過KEGG分析的主要通路

表4 UNTREATED vs ATRA72H 上、下調(diào)DEGs通過KEGG分析的主要通路

圖9 UNTREATED vs ATRA48H 上調(diào)DEGs通過GO富集的主要通路

圖10 UNTREATED vs ATRA48H 下調(diào)DEGs通過GO富集的主要通路

圖11 UNTREATED vs ATRA72H 上調(diào)DEGs通過GO富集的主要通路

3.4 PPI網(wǎng)絡(luò)構(gòu)建及分析

通過 STRING 數(shù)據(jù)庫構(gòu)建 PPI網(wǎng)絡(luò)(圖13和14)以及Cytoscape軟件中的CytoHubba插件進一步明確ATRA48H后上、下調(diào)DEGs所涉及的蛋白相互作用及其中參與調(diào)控的關(guān)鍵基因(圖15和16)，結(jié)果顯示，上調(diào)DEGs的PPI網(wǎng)絡(luò)中包括148個節(jié)點蛋白的145種關(guān)系，下調(diào)DEGs的PPI網(wǎng)絡(luò)中包括115個節(jié)點蛋白的241種關(guān)系，其中多數(shù)蛋白間都存在相互作用關(guān)系，篩選出下調(diào)DEGs排名前十的關(guān)鍵基因有：PSMB8、CXCR2、GNG2、CXCL3、PSMB9、CXCL6、P2RY13、TAS2R10、TAS2R7、PSMA3。用同樣方法得出ATRA72H后上、下調(diào)DEGs所涉及的蛋白相互作用圖(圖17和18)及其中參與調(diào)控的關(guān)鍵基因(圖19和20)，結(jié)果顯示，上調(diào)DEGs的PPI網(wǎng)絡(luò)中包括205個節(jié)點蛋白的261種關(guān)系，下調(diào)DEGs的PPI網(wǎng)絡(luò)中包括166個節(jié)點蛋白的508種關(guān)系，其中多數(shù)蛋白間都存在相互作用關(guān)系，篩選出下調(diào)DEGs排名前十的關(guān) 鍵 基 因 有：ITGAM、ITGB2、ITGAX、IRF1、OAS1、OAS3、IFI35、LILRB2、IRF2、IFI6。

圖12 UNTREATED vs ATRA72H 下調(diào)DEGs通過GO富集的主要通路

圖13 UNTREATED vs ATRA48H 上調(diào)DEGs的PPI網(wǎng)絡(luò)圖

圖14 UNTREATED vs ATRA48H 下調(diào)DEGs的PPI網(wǎng)絡(luò)圖

圖15 UNTREATED vs ATRA48H 上調(diào)DEGs的前10個關(guān)鍵基因

圖16 UNTREATED vs ATRA48H 下調(diào)DEGs的前10個關(guān)鍵基因

圖17 UNTREATED vs ATRA72H 上調(diào)DEGs的PPI網(wǎng)絡(luò)圖

圖18 UNTREATED vs ATRA72H 下調(diào)DEGs的PPI網(wǎng)絡(luò)圖

圖19 UNTREATED vs ATRA72H 上調(diào)DEGs的前10個關(guān)鍵基因

4 討論

本文分析了ATRA處理NB4細胞后的差異基因，GO分析結(jié)果表明上調(diào)DEGs主要與一氧化氮的生物合成過程、細胞周期調(diào)控、葡萄糖進入細胞的負調(diào)控、異源代謝過程、RNA聚合酶II轉(zhuǎn)錄預起始復合物裝配、建立平面極性、宿主共生體生長的負調(diào)控；而下調(diào)的DEGs主要與炎癥反應、中性粒細胞趨化途徑、細胞粘附、細胞凋亡、免疫反應、蛋白質(zhì)合成相關(guān)。同時在 KEGG 分析發(fā)現(xiàn)差異基因主要富集在趨化因子信號通路、細胞粘附分子、NOD樣受體信號轉(zhuǎn)導通路、白細胞跨內(nèi)皮遷移、PI3K-Akt信號傳導途徑。另外，在DEGs形成的PPI網(wǎng)絡(luò)進行分析，可觀察到下調(diào)基因參與的通路更加具有相關(guān)性，設(shè)計的通路和關(guān)系節(jié)點更加密集，對ATRA48H和ATRA72H后下調(diào)DEGs的PPI網(wǎng)絡(luò)中分別篩選出前10 位與ATRA處理后NB4細胞表達下調(diào)的關(guān)鍵基因分別繪制生存曲線(組圖23)，發(fā)現(xiàn)6個高表達提示預后不良，包括：OAS1、OAS3、PSMB9、PSMB8、ITGAM、IFI35可作為AML的危險因素，為AML的治療提供了新靶點。

圖20 UNTREATED vs ATRA72H 下調(diào)DEGs的前10個關(guān)鍵基因

圖21 KEGG趨化因子通路，紅色代表上調(diào)DEGs,藍色代表下調(diào)DEGs,

圖22 KEGG白細胞跨內(nèi)皮遷移通路，紅色代表上調(diào)DEGs,藍色代表下調(diào)DEGs

圖23 ATRA處理后NB4細胞表達下調(diào)的關(guān)鍵基因在AML中的生存曲線