基于全血基因表達譜的腦組織基因表達量預測模型的建立

徐文劍 李 巍

(國家兒童醫學中心 首都醫科大學附屬北京兒童醫院 遺傳與出生缺陷防治中心；北京市兒科研究所 出生缺陷遺傳學研究北京市重點實驗室；兒科重大疾病研究教育部重點實驗室，北京 100045)

基因表達實驗通過刻畫細胞轉錄水平變化情況來闡釋生物學表型，是輔助疾病機利解析、藥效評價和臨床疾病亞型分類等工作的有力工具[1-5]。關于轉錄組測序和基因芯片表達譜的基因表達綜合庫(gene expression omnibus，GEO)、基于網絡的細胞反應印記整合圖書館(library of integrated network-based cellular signatures，LINCS)和基因型-組織表達(genotype-tissue expression，GTEx)項目等公共數據庫和基因表達數據分析算法應運而生[6-9]。隨著實驗分析成本下降和組學數據分析流程完善，基于轉錄組測序的基因表達譜分析最終將從基礎研究走向臨床應用。

疾病通常有明確的受累組織，需要采集疾病相關的組織樣本，使用基于轉錄組測序或微陣列基因表達譜分析實驗以便精準辨識出疾病相關組織細胞內中基因異常表達情況。例如孤獨癥譜系障礙(autism spectrum disorder，ASD)的遺傳學分子機制可以通過基因表達分析來輔助解析，常用組織包括全血、淋巴母細胞系、小腦、額葉皮質、前額葉皮質、顳葉皮質、尾狀核等[10-12]。最近Quesnel-Vallières等[13]報道至少有12篇文獻采用前額葉皮質、顳上皮質和小腦等腦區組織樣本來進行基因表達分析來研究ASD遺傳基因作用機制和相關基因表達標志物。

腦組織樣品的基因表達實驗樣品收集工作難度大、風險高、成本高，限制了組織特異性基因表達譜分析在該領域的發展，因此亟需一種能規避腦組織取樣的表達譜檢測替代方案。全血樣本是一個無創、易采集的替代性組織，并且與腦組織有基因表達模式的相似性。例如Mundalil Vasu等[14]報道血清中的一組microRNA(miRNA)分子有可能成為ASD的非侵入性生物標志物。Sullivan等[15]報道了基因在血液與17個大腦組織的表達量成對斯皮爾曼(Spearman)相關系數中位數為0.5，血液與中樞神經系統組織有中等強度的相關性。Bosker等[16]報道了精神分裂癥(schizophrenia)動物模型在不給藥且不恐懼刺激條件下前扣帶回皮質和血液白細胞的基因表達量有弱相關性，在給藥并且恐懼刺激條件下有強相關性。Hensman等[17]利用亨廷頓癥(Huntington’s disease，HD)的兩個人群隊列中血液樣本，發現這兩組HD患者血液與公共數據庫中HD患者尾狀核的基因表達模式有顯著的一致性。

從血液組織表達數據推斷另一個組織中的基因表達量已被證明是可行的。Halloran等[18]曾報道利用血液與肺組織的基因表達的相關性建立過肺組織基因表達量預測模型，使用GTEx項目31對全血和肺組織樣本轉錄組數據為每個基因分別構建了以年齡、性別、該基因的血液表達量為輸入變量預測該基因的肺組織表達量的線性回歸模型，結果顯示18%的基因的血液表達量與肺組織表達量顯著相關。本研究立足于全血轉錄組基因表達量與腦組織中基因表達量的潛在相關性，挖掘了基因表達量跨組織多對多關系，構建了一個基于全血轉錄組表達量數據的未取樣腦組織中基因表達量的計算預測模型。

1 材料與方法

1.1 數據來源及處理

1)數據下載：截至2019年4月，GTEx研究聯盟收集并以高通量轉錄組測序(RNAseq)研究了來自714名生前健康的人類捐獻者的11 688份尸檢組織樣本，涵蓋53個組織，包括實體器官組織、腦分區、全血、兩種來自血液和皮膚的細胞系，這些樣本為研究基因表達的組織特異性和個體特異性提供了關鍵材料。該項目產生的所有分析數據保存在公共數據庫GTEx網站上，每位捐獻者的一個生物組織樣本對應著表達量數據集中的一條記錄。

從GTEx數據庫(https://gtexportal.org/)下載相關RNAseq基因原始表達量數據文件，主要包括2016-01-15_v7版本的每百萬讀序列的單位長度轉錄本(transcripts per kilobase million，TPM)表達量文件(GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz，841M)和GTEx數據庫樣本屬性的注釋信息(GTEx_v7_Annotations_SampleAttributesDS.txt，7.9M)。

2)數據分割：本研究期望建立的每個腦組織基因表達量預測的模型，要求同一個組織捐獻者同時有全血樣本和腦組織基因表達量數據。分別為每個腦組織從原始數據文件gct中提取出全血樣本與靶組織樣本的捐獻者編號(subject ID)完全匹配的數據條目，分別命名為靶組織基因表達量預測任務對應的全血表達量特征數據集(例如Brain-Cortex_features.txt)和靶組織基因表達量目標數據集(例如Brain-Cortex_targets.txt)，二者的數據條目一一對應，是待構建模型的實例數據。

3)數據清洗：在匹配全血表達量特征數據集中的每個基因表達量特征在部分捐獻者條目存在TPM表達量為 0的情況，不能區分是測序建庫的誤差還是樣本確實未表達該基因。為減小數據噪聲，剔除存在表達量TPM為0情況的特征，只保留所有捐獻者條目中表達量值完整的特征。同理，對靶組織基因目標數據集進行處理。各預測任務的特征維度在14 000左右，目標維度在19 000左右。

4)數據標準化：各組織表達量子數據集為二維矩陣，一個維度代表樣本，另一個維度代表基因。基因的表達量差異有統計學意義，表達量較高和較低的基因表達量數值的平均值可相差若干個數量級，無法直接用于回歸建模。因此使用python sklearn模塊中preprocessing.StandardScaler()函數依次在每個基因的維度上進行樣本表達量數據的Z-score標準化，得到均值為0，標準差為1的新的表達量值矩陣，使得所有基因的表達量達到同一個數量級且保留樣本間的差異信息。

1.2 基于線性模型的特征選擇

本研究構建基于最小二乘回歸(least squares regression, LSR)的線性回歸預測模型。本研究的樣本量區間為[50，101]，特征數均為104數量級。本研究采用經驗法則：要求特征維度必須小于樣本量，樣本量是特征數4倍及以上。彈性網絡(elastic net)是使用L1，L2范數作為正則項的線性模型，能學習出只有少量非零系數的稀疏模型，作為特征選擇的方法已報道[19-21]可應用到基于生物醫學文本的重癥監護病房(intensive care unit，ICU)的患者病情分級、基于基因表達譜和基因變異信息的藥物敏感性預測、腦區靜息態功能性磁共振成像[resting state functional magnetic resonance imaging，(rs)fMRI]等相關領域，因此本研究選用彈性網絡作為特征選擇方法。

預實驗階段先以基于彈性網絡模型的特征選擇方法分別為每個腦組織的前200個基因預測任務提出5、10、15和20個最相關的全血基因表達量特征，調用函數為sklearn模塊中的linear_model.ElasticNet(alpha=0.02, max_iter=10 000)。之后綜合考慮預測模型性能和可用樣本量的限制條件，從4種方案中選擇最佳關鍵特征個數完成剩余1萬多個基因預測任務的特征選擇。

以一個特定腦組織為例，對建模過程予以說明。設預處理后的全血表達量數據矩陣為Dblood，樣本量為S，共有M個基因的表達量信息，樣本i的全血表達譜數據記為向量xi=(xi(1),xi(2),…，xix(m))其中i∈(0,S),該腦組織表達譜數據矩陣為Dbrain，樣本量同為S，共有N個基因，樣本i的腦組織表達譜設為向量yi=(yi(1),yi(2),…，yix(N))其中i∈(0,S),針對以xi為自變量預測靶基因t的表達量yi(t)的基本預測任務，在特征選擇階段，選取彈性網絡模型為基因t的目標函數，求解最優參數

(1)

其中模型參數w(t)∈RM,b(t)∈R,正則項比例系數a∈(0,1)。

1.3 基于全血基因表達量的多組織基因表達量預測模型構建

(2)

1.4 統計學方法

1.4.1 預測模型性能評估指標

1.4.2 預測值與真實值的決定系數

2 結 果

2.1 腦組織與全血組織配對的可用樣本量統計

GTEx表達量數據集中包含13個腦組織，根據捐獻者樣本編號將全血樣本和腦組織樣本進行配對，經數據分割、清洗、標準化，分別將全血、每個腦組織的表達量數據提取為單獨數據文件。統計可用樣本量和可用基因個數，詳見表 1。

2.2 腦組織不同特征選擇方案下前200個基因表達量預測模型性能

各腦組織與全血配對可用的樣本量最小只有50。為確定特征選擇方案，將前200個基因預測任務作為預實驗。每個腦組織的前200個基因表達量預測任務分別提取5、10、15和20個最相關的全血基因表達量特征，基于全血基因表達量的LSR腦組織基因表達量預測模型的MAE結果匯總詳見表2。

表1 GTEx數據集中腦組織與全血樣本配對可用樣本統計信息Tab.1 Summary of available sample of brain tissue matched with blood in GTEx dataset

GTEx：genotype-tissue expression；BA:brodmann area; .

表2 各腦組織前200個基因表達量預測模型平均絕對誤差Tab.2 MAE of first 200 gene expression prediction task of each brain tissue

*: MAE of optimal feature sets;△: MAE of used feature sets；MAE:mean absolute error;BA:brodmann area.

對于杏仁核(amygdala)和黑質(substantia nigra)兩個樣本量最小(n=50，表 1)腦組織，提取25個最相關特征的預測模型達到最優預測性能，提取少于25個特征處在欠擬合區域，提取多于25個特征則處在過擬合區域；其他10個腦組織在30個最相關特征時達到最優性能(表 2中以*表示)，仍處于欠擬合區域。為保證腦組織最小樣本量是特征數的4倍左右，統一選擇提取15個最相關特征的方案(表2中以Δ表示)。重新用全部樣本訓練線性回歸模型，此線性模型得到的腦組織前200個基因表達量的預測值與真實值的擬合度良好(圖 1)。

2.3 腦組織全基因表達譜預測模型性能

根據預實驗結果，提取包含15個最相關的全血基因表達量特征構成低維度新特征數據集，并構建13個腦組織所有基因表達量線性回歸預測模型，預測模型在交叉驗證集上的性能總結詳見表3。所得的一系列線性回歸預測模型由式子(2)中最優參數表示。預測模型MAE和RMSE結果趨勢一致。其中杏仁核(amygdala)組織預測性能最佳，尾狀核(caudate (basal ganglia))組織預測性能最差。

圖1 13個腦組織前200個基因表達量預測值與真實值的線性擬合度Fig.1 Goodness of fit of model predicted value and true value of 200 genes in 13 tissues

A:cortex;B:cerebellum;C:hippocampus;D:substantia nigra;E:anterior cingulate cortex;F:frontal cortex;G:cerebellar hemisphere;H:caudate (basal ganglia);I:nucleus accumbens (basal ganglia);J:putamen (basal ganglia);K:hypothalamus;L:spinal cord (cervical c-1);M:amygdala.

表3 各腦組織基因表達量預測模型平均絕對誤差和均方根誤差Tab.3 MAE and RMSE of gene expression prediction model of brain tissues

3 討論

腦組織樣本手術取樣限制了相關組織內基因表達分析的大規模開展。本研究針對臨床研究中的腦組織樣本收集困難這一實際問題，從GTEx轉錄組數據集中挖掘腦組織轉錄組與全血轉錄組的內在數值相關性，期望摸索出一個利用有限樣本集構建具有一定準確度的轉錄組預測模型的通用設計模式。相比組織轉錄組的穩定性，全血轉錄組信息動態變化速度較快。受到GTEx數據集樣本量的規模所限制，從全血轉錄組中能夠間接推測出組織轉錄組真實信息的比例及線性預測模型可靠性依然是未知數。采用特征選擇方法對每個單元預測任務進行高維特征預處理，一是為提取出有利于提高預測性能的關鍵特征，二是起到過濾全血轉錄組中無關基因表達數據的干擾。

本研究應用常規線性回歸模型建立了基于全血轉錄組表達量數據的腦組織基因表達量的一組預測模型，表明僅用全血表達譜數據能比較準確地預測出未取樣腦組織基因表達量，將來進一步發展成熟后，或許可以在轉錄組研究中規避腦組織樣本的手術取樣。除了孤獨癥譜系障礙，腦組織基因表達量預測模型也能為帕金森病和精神分裂癥、雙相情感障礙、抑郁癥和酗酒等重要精神障礙疾病的基因表達譜研究提供一種備選工具[22-27]。

本研究的總體目標是針對每個靶組織獨立地構建一個基于全血樣本基因表達量的組織內基因表達量回歸預測模型，預測模型的輸入特征為個人全血樣本基因表達量，預測值為靶組織內各基因表達量。回歸預測模型性能評估環節計算預測表達量與真實表達量的差異程度。將一個靶組織預測模型構建任務化整為零，分解為靶組織所有基因表達量預測和靶組織的單基因表達量預測模型兩個層次。預測模型的最小單元任務為基于全血基因表達量數據的靶組織中的單基因表達量預測模型構建問題。

全血表達量數據包含約104個轉錄本，平均每個組織只有少于100個數據樣本，可知輸入特征的維度將高出樣本量2個數量級，原輸入特征直接用于回歸建模會帶來過擬合問題。因此本研究用特征選擇方法計算出每個輸入特征與預測目標的關聯程度表征數值，經排序選擇出適合訓練集樣本量規模和預測問題性質的一組最優的輸入特征。生物醫學數據集通常具有特征維度高和樣本量小的特點，阻礙了挖掘其中所蘊含的海量生物學信息和解碼深層次基因功能的基礎數據研究，特征提取日益成為生物醫學大數據發展和應用的一種專門技術[28]。本研究初步研究了應對這一問題的一種生物信息學分析思路，本預測模型為組織轉錄組分析建模工作提供了新的視野。