李建偉，岳欣蕾，胡和智

(河北工業(yè)大學(xué) 人工智能與數(shù)據(jù)科學(xué)學(xué)院，天津 300401)

急性T淋巴細(xì)胞白血病(t-cell acute lymphoblastic leukemia，T-ALL)是最常見的兒童惡性腫瘤疾病之一，占兒童惡性腫瘤的50%左右[1]。這種成年人罕見的血液系統(tǒng)惡性疾病，以T細(xì)胞異常增生、積聚和組織浸潤為特點，且危險等級高、復(fù)發(fā)率高[2]。臨床上，T-ALL的患者通常表現(xiàn)為貧血、造血功能衰竭、縱隔胸腺腫塊等，并伴有外周血的白細(xì)胞計數(shù)偏高、中性粒細(xì)胞減少和血小板減少[3]。該病發(fā)病機制復(fù)雜，易浸潤中樞神經(jīng)系統(tǒng)，死亡率相對較高，病因至今尚未完全清楚。T-ALL對化療敏感，雖緩解率高，但極易廣泛轉(zhuǎn)移和復(fù)發(fā)，預(yù)后效果不理想[4]。

近年來，隨著高通量測序技術(shù)的廣泛應(yīng)用和組學(xué)數(shù)據(jù)的指數(shù)型增長，研究者們從不同的組學(xué)層面對T-ALL進行了深入研究。有研究發(fā)現(xiàn)，表觀遺傳的改變會導(dǎo)致與細(xì)胞生長、凋亡及調(diào)控細(xì)胞周期相關(guān)的基因發(fā)生變化，并參與白血病的發(fā)病進程[5]。DNA甲基化作為表觀遺傳修飾的重要組成部分，其在基因的表達調(diào)控、基因組印記等生命活動中均發(fā)揮重要作用，與細(xì)胞的生長發(fā)育、基因的選擇性表達密切相關(guān)[6]。此外，CCCTC結(jié)合因子(CCCTC-binding factor，CTCF)通過結(jié)合特定的DNA序列實現(xiàn)對真核生物基因的轉(zhuǎn)錄調(diào)控[5]。DNA甲基化和轉(zhuǎn)錄因子CTCF結(jié)合，構(gòu)成了一個復(fù)雜的調(diào)控系統(tǒng)，其對T-ALL中的基因表達起著關(guān)鍵性的調(diào)控作用[7-8]。因此，對T-ALL中的DNA甲基化和轉(zhuǎn)錄因子CTCF調(diào)控功能進行研究具有重要意義，為深入闡釋T-ALL中的基因調(diào)控機制提供一種新的思路。

本研究整合了T-ALL的全基因組RNA-seq、CTCF ChIP-seq以及DNA甲基化數(shù)據(jù)，采用生物信息學(xué)方法對多組學(xué)數(shù)據(jù)進行綜合分析，篩選T-ALL患者與健康人T細(xì)胞之間的差異基因，并分別通過基因相似性融合網(wǎng)絡(luò)和PPI網(wǎng)絡(luò)篩選出其中的關(guān)鍵基因與核心基因。這些核心基因有成為T-ALL生物標(biāo)志物的潛力，為深入探索T-ALL的發(fā)病機制及研發(fā)靶向藥物提供幫助。

1 材料與方法

1.1 材料

GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)是由NCBI(美國國立生物技術(shù)信息中心)創(chuàng)建并維護的保存各種高通量實驗數(shù)據(jù)的公共存儲數(shù)據(jù)庫[9]。SRA數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/sra/)也由NCBI于2007年創(chuàng)建，主要用于存儲二代測序的原始序列數(shù)據(jù)[10-12]。研究從GEO數(shù)據(jù)庫中下載T-ALL的相關(guān)數(shù)據(jù)集，從GSE115895、GSE141140中得到RNA-seq數(shù)據(jù)，從GSE115893中得到CTCF ChIP-seq數(shù)據(jù)，從GSE42079中得到DNA methylation數(shù)據(jù)，并在SRA數(shù)據(jù)庫中下載其相應(yīng)的原始序列數(shù)據(jù)。T-ALL多組學(xué)數(shù)據(jù)分析流程如圖1所示。

(a)為不同組學(xué)的患者(T-all)與正常人(T-cell)樣本數(shù)目的比例[不同組學(xué)的數(shù)據(jù)在GEO與SRA數(shù)據(jù)庫的對應(yīng)關(guān)系為RNA-seq：SRR9822189-SRR9822203(GSE115895，T-all)與SRR10550198-SRR10550201(GSE141140，T-cell)，CTCF ChIP-seq：SRR9822138-SRR9822147(GSE115893，T-all)與SRR9822126-SRR9822131(GSE115893，T-cell)，DNA methylation：GSE42079]；(b)為RNA-seq、CTCF ChIP-seq和DNA methylation數(shù)據(jù)分析流程；(c)為不同組學(xué)整合后得到的共有差異基因；(d)為基因相似性網(wǎng)絡(luò)融合過程(其中1代表RNA-seq，2代表CTCF ChIP-seq，3代表DNA甲基化)。

1.2 數(shù)據(jù)預(yù)處理

參考基因組為人類基因組hg19，采用BWT算法對序列進行比對[13-14]。此外，使用htseq-count軟件[15]量化RNA-seq數(shù)據(jù)中基因的表達值，生成表達值矩陣；使用SAMtools軟件[16]對CTCF ChIP-seq數(shù)據(jù)進行PCR序列去重，利用MACS軟件[17]尋找可能的結(jié)合位點(即富集到基因組上的區(qū)域，又稱peaks區(qū)域)，應(yīng)用R語言的ChIPseeker包[18]提取peaks區(qū)域附近的基因；使用ChAMP工具包[19]對DNA甲基化芯片數(shù)據(jù)集進行原始數(shù)據(jù)處理，采用BMIQ算法進行表達值處理，如數(shù)據(jù)標(biāo)準(zhǔn)化、數(shù)據(jù)統(tǒng)計等。

1.3 差異基因篩選

為研究方便，本次采取較為寬松的閾值篩選差異基因，以期獲取更多的差異基因。對RNA-seq數(shù)據(jù)集分別應(yīng)用DESeq2、edgeR兩種統(tǒng)計模型進行統(tǒng)計分析[20]；使用R語言DiffBind包分析CTCFChIP-seq數(shù)據(jù)的差異結(jié)合位點，定義P<0.05為有效差異結(jié)合位點，規(guī)定peaks相對于轉(zhuǎn)錄起始位點(TSS)的距離在3 kb以內(nèi)，此類peaks內(nèi)的CTCF靶基因為有效基因[21-22]；調(diào)用limma函數(shù)對DNA甲基化數(shù)據(jù)集進行差異分析，通過BH方法進行多重假設(shè)檢驗的校正，當(dāng)P<0.01時為DNA甲基化差異探針，對其進行注釋，得到有效基因[19]。篩選3組數(shù)據(jù)集中共有的有效基因作為T-ALL的最終差異基因集，并進行后續(xù)研究分析。

1.4 構(gòu)建基因相似性融合網(wǎng)絡(luò)

RNA-seq相似性網(wǎng)絡(luò)、CTCF ChIP-seq相似性網(wǎng)絡(luò)和DNA甲基化相似性網(wǎng)絡(luò)是3個單一數(shù)據(jù)源的異構(gòu)相似性網(wǎng)絡(luò)，它們的特征信息各不相同。為更加準(zhǔn)確地篩選T-ALL關(guān)鍵基因，需對這些相似性網(wǎng)絡(luò)進行深度融合。研究選用相似性網(wǎng)絡(luò)融合方法(similarity network fusion，SNF)對以上3種相似性網(wǎng)絡(luò)進行融合[23-24]。SNF方法依次構(gòu)建單組學(xué)數(shù)據(jù)下的基因相似性網(wǎng)絡(luò)，經(jīng)迭代后得到多組學(xué)基因相似性融合網(wǎng)絡(luò)，最終篩選出關(guān)鍵基因。依據(jù)差異基因分析的結(jié)果，如log2(FC)、P值等這些特征值，構(gòu)建各組學(xué)的特征矩陣。首先，針對每個組學(xué)分析得到的差異基因，計算差異基因ri和rj的歐幾里得距離ρ(ri,rj)[25]。

(1)

(2)

通過式(2)計算差異基因間的相似性，得到相似性矩陣W，μ是一個超參數(shù)。矩陣ε用來消除數(shù)據(jù)的幅度差異性，矩陣ε中每個元素εi,j的計算公式如下所示：

(3)

其中，Mi表示基因ri的鄰居集合，mean(·)表示求d(ri,Mi)的均值；d(ri,Mi)表示距離向量，向量中的每個元素表示基因ri與其鄰居集Mi中每個基因間的距離。

對W(i,j)進行如下標(biāo)準(zhǔn)化：

(4)

同時計算基因ri和rj的局部親密度Y(i,j)：

(5)

當(dāng)僅有兩組數(shù)據(jù)時，基于迭代的思想，將不同組學(xué)得到的基因間相似矩陣進行融合：

(6)

(7)

當(dāng)組學(xué)數(shù)據(jù)類型為m(m>2)時，經(jīng)過迭代得到最終的相似性融合矩陣O(v)：

v=1,2,…,m(m>2)

(8)

(9)

A′k=D-1/2AkD1/2,Dii=∑jAijk

(10)

對每個基因，用F=[F[1],…,F[M]]獲得每個網(wǎng)絡(luò)中的等級，再對每個基因等級F[i]計算z-score值。最終對所有網(wǎng)絡(luò)中基因的z-score求平均值，獲得該基因的排名，并將其作為基因的重要性進行后續(xù)研究。若對100個基因進行排序，排名第一的基因?qū)?yīng)重要性的值為100，以此類推。

提取轉(zhuǎn)錄組的基因表達矩陣與DNA甲基化的矩陣，對基因進行共識聚類，分別確定轉(zhuǎn)錄組與DNA甲基化的最佳聚類數(shù)r與d，最終的聚類數(shù)k=r×d，并應(yīng)用以下公式所示的覆蓋度函數(shù)進行分析。

(11)

(12)

針對k類基因數(shù)據(jù)，統(tǒng)計在兩組學(xué)數(shù)據(jù)中的相同基因與所有基因。Sk為聚類k中的相同基因個數(shù)，Lk為聚類k中的所有基因個數(shù)，SL代表所有聚類相同基因的總個數(shù)。Vk為聚類k中的所有相同基因重要性總和，Lk為聚類k中的所有基因重要性總和，VI代表所有聚類相同基因重要性總和。最后通過SCORE=Mk+Nk來衡量聚類k，并取聚類k中的相同基因作為關(guān)鍵基因進行后續(xù)分析。

1.5 富集分析

為了確定關(guān)鍵基因富集的生物過程、細(xì)胞組分、分子功能以及生物途徑，使用在線網(wǎng)絡(luò)工具DAVID(Database for Annotation，Visualization and Integrated Discovery)進行GO(Gene Ontology，http://www.geneontology.org)功能富集分析和KEGG(Kyoto Encyclopedia of Genes and Genomes，https://www.kegg.jp/kegg/)通路富集分析[26-28]，篩選標(biāo)準(zhǔn)為P<0.05。

1.6 PPI網(wǎng)絡(luò)分析

利用STRING數(shù)據(jù)庫構(gòu)建關(guān)鍵基因的蛋白質(zhì)相互作用網(wǎng)絡(luò)，篩選標(biāo)準(zhǔn)為score>0.4[29]。應(yīng)用Cytoscape軟件可視化網(wǎng)絡(luò)，并利用MCODE插件篩選出score值最高的子網(wǎng)絡(luò)，網(wǎng)絡(luò)中的節(jié)點作為核心基因[30]。為了進一步驗證篩選得到的核心基因，通過癌癥基因普查數(shù)據(jù)庫CGC(Cancer Gene Census，https://cancer.sanger.ac.uk/census/)和比較毒理基因組學(xué)數(shù)據(jù)庫CTD(Comparative Toxicogenomics Database，http://ctdbase.org/)獲取與T-ALL相關(guān)的基因[31-32]。

2 結(jié)果與分析

2.1 差異基因的篩選

基于上述3種類型的組學(xué)數(shù)據(jù)，通過生物信息學(xué)方法進行分析，篩選出正常組和T-ALL疾病組之間差異基因進行后續(xù)分析。

2.1.1 RNA-seq差異基因

正常組和疾病組的RNA-seq數(shù)據(jù)經(jīng)過DESeq2處理后，獲得6 166個差異基因，其中上調(diào)4 461個，下調(diào)1 705個，這些差異基因聚類分析的結(jié)果如圖2所示。該RNA-seq數(shù)據(jù)經(jīng)過edgeR處理后，獲得7 790個差異基因，其中上調(diào)5 102個，下調(diào)2 688個，差異基因聚類分析結(jié)果如圖3所示。篩選標(biāo)準(zhǔn)均為|log2(FC)|>4且FDR<0.01。經(jīng)統(tǒng)計分析，兩種方法均包含的差異基因為5 887個，兩種方法得到的差異基因均能有效區(qū)分疾病組和正常組。從聚類結(jié)果中發(fā)現(xiàn)疾病組的第1與6、7組明顯不同于其他組，由此判斷可能存在T-ALL的亞型ETP-ALL(早期前體T細(xì)胞急性淋巴細(xì)胞白血病)。ETP細(xì)胞起源于造血干細(xì)胞，是由骨髓遷移到胸腺的細(xì)胞亞群，但因下載的數(shù)據(jù)中缺少相關(guān)臨床數(shù)據(jù)，未對免疫分型進行分析。在今后具備相關(guān)數(shù)據(jù)后，可對該方面進行深入分析。

圖2 DESeq2篩選差異基因的聚類分析

圖3 edgeR篩選差異基因的聚類分析

2.1.2 CTCF ChIP-seq差異基因

分別采用DiffBind中的DESeq2與edgeR兩種統(tǒng)計模型鑒定CTCF ChIP-seq數(shù)據(jù)的差異結(jié)合位點，基于結(jié)合親和力鑒定具有統(tǒng)計顯著性的差異結(jié)合位點。其中，DESeq2獲得差異結(jié)合位點25 664個，edgeR獲得差異結(jié)合位點19 095個。經(jīng)統(tǒng)計分析，兩種方法均包含的差異結(jié)合位點為11 607個，在peaks區(qū)域內(nèi)受CTCF調(diào)控的靶基因為5 315個。

2.1.3 DNA甲基化差異基因

對DNA甲基化數(shù)據(jù)分析后得到4 240個差異甲基化探針，注釋后共得到2 196個差異基因。對RNA-seq、CTCF ChIP-seq與DNA甲基化等3種組學(xué)數(shù)據(jù)得到差異基因集合取交集，得到119個共有差異基因，如圖4所示。

RNA1代表RNA-seq數(shù)據(jù)使用DESeq2方法分析；RNA2代表edgeR分析；ChIP1、ChIP2代表CTCF ChIP-seq數(shù)據(jù)使用DESeq2、edgeR方法分析；DNA代表DNA甲基化數(shù)據(jù)分析得到的差異基因。

2.2 基因相似性融合網(wǎng)絡(luò)與關(guān)鍵基因的篩選

利用2.1節(jié)得到的119個共有的差異基因構(gòu)建相似性融合網(wǎng)絡(luò)，采用DESeq2和edgeR分析RNA-seq和CTCF ChIP-seq數(shù)據(jù)，得到兩個融合網(wǎng)絡(luò)矩陣ODESeq2和OedgeR。整合兩個矩陣并對119個基因進行排序，排名第一的基因?qū)?yīng)重要性的值為119，以此類推?；蜻M行共識聚類結(jié)果顯示，轉(zhuǎn)錄組最佳聚類數(shù)為4，甲基化最佳聚類數(shù)為8。統(tǒng)計最終聚類中相同基因數(shù)目大于10的基因數(shù)據(jù)，得到4類不同的基因數(shù)據(jù)，通過覆蓋度函數(shù)計算得出相同基因數(shù)為48的數(shù)據(jù)SCORE最高。

2.3 關(guān)鍵基因的GO和KEGG分析

應(yīng)用DAVID對篩選得到的48個關(guān)鍵基因進行GO分析，結(jié)果顯示與生物學(xué)過程相關(guān)的有11個，主要包括生物調(diào)節(jié)、細(xì)胞過程、代謝過程以及對刺激的反應(yīng)等；與細(xì)胞成分相關(guān)的有5個，主要涉及細(xì)胞連接、細(xì)胞部分等；與分子功能有3個，主要涉及催化活性等，部分結(jié)果如表1所示。48個關(guān)鍵基因的KEGG分析結(jié)果顯示，主要參與FAS信號通路、代謝型谷氨酸受體II型通路、趨化因子和細(xì)胞因子信號傳導(dǎo)途徑介導(dǎo)的炎癥通路等。

表1 關(guān)鍵基因的Gene Ontology部分結(jié)果

2.4 構(gòu)建 PPI 網(wǎng)絡(luò)

利用STRING數(shù)據(jù)庫對48個關(guān)鍵基因構(gòu)建PPI網(wǎng)絡(luò)，通過Cytoscape的MCODE插件篩選出score分值最高(分值為6.286)的子網(wǎng)絡(luò)，如圖5所示。該網(wǎng)絡(luò)含有8個關(guān)鍵基因(CD7：CD7 molecule；GPR29：CCR6，C-C motif chemokine receptor 6；CTLA4：cytotoxic T-lymphocyte associated protein 4；CD5：CD5 molecule；CD274：CD274 molecule；IL2RB：interleukin 2 receptor subunit beta；FASLG：Fas ligand；CD247：CD247 molecule)，均作為核心基因。通過檢索CGC與CTD數(shù)據(jù)庫，表明篩選得到的8個核心基因確實與T-ALL有關(guān)，結(jié)果如表2所示。這8個核心基因所涉及的通路主要富集在造血細(xì)胞譜系、細(xì)胞因子受體相互作用等信號通路。

圖5 核心基因互作關(guān)系網(wǎng)絡(luò)

表2 CGC與CTD中8個核心基因的驗證結(jié)果

3 討論與結(jié)論

目前隨著生物數(shù)據(jù)大量涌現(xiàn)，有效地融合多組學(xué)數(shù)據(jù)構(gòu)建基因相互作用網(wǎng)絡(luò)來探索人類復(fù)雜疾病的致病機理具有重要的學(xué)術(shù)意義和廣泛的應(yīng)用價值[33]。近年來，在BELVER 等[2]的研究中發(fā)現(xiàn)，NOTCH信號通路和PI3K-Akt信號通路在T-ALL的發(fā)病機制中起著重要作用。NOTCH1可以通過合成代謝途徑(包括核糖體生物合成、蛋白質(zhì)翻譯以及核苷酸和氨基酸代謝)的轉(zhuǎn)錄上調(diào)促進白血病細(xì)胞的增殖和生存；PI3Kγ和PI3Kδ的活性對胸腺細(xì)胞的生長、增殖和存活也有直接調(diào)控作用，且T細(xì)胞的信號通路突變能激活PI3K-AKT-mTOR路徑。但是T-ALL等白血病的病因和發(fā)病機制目前尚未明確，它的發(fā)生和發(fā)展是多種因素共同作用的結(jié)果。已有研究表明，轉(zhuǎn)錄因子CTCF在T-ALL中發(fā)揮著關(guān)鍵的調(diào)控作用，其與基因結(jié)合的活性能夠被DNA甲基化所影響[34]。隨著生物信息學(xué)的快速發(fā)展，使得研究轉(zhuǎn)錄因子CTCF與DNA甲基化的關(guān)聯(lián)以及它們共同對T-ALL發(fā)生、發(fā)展的影響成為可能。

本文針對3種不同組學(xué)的T-ALL數(shù)據(jù)進行了綜合分析，篩選出共同差異基因及通路。為了使找到的差異基因更為準(zhǔn)確，研究對T-ALL的全基因組RNA-seq數(shù)據(jù)和CTCF ChIP-seq數(shù)據(jù)分別采用兩種基于統(tǒng)計學(xué)的生物信息學(xué)方法(DESeq2和edgeR)進行分析。兩種方法均使用負(fù)二項分布對讀段計數(shù)進行建模，但在離散度參數(shù)的選擇上有所不同。此外，運用CHAMP方法從DNA甲基化數(shù)據(jù)篩選差異基因。與各自正常的對照組相比，T-ALL中RNA-seq、CTCF ChIP-seq和DNA甲基化的差異基因數(shù)分別為5 887、5 315和2 196個。結(jié)果表明不同組學(xué)的差異基因數(shù)目差別很大，即在不同組學(xué)背景下，參與T-ALL發(fā)病的基因可能存在很大差異。當(dāng)然，這種改變有可能是由于數(shù)據(jù)較少、準(zhǔn)確性較低、抽樣誤差等因素引起。之后，對3組差異基因集取交集，用得到的共有差異基因構(gòu)建基因相似性融合網(wǎng)絡(luò)，從中篩選得到48個關(guān)鍵基因。這些關(guān)鍵基因主要對生物調(diào)節(jié)、免疫系統(tǒng)過程、代謝過程以及對刺激的反應(yīng)等生物學(xué)過程有一定影響，對其深入研究，有助于進一步加深對T-ALL分子機制的認(rèn)識。KEGG信號通路富集分析結(jié)果顯示，在FAS信號通路、代謝型谷氨酸受體II型通路及趨化因子和細(xì)胞因子信號傳導(dǎo)途徑介導(dǎo)的炎癥等信號通路有重要作用。研究表明Akt信號通路、T細(xì)胞受體信號傳導(dǎo)途徑與T-ALL的轉(zhuǎn)移、復(fù)發(fā)和預(yù)后相關(guān)，該通路有望成為T-ALL的治療靶點[35-36]。然而，其他信號通路在T-ALL發(fā)病過程中的具體作用機制目前很少報道，這些信號通路為T-ALL的分子機制研究提供了新的方向。研究分析參與FAS信號通路的基因有MAP3K5、FASLG、PARP3，它們能夠結(jié)合死亡受體TNFRSF6/FAS的細(xì)胞因子，在T-cell的發(fā)育中介導(dǎo)其由于細(xì)胞毒性引起的凋亡[37]。這些研究成果為研究T-ALL的發(fā)病機制及生物標(biāo)志物的篩選提供了理論依據(jù)。

為了進一步明確關(guān)鍵基因之間的相互作用，通過PPI分析篩選出CD7、GPR29、CTLA4、CD5、CD274、IL2RB、FASLG、CD247等8個核心基因。它們均為蛋白質(zhì)編碼基因，前4個基因在相關(guān)的研究中均已被證實。其中，CD7在T-ALL中高度表達，在T淋巴細(xì)胞發(fā)育成熟過程中與配體K12/SECTM1結(jié)合發(fā)揮協(xié)同刺激的作用，可作為治療T-ALL的有效靶標(biāo)[38]；T-ALL中參與趨化因子和細(xì)胞因子信號傳導(dǎo)途徑介導(dǎo)的炎癥相關(guān)信號通路被高度激活，GPR29、GNB3等炎癥關(guān)鍵基因高度表達，能夠吸引免疫細(xì)胞到達炎癥部位[39]；基因多態(tài)性導(dǎo)致的CTLA4異常表達與兒童急性T淋巴細(xì)胞白血病有關(guān)，CTLA4異常表達往往導(dǎo)致T細(xì)胞活化異常，從而影響機體免疫功能[40]；CD5在T-ALL中高表達，在T細(xì)胞受體信號轉(zhuǎn)導(dǎo)過程中起負(fù)調(diào)控作用，促進惡性T淋巴細(xì)胞的存活[41]。其他4個核心基因尚未有相關(guān)報道，提示相關(guān)實驗人員可對其進行深入研究。本文篩選出的關(guān)鍵基因為T-ALL診療提供了全新視角。研究的不足之處在于研究的原始數(shù)據(jù)來源于不同患者的組學(xué)數(shù)據(jù)，鑒于不同的人在基因遺傳背景上存在一定差異，研究取得的結(jié)果有待于在采集同一組T-ALL患者樣本上進一步驗證。