李建偉，岳欣蕾，胡和智

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

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

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

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

1 材料與方法

1.1 材料

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

(a)為不同組學的患者(T-all)與正常人(T-cell)樣本數目的比例[不同組學的數據在GEO與SRA數據庫的對應關系為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數據分析流程；(c)為不同組學整合后得到的共有差異基因；(d)為基因相似性網絡融合過程(其中1代表RNA-seq，2代表CTCF ChIP-seq，3代表DNA甲基化)。

1.2 數據預處理

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

1.3 差異基因篩選

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

1.4 構建基因相似性融合網絡

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

(1)

(2)

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

(3)

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

對W(i,j)進行如下標準化：

(4)

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

(5)

當僅有兩組數據時，基于迭代的思想，將不同組學得到的基因間相似矩陣進行融合：

(6)

(7)

當組學數據類型為m(m>2)時，經過迭代得到最終的相似性融合矩陣O(v)：

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

(8)

(9)

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

(10)

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

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

(11)

(12)

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

1.5 富集分析

為了確定關鍵基因富集的生物過程、細胞組分、分子功能以及生物途徑，使用在線網絡工具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]，篩選標準為P<0.05。

1.6 PPI網絡分析

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

2 結果與分析

2.1 差異基因的篩選

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

2.1.1 RNA-seq差異基因

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

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

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

2.1.2 CTCF ChIP-seq差異基因

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

2.1.3 DNA甲基化差異基因

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

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

2.2 基因相似性融合網絡與關鍵基因的篩選

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

2.3 關鍵基因的GO和KEGG分析

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

表1 關鍵基因的Gene Ontology部分結果

2.4 構建 PPI 網絡

利用STRING數據庫對48個關鍵基因構建PPI網絡，通過Cytoscape的MCODE插件篩選出score分值最高(分值為6.286)的子網絡，如圖5所示。該網絡含有8個關鍵基因(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數據庫，表明篩選得到的8個核心基因確實與T-ALL有關，結果如表2所示。這8個核心基因所涉及的通路主要富集在造血細胞譜系、細胞因子受體相互作用等信號通路。

圖5 核心基因互作關系網絡

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

3 討論與結論

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

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

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