基于全血基因表達(dá)譜的腦組織基因表達(dá)量預(yù)測(cè)模型的建立

徐文劍 李 巍

(國(guó)家兒童醫(yī)學(xué)中心 首都醫(yī)科大學(xué)附屬北京兒童醫(yī)院 遺傳與出生缺陷防治中心；北京市兒科研究所 出生缺陷遺傳學(xué)研究北京市重點(diǎn)實(shí)驗(yàn)室；兒科重大疾病研究教育部重點(diǎn)實(shí)驗(yàn)室，北京 100045)

基因表達(dá)實(shí)驗(yàn)通過(guò)刻畫(huà)細(xì)胞轉(zhuǎn)錄水平變化情況來(lái)闡釋生物學(xué)表型，是輔助疾病機(jī)利解析、藥效評(píng)價(jià)和臨床疾病亞型分類(lèi)等工作的有力工具[1-5]。關(guān)于轉(zhuǎn)錄組測(cè)序和基因芯片表達(dá)譜的基因表達(dá)綜合庫(kù)(gene expression omnibus，GEO)、基于網(wǎng)絡(luò)的細(xì)胞反應(yīng)印記整合圖書(shū)館(library of integrated network-based cellular signatures，LINCS)和基因型-組織表達(dá)(genotype-tissue expression，GTEx)項(xiàng)目等公共數(shù)據(jù)庫(kù)和基因表達(dá)數(shù)據(jù)分析算法應(yīng)運(yùn)而生[6-9]。隨著實(shí)驗(yàn)分析成本下降和組學(xué)數(shù)據(jù)分析流程完善，基于轉(zhuǎn)錄組測(cè)序的基因表達(dá)譜分析最終將從基礎(chǔ)研究走向臨床應(yīng)用。

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

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

從血液組織表達(dá)數(shù)據(jù)推斷另一個(gè)組織中的基因表達(dá)量已被證明是可行的。Halloran等[18]曾報(bào)道利用血液與肺組織的基因表達(dá)的相關(guān)性建立過(guò)肺組織基因表達(dá)量預(yù)測(cè)模型，使用GTEx項(xiàng)目31對(duì)全血和肺組織樣本轉(zhuǎn)錄組數(shù)據(jù)為每個(gè)基因分別構(gòu)建了以年齡、性別、該基因的血液表達(dá)量為輸入變量預(yù)測(cè)該基因的肺組織表達(dá)量的線(xiàn)性回歸模型，結(jié)果顯示18%的基因的血液表達(dá)量與肺組織表達(dá)量顯著相關(guān)。本研究立足于全血轉(zhuǎn)錄組基因表達(dá)量與腦組織中基因表達(dá)量的潛在相關(guān)性，挖掘了基因表達(dá)量跨組織多對(duì)多關(guān)系，構(gòu)建了一個(gè)基于全血轉(zhuǎn)錄組表達(dá)量數(shù)據(jù)的未取樣腦組織中基因表達(dá)量的計(jì)算預(yù)測(cè)模型。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源及處理

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

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

2)數(shù)據(jù)分割：本研究期望建立的每個(gè)腦組織基因表達(dá)量預(yù)測(cè)的模型，要求同一個(gè)組織捐獻(xiàn)者同時(shí)有全血樣本和腦組織基因表達(dá)量數(shù)據(jù)。分別為每個(gè)腦組織從原始數(shù)據(jù)文件gct中提取出全血樣本與靶組織樣本的捐獻(xiàn)者編號(hào)(subject ID)完全匹配的數(shù)據(jù)條目，分別命名為靶組織基因表達(dá)量預(yù)測(cè)任務(wù)對(duì)應(yīng)的全血表達(dá)量特征數(shù)據(jù)集(例如Brain-Cortex_features.txt)和靶組織基因表達(dá)量目標(biāo)數(shù)據(jù)集(例如Brain-Cortex_targets.txt)，二者的數(shù)據(jù)條目一一對(duì)應(yīng)，是待構(gòu)建模型的實(shí)例數(shù)據(jù)。

3)數(shù)據(jù)清洗：在匹配全血表達(dá)量特征數(shù)據(jù)集中的每個(gè)基因表達(dá)量特征在部分捐獻(xiàn)者條目存在TPM表達(dá)量為 0的情況，不能區(qū)分是測(cè)序建庫(kù)的誤差還是樣本確實(shí)未表達(dá)該基因。為減小數(shù)據(jù)噪聲，剔除存在表達(dá)量TPM為0情況的特征，只保留所有捐獻(xiàn)者條目中表達(dá)量值完整的特征。同理，對(duì)靶組織基因目標(biāo)數(shù)據(jù)集進(jìn)行處理。各預(yù)測(cè)任務(wù)的特征維度在14 000左右，目標(biāo)維度在19 000左右。

4)數(shù)據(jù)標(biāo)準(zhǔn)化：各組織表達(dá)量子數(shù)據(jù)集為二維矩陣，一個(gè)維度代表樣本，另一個(gè)維度代表基因?；虻谋磉_(dá)量差異有統(tǒng)計(jì)學(xué)意義，表達(dá)量較高和較低的基因表達(dá)量數(shù)值的平均值可相差若干個(gè)數(shù)量級(jí)，無(wú)法直接用于回歸建模。因此使用python sklearn模塊中preprocessing.StandardScaler()函數(shù)依次在每個(gè)基因的維度上進(jìn)行樣本表達(dá)量數(shù)據(jù)的Z-score標(biāo)準(zhǔn)化，得到均值為0，標(biāo)準(zhǔn)差為1的新的表達(dá)量值矩陣，使得所有基因的表達(dá)量達(dá)到同一個(gè)數(shù)量級(jí)且保留樣本間的差異信息。

1.2 基于線(xiàn)性模型的特征選擇

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

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

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

(1)

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

1.3 基于全血基因表達(dá)量的多組織基因表達(dá)量預(yù)測(cè)模型構(gòu)建

(2)

1.4 統(tǒng)計(jì)學(xué)方法

1.4.1 預(yù)測(cè)模型性能評(píng)估指標(biāo)

1.4.2 預(yù)測(cè)值與真實(shí)值的決定系數(shù)

2 結(jié) 果

2.1 腦組織與全血組織配對(duì)的可用樣本量統(tǒng)計(jì)

GTEx表達(dá)量數(shù)據(jù)集中包含13個(gè)腦組織，根據(jù)捐獻(xiàn)者樣本編號(hào)將全血樣本和腦組織樣本進(jìn)行配對(duì)，經(jīng)數(shù)據(jù)分割、清洗、標(biāo)準(zhǔn)化，分別將全血、每個(gè)腦組織的表達(dá)量數(shù)據(jù)提取為單獨(dú)數(shù)據(jù)文件。統(tǒng)計(jì)可用樣本量和可用基因個(gè)數(shù)，詳見(jiàn)表 1。

2.2 腦組織不同特征選擇方案下前200個(gè)基因表達(dá)量預(yù)測(cè)模型性能

各腦組織與全血配對(duì)可用的樣本量最小只有50。為確定特征選擇方案，將前200個(gè)基因預(yù)測(cè)任務(wù)作為預(yù)實(shí)驗(yàn)。每個(gè)腦組織的前200個(gè)基因表達(dá)量預(yù)測(cè)任務(wù)分別提取5、10、15和20個(gè)最相關(guān)的全血基因表達(dá)量特征，基于全血基因表達(dá)量的LSR腦組織基因表達(dá)量預(yù)測(cè)模型的MAE結(jié)果匯總詳見(jiàn)表2。

表1 GTEx數(shù)據(jù)集中腦組織與全血樣本配對(duì)可用樣本統(tǒng)計(jì)信息Tab.1 Summary of available sample of brain tissue matched with blood in GTEx dataset

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

表2 各腦組織前200個(gè)基因表達(dá)量預(yù)測(cè)模型平均絕對(duì)誤差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.

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

2.3 腦組織全基因表達(dá)譜預(yù)測(cè)模型性能

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

圖1 13個(gè)腦組織前200個(gè)基因表達(dá)量預(yù)測(cè)值與真實(shí)值的線(xiàn)性擬合度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 各腦組織基因表達(dá)量預(yù)測(cè)模型平均絕對(duì)誤差和均方根誤差Tab.3 MAE and RMSE of gene expression prediction model of brain tissues

3 討論

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

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

本研究的總體目標(biāo)是針對(duì)每個(gè)靶組織獨(dú)立地構(gòu)建一個(gè)基于全血樣本基因表達(dá)量的組織內(nèi)基因表達(dá)量回歸預(yù)測(cè)模型，預(yù)測(cè)模型的輸入特征為個(gè)人全血樣本基因表達(dá)量，預(yù)測(cè)值為靶組織內(nèi)各基因表達(dá)量?；貧w預(yù)測(cè)模型性能評(píng)估環(huán)節(jié)計(jì)算預(yù)測(cè)表達(dá)量與真實(shí)表達(dá)量的差異程度。將一個(gè)靶組織預(yù)測(cè)模型構(gòu)建任務(wù)化整為零，分解為靶組織所有基因表達(dá)量預(yù)測(cè)和靶組織的單基因表達(dá)量預(yù)測(cè)模型兩個(gè)層次。預(yù)測(cè)模型的最小單元任務(wù)為基于全血基因表達(dá)量數(shù)據(jù)的靶組織中的單基因表達(dá)量預(yù)測(cè)模型構(gòu)建問(wèn)題。

全血表達(dá)量數(shù)據(jù)包含約104個(gè)轉(zhuǎn)錄本，平均每個(gè)組織只有少于100個(gè)數(shù)據(jù)樣本，可知輸入特征的維度將高出樣本量2個(gè)數(shù)量級(jí)，原輸入特征直接用于回歸建模會(huì)帶來(lái)過(guò)擬合問(wèn)題。因此本研究用特征選擇方法計(jì)算出每個(gè)輸入特征與預(yù)測(cè)目標(biāo)的關(guān)聯(lián)程度表征數(shù)值，經(jīng)排序選擇出適合訓(xùn)練集樣本量規(guī)模和預(yù)測(cè)問(wèn)題性質(zhì)的一組最優(yōu)的輸入特征。生物醫(yī)學(xué)數(shù)據(jù)集通常具有特征維度高和樣本量小的特點(diǎn)，阻礙了挖掘其中所蘊(yùn)含的海量生物學(xué)信息和解碼深層次基因功能的基礎(chǔ)數(shù)據(jù)研究，特征提取日益成為生物醫(yī)學(xué)大數(shù)據(jù)發(fā)展和應(yīng)用的一種專(zhuān)門(mén)技術(shù)[28]。本研究初步研究了應(yīng)對(duì)這一問(wèn)題的一種生物信息學(xué)分析思路，本預(yù)測(cè)模型為組織轉(zhuǎn)錄組分析建模工作提供了新的視野。