基因大數(shù)據(jù)的集成分析

胡湘紅¹, 彭衡², 楊燦³, 張縱輝¹, 萬翔¹, 羅智泉¹

1 深圳市大數(shù)據(jù)研究院，廣東 深圳 518172

2 香港浸會大學數(shù)學系，香港 999077

3 香港科技大學數(shù)學系，香港 999077

摘要：隨著生物科技（如基因芯片和測序技術(shù)）的飛速發(fā)展，全世界已經(jīng)積累了海量的數(shù)據(jù)。有效地整合和集成多層面和多維度的基因大數(shù)據(jù)，對于全方位解析從遺傳變異到疾病發(fā)生的整個因果鏈條具有關(guān)鍵作用，可為個性化、精準醫(yī)療服務(wù)奠定科學的基礎(chǔ)。從3個方面對基因大數(shù)據(jù)的集成分析進行綜述：檢測風險位點及其功能分析、基因多效性的分析、基于孟德爾隨機化的因果推斷。進一步結(jié)合具體的應(yīng)用案例進行了闡述，最后對基因大數(shù)據(jù)的集成分析研究進行了總結(jié)以及展望。

關(guān)鍵詞：GWAS ; 集成分析 ; 多基因效應(yīng) ; 基因多效性 ; 孟德爾隨機化

論文引用格式：

胡湘紅, 彭衡, 楊燦, 張縱輝, 萬翔, 羅智泉.基因大數(shù)據(jù)的集成分析. 大數(shù)據(jù)[J], 2019, 5(4):67-88

HU X H, PENG H, YANG C, ZHANG Z H, WAN X, LUO Z Q.Integrative analysis for big data in genomics. Big Data Research[J], 2019, 5(4):67-88

1 引言

人類基因組計劃（human genome project，HGP）以及人類遺傳變異圖譜在21世紀初宣告完成。這一里程碑式的事件拉開了大規(guī)模利用生物醫(yī)療數(shù)據(jù)研究復(fù)雜人類疾病的序幕。隨著生物科技（如基因芯片和測序技術(shù)）的飛速發(fā)展，全世界的研究者已經(jīng)積累了多層面、多維度的基因大數(shù)據(jù)。這些數(shù)據(jù)覆蓋從遺傳變異到生命體表征的各個層面的數(shù)據(jù)，包括基因組（genome）層面、表觀基因組（epigenome）層面、轉(zhuǎn)錄組（transcriptome）層面、蛋白質(zhì)組（proteome）層面、代謝組（metabolome）層 面 以 及 生 物 體 層 面 的 表 型 特 征（phenome），成為科學家研究復(fù)雜疾病的寶貴資源。無疑，基因大數(shù)據(jù)將對現(xiàn)有醫(yī)學研究、個體化醫(yī)療產(chǎn)生顛覆性的影響。

全基因組關(guān)聯(lián)分析（genome-wide association studies，GWAS）在基因大數(shù)據(jù)的研究中處于基礎(chǔ)性地位。以GWAS為基礎(chǔ)的大量研究項目的開展，開啟了人類對各種復(fù)雜疾病的遺傳結(jié)構(gòu)的探索。GWAS采用高通量的方式獲取全基因組的上百萬個遺傳變異位點——單核苷酸變異（single-nucleotide polymorphsim， SNP），并對這些變異點和復(fù)雜疾病或性狀進行關(guān)聯(lián)分析，找出風險變異點。截至2019年4月，已經(jīng)公開發(fā)表的GWAS研究達3 923個，發(fā)現(xiàn)了約134 705個與疾病或性狀顯著關(guān)聯(lián)的變異位點（顯著性檢驗的p值<5×10^-8）。大部分的研究結(jié)果可在全基因組關(guān)聯(lián)分析數(shù)據(jù)庫中獲得。這些結(jié)果將基因組和生物體表型組聯(lián)系起來，極大地增進了人們對生物體表型的遺傳學結(jié)構(gòu)的認識。隨著GWAS研究越來越受到重視，大規(guī)模的研究也在不斷推進。例如，英國生物樣本庫（UK Biobank）成立于2006年，迄今為止，UK Biobank收集了約50萬人的DNA測序數(shù)據(jù)以及大量與疾病和健康相關(guān)的表征數(shù)據(jù)，包括生理指標、生活方式、血液和尿液采樣、大腦和身體成像數(shù)據(jù)以及參與者的醫(yī)療記錄、電子病歷數(shù)據(jù)等。這種從基因測序數(shù)據(jù)、表征數(shù)據(jù)到醫(yī)療成像數(shù)據(jù)和電子病歷數(shù)據(jù)的多維度積累為科學家進行復(fù)雜疾病的研究提供了數(shù)據(jù)支撐，將助力于個性化的精準醫(yī)療服務(wù)。

盡管GWAS的研究取得了令人矚目的成績，然而單獨的GWAS研究往往只考慮基因組層面與表型特征的信息，忽略了生命過程中其他層面的相關(guān)信息，無法貫穿從遺傳變異到表型特征的整條證據(jù)鏈。因此需要整合和集成從遺傳變異到生命體表征的多組學層面的信息并進行分析，從而全方位解析從遺傳變異到表型特征的整個因果鏈條。這里的多組學層面包括基因組層面、表觀基因組層面、轉(zhuǎn)錄組層面、蛋白質(zhì)組層面、代謝組層面以及生命體表征組層面的表型特征（如圖1所示）。伴隨著各項大型生物醫(yī)療國際項目的開展，這些多層面、多維度、高質(zhì)量的基因大數(shù)據(jù)正在不斷產(chǎn)生和積累。2012年9月完成的DNA元件百科全書項目（encyclopedia of DNA elements project，ENCODE項目）是繼人類基因組計劃后又一重要的突破性工程。該項目對大部分非編碼序列（約占全基因組的98%）的功能進行了注釋，例如“這部分序列與一種蛋白質(zhì)結(jié)合”“這部分序列常被甲基基團標記”“這部分序列通常隱藏在組蛋白的包圍中”。隨之開展的表觀基因組的圖譜計劃（roadmap epigenomics project）直接從人體的細胞系或組織中取樣，并將其數(shù)據(jù)向公眾開放。這些結(jié)果將為研究人類細胞系和組織的表觀基因組功能發(fā)揮重要的作用。與此同時，2010年開始的基因型-組織表達（genotype-tissue expression， GTEX）項目于2013公開數(shù)據(jù)庫，到目前為止收集了約714個捐獻者53個人體組織的11 688個樣本數(shù)據(jù)，用于研究不同組織中基因型與表達型的關(guān)系，完成了基因組到轉(zhuǎn)錄組的跨越，為科學家提供了寶貴的資源庫。此外，現(xiàn)代核磁共振技術(shù)、質(zhì)譜和色譜等技術(shù)的發(fā)展使得對大規(guī)模的生物體小分子的定量研究成為可能，這也大大促進了蛋白質(zhì)組層面和代謝組層面研究的發(fā)展。目前已公布的蛋白質(zhì)層面和代謝組層面的研究達80多項，產(chǎn)生了大量可公開獲取的數(shù)據(jù)。

圖1 多組學層面的數(shù)據(jù)

這些多層面、多維度、高質(zhì)量的基因大數(shù)據(jù)為刻畫完整的疾病產(chǎn)生的因果鏈條創(chuàng)造了條件，同時也促進了該領(lǐng)域集成分析方法研究的發(fā)展。本文總結(jié)了當前基因大數(shù)據(jù)的集成分析的一些研究進展，這些研究大部分基于多基因遺傳結(jié)構(gòu)（多個基因影響一個性狀）的共識以及以下兩方面的科學發(fā)現(xiàn)：遺傳變異的多效性（一個變異影響多種表型）以及非編碼遺傳變異的調(diào)控功能。具體地，本文對基因大數(shù)據(jù)的集成分析從以下3個方面進行綜述：檢測風險位點及其功能分析；基因多效性的分析；基于孟德爾隨機化的因果推斷。本文進一步結(jié)合具體的應(yīng)用案例進行分析，最后對基因大數(shù)據(jù)的集成分析研究進行了總結(jié)以及展望。

2 基因大數(shù)據(jù)集成分析的科學依據(jù)

2.1 多基因遺傳結(jié)構(gòu)

復(fù)雜性狀或疾病的多基因遺傳結(jié)構(gòu)是指復(fù)雜性狀和疾病的變異是由多個基因共同影響的。盡管科學家在GWAS研究中發(fā)現(xiàn)了很多顯著的變異位點，但是這些顯著的變異位點只能解釋性狀方差的小部分，這個現(xiàn)象通常被稱為“丟失的遺傳率（missing heritability）”。以人體身高這一性狀為例，根據(jù)以往對家庭譜系的研究，已經(jīng)存在的廣泛認識是人類身高的遺傳率達70%~80%等。多基因遺傳結(jié)構(gòu)的存在使得GWAS研究面臨著諸多挑戰(zhàn)：第一，多基因結(jié)構(gòu)意味著單個基因?qū)π誀畹呢暙I是微弱的，受GWAS樣本量的限制，單個GWAS研究很難找到遺傳變異中微弱的信號點；第二，常用的基于稀疏和強信號假設(shè)的建模方法在這里不再適用。圖2中橫軸為染色體編號，每個SNP按染色體上的位置排列，縱軸為對應(yīng)SNP的GWAS統(tǒng)計檢驗的-lg(p)。圖中頂部的點對應(yīng)-lg(p)>30的SNP。

圖2 身高的曼哈頓圖

基于多基因遺傳結(jié)構(gòu)的假設(shè)，線性混合模型（LMM）開始受到關(guān)注。以遺傳率的估算為例，假設(shè)已經(jīng)獲取n個個體的M個SNP的標準化后的基因測序數(shù)據(jù)G∈R^n×M，對應(yīng)的表型（如身高）數(shù)據(jù)y∈R^n×1，X∈R^n×p表示 p個協(xié)變 量信息（如年齡、性別和一些用以控制群體分層的主成分數(shù)據(jù)），它們之間的關(guān)系可以由線性混合模型來建立，具體如下：

其中，β∈R^p×1是固定效應(yīng)，u∈R^M×1是隨機效應(yīng)，e是由環(huán)境因素引起的隨機誤差。這里假設(shè)u和e都是正態(tài)分布的，即和。式（1）建立起了基因型與表型的關(guān)系。由遺傳率的定義（基因型方差在表型方差中的占比）及以上線性混合模型，可以得到遺傳率的計算式為：

限制極大似然估計（REML）和最大期望算法（expectation-maximization algorithm，EM）常被用于估算遺傳率，如全基因組復(fù)雜性狀分析（genome-wide complex trait analysis，GCTA）工具。

2.2 基因的多效性

基因的多效性（pleiotropy）是指一個基因影響著多種表型。基因的多效性廣泛存在于復(fù)雜性狀中，具體的例子有：與維生素D缺乏癥強相關(guān)的CYP2R1基因會同時導致多發(fā)性硬化癥；CLPTM1基因同時影響神經(jīng)膠質(zhì)瘤、膀胱癌和肺癌；基因多效性大量存在于自身免疫綜合征和精神疾病中。根據(jù)遺傳變異影響性狀的機制，基因多效性主要可分為兩大類：一類是基因的生物多效性（biological pleiotropy），即基因直接影響多種表型；另一類是基因的中介多效性（mediated pleiotropy），即基因通過影響一種表型對另一種表型產(chǎn)生影響（如圖3所示）。以肥胖基因（FTO基因）為例，研究者發(fā)現(xiàn)，該基因?qū)顷P(guān)節(jié)炎的影響是通過影響體質(zhì)指數(shù)（body mass index，BMI）來傳遞的?；虻闹薪槎嘈栽谔剿鲝?fù)雜疾病或性狀之間的因果關(guān)系中有重要的作用。

圖3 基因多效性

基因多效性的存在引導人們集成多個GWAS數(shù)據(jù)一起研究。一方面，基因的多效性可以幫助人們探索復(fù)雜疾病之間的關(guān)系，如遺傳相關(guān)性、基于孟德爾隨機化的復(fù)雜疾病的因果推斷等；另一方面，由于多基因遺傳結(jié)構(gòu)的存在，單獨分析一個GWAS會導致具備微弱效應(yīng)的遺傳位點很難被檢測到。集成多個與遺傳相關(guān)的復(fù)雜性狀的GWA S數(shù)據(jù)，通過建立有效的統(tǒng)計模型，可以幫助檢測到更多的多效性位點，從而提高統(tǒng)計分析的效率。

2.3 基因的調(diào)控功能

作為基本的遺傳單位，基因是一段有功能性的DNA序列?；蛑写蟛糠治稽c（如啟動子和增強子）處于非編碼區(qū)，不參與基因編碼，但對基因的轉(zhuǎn)錄、翻譯和表達起著重要的調(diào)控作用。對于整體的DNA序列，只有少部分位于基因的調(diào)控區(qū)，參與基因編碼，能夠轉(zhuǎn)錄為信使RNA，進而指導蛋白質(zhì)的合成。ENCODE項目對人類基因序列進行了詳細分析，研究結(jié)果顯示，在98%人類非編碼區(qū)的序列中，81%保持著生物活性。大量證據(jù)顯示，在GWAS發(fā)現(xiàn)的與眾多復(fù)雜性狀相關(guān)的位點中，大部分位點處在基因的非編碼區(qū)。比如，科學家調(diào)查5 654個非編碼區(qū)位點（這些位點與654種復(fù)雜疾病相關(guān)聯(lián)），發(fā)現(xiàn)76.6%的位點處于脫氧核糖核酸酶I高敏區(qū)（DHS），這一區(qū)域與基因的轉(zhuǎn)錄功能高度相關(guān)?；谵D(zhuǎn)錄組層面的研究發(fā)現(xiàn)，大量與疾病關(guān)聯(lián)的位點位于表達數(shù)量性狀基因座（expression quantitative trait loci， eQTL）（如圖4所示）以及轉(zhuǎn)錄因子結(jié)合位點（transcription factor binding sites，TF）。在最新的一項關(guān)于精神分裂癥（schizophrenia，SCZ）的研究中發(fā)現(xiàn)，在與SCZ相關(guān)的100多個基因位點中，有超過20%的位點與基因的表達和調(diào)控機制相關(guān)。這些研究結(jié)果表明，與疾病關(guān)聯(lián)的風險變異點可能不直接參與基因編碼，而是通過調(diào)控基因的轉(zhuǎn)錄、表達或翻譯中的任何一個環(huán)節(jié)發(fā)揮作用。因此將疾病的GWAS數(shù)據(jù)和功能型數(shù)據(jù)集成并進行分析，能夠深化人們對疾病的遺傳機理的認識。

圖4 與疾病關(guān)聯(lián)的位點常富集于表達數(shù)量性狀基因座

3 基因大數(shù)據(jù)的集成分析方法

基于上述科學進展，本節(jié)將從以下3個方面對基因大數(shù)據(jù)的集成分析進行綜述。

● 檢測風險位點及其功能分析：該類型方法集成GWAS數(shù)據(jù)和多組學數(shù)據(jù)（如變異位點功能型數(shù)據(jù)、轉(zhuǎn)錄組的基因表達數(shù)據(jù)等），力圖找到更多的危險變異點，并對變異點的生理功能進行分析。

● 基因多效性的分析：該類方法集成多個GWAS數(shù)據(jù)以及多組學數(shù)據(jù)，從而探索復(fù)雜疾病遺傳上的相關(guān)關(guān)系以及共同的致病基因。

● 基于孟德爾隨機化的因果推斷：該類型方法集成多種風險因素和疾病的GWAS數(shù)據(jù)，探討風險因素（如血壓、體重、代謝物等）對疾病的影響。

為描述方便，先對使用的數(shù)據(jù)類型進行介紹。在基因大數(shù)據(jù)的集成分析中使用的數(shù)據(jù)主要包括兩種類型：第一類是個體樣本層面的數(shù)據(jù)，該類型數(shù)據(jù)包含每個樣本的基因型、每個樣本的疾病狀態(tài)（糖尿病、高血壓或者健康）等；第二類是概括型數(shù)據(jù)，該類型數(shù)據(jù)是指對于個體樣本層面的數(shù)據(jù)通過關(guān)聯(lián)分析得到每一個SNP與疾病的概括性統(tǒng)計量，包括回歸效應(yīng)的估計值、標準差、檢驗統(tǒng)計量、p值等。因為使用的數(shù)據(jù)類型不同，集成分析方法也會存在差異。

3.1 風險遺傳變異的檢測及功能

由于與疾病相關(guān)的位點總是具備某些調(diào)控功能，將相關(guān)位點的調(diào)控功能信息植入某種疾病的GWAS數(shù)據(jù)分析中，能幫助人們找到真正的風險遺傳變異點，深入地研究遺傳變異的功能。在這個研究領(lǐng)域里，貝葉斯方法被大量運用到GWAS數(shù)據(jù)與調(diào)控信息數(shù)據(jù)的集成分析中。該類方法對疾病的GWAS個體層面數(shù)據(jù)或概括性數(shù)據(jù)（用G表示）進行分析，同時引入功能型數(shù)據(jù)或基因表達數(shù)據(jù)（用A表示），推斷遺傳變異位點L（可以是一個基因也可以是單個SNP）是否為疾病發(fā)生的風險位點（見計算式（3））。基于模型分析的結(jié)果可以進一步對變異點進行功能分析。在這里，筆者介紹兩類風險遺傳變異檢測的集成分析方法：第一類是結(jié)合功能型數(shù)據(jù)的集成分析；第二類是結(jié)合基因表達數(shù)據(jù)的集成分析。

P(L是危險變異點│G WA S數(shù)據(jù)G，功能型數(shù)據(jù)A) （3）

3.1.1 結(jié)合功能型數(shù)據(jù)的集成分析

假設(shè)已經(jīng)收集了M個位點的D列功能型數(shù)據(jù)，并存放在矩陣A中，A中對應(yīng)的第i行、第j列的元素記為A_{i j}，代表第i個SNP的第j種功能型的取值。以eQTL數(shù)據(jù)為例，假設(shè)A中的每一列對應(yīng)人體的不同組織器官（如心臟、大腦、肝臟等）。如果第i個SNP能夠調(diào)節(jié)第j個組織中某個基因的表達，那么這個SNP是第j個組織的一個eQTL，則A_{i j}取1，否則取0。將需要分析的SNP與已有的功能型數(shù)據(jù)庫進行匹配，就可以產(chǎn)生矩陣A。下面以一種常用的統(tǒng)計模型為例來闡述GWAS數(shù)據(jù)與功能型數(shù)據(jù)的集成分析。

假設(shè)觀測到n個個體樣本的表型數(shù)據(jù)（記為y∈R_n）及對應(yīng)的基因型數(shù)據(jù)（用矩陣X^n×Mn×M表示），那么，對于第i個樣本個體，其表型數(shù)據(jù)yi與其基因型數(shù)據(jù)X_i1,?,X_iM可以建立以下線性關(guān)系：

其中，β₁?β_M為回歸系數(shù)，εi 為殘差項，且服從正態(tài)分布N(0,σ2)。若回歸系數(shù)β_j不為0，則代表第i個變異位點與疾病相關(guān)，因此進一步假設(shè)存在一個指示變量γ_j，γ_j可取0和1，分別代表β_j等于0和β_j不等于0。具體來說，β_j與γ_j的關(guān)系可以通過以下模型建立：

接下來通過一個Logistic模型建立γ_j與第i個SNP的功能型數(shù)據(jù)信息A_j（即矩陣A中的第j行）的關(guān)系：

其中，θ∈R^K和θ₀是需要估計的參數(shù)。若θ中的某個元素不為0，代表疾病的致病基因可能與對應(yīng)的功能組相關(guān)。也就是說，如果知道某個SNP的功能信息，引入這些信息能夠調(diào)整γ_j的先驗分布，模型會自動給與疾病相關(guān)性高的某種功能型的位點賦予更高的權(quán)重，從而提高統(tǒng)計分析的效率。若模型的參數(shù)估計結(jié)果已經(jīng)得到，通過計算在給定y、G、A下γ_j取1或0的后驗概率，可以推斷γ_j取1或0的概率，從而可以檢測該變異位點是否屬于風險變異點。對模型的參數(shù)進行統(tǒng)計推斷可以對該變異位點進行功能分析。

上述模型的框架可以擴展到GWAS概括性數(shù)據(jù)建模分析中，如可以對概括性統(tǒng)計量的假設(shè)檢驗的p值進行建模，通過假設(shè)p值來自一個由均勻分布和Beta分布組成的混合分布，建立起p值（如p_i）與指示變量γ_i的聯(lián)系。具體來說，假設(shè)觀測得到的p值用p_i,?,p_n表示，即與疾病不相關(guān)的SNP來自[0,1]區(qū)間的均勻分布U（0,1），對應(yīng)的指示變量γ_j此時取1，而與疾病相關(guān)的SNP的p值來自Beta分布Betta(a,1)，對應(yīng)的指示變量取0，具體如下：

類似地，可以通過Logistic模型（見計算式（6））建立起γ_j與第j個SNP的功能型數(shù)據(jù)信息A_j的聯(lián)系。最大期望算法常用于求解上述模型的參數(shù)，目前的研究已經(jīng)可以滿足集成大規(guī)模功能型數(shù)據(jù)的分析和計算的需求。

近年來，深度學習作為一種機器學習技術(shù)開始得到越來越多的關(guān)注，在圖像識別和自然語言處理等領(lǐng)域取得了突破性進展。因其優(yōu)異的數(shù)據(jù)處理能力，深度學習也開始被逐步應(yīng)用于探索人類基因組密碼，尤其是功能基因（functional genomics）組學。例如，2015年提出的DeepSEA，該方法將公開的多組學功能型數(shù)據(jù)（ENCODE功能型數(shù)據(jù)，表觀基因組的圖譜計劃）作為輸入，搭建深度學習模型，從而實現(xiàn)對非編碼區(qū)域遺傳變異位點（占人體基因組的98%）的功能分析。2016年提出的DeepWAS首先基于DeepSEA找到潛在的影響疾病的某一類型的功能型位點，然后再構(gòu)造帶懲罰函數(shù)的多元回歸模型（LASSO），實現(xiàn)風險變異位點的檢測。單獨的GWAS分析只能確定某個變異點與疾病的關(guān)聯(lián)性，DeepWAS則能夠更加綜合地分析致病基因突變，幫助尋找其致病機理。實際應(yīng)用顯示，集成分析GWAS數(shù)據(jù)和功能型數(shù)據(jù)不僅可以幫助檢測新的危險變異，而且對探索疾病的發(fā)病機制有非常重要的作用。

3.1.2 結(jié)合基因表達數(shù)據(jù)的集成分析

結(jié)合基因表達數(shù)據(jù)的集成分析方法有PrediXcan、MetaXcan、TWAS以及CoMM。該類型方法的基本思想是以已有的基因表達樣本庫的數(shù)據(jù)為參考面，估算GWA S的個體的基因表達，然后檢測基因表達和表征是否關(guān)聯(lián)，其基本原理如圖5所示，此類方法的目的在于對基因調(diào)控的表達部分和性狀進行關(guān)聯(lián)分析（圖5中實線箭頭）。以CoMM為例。假設(shè)參考面的某個基因的基因表達數(shù)據(jù)以及基因型數(shù)據(jù)為y_1g和X_1g，感興趣的疾病或性狀的表型和基因型數(shù)據(jù)為y₂和X _2g，這里基因的表達數(shù)據(jù)和GWAS個體層面的數(shù)據(jù)是不同的樣本，基因表達型數(shù)據(jù)和GWAS數(shù)據(jù)的關(guān)系通過共享參數(shù)u來建立，具體如下：

圖5 基因表達的數(shù)據(jù)和疾病的GWAS數(shù)據(jù)的集成分析基本原理

其中，α_g為未知參數(shù)，通過檢驗α_g是否為0，可以推斷基因表示與表征的關(guān)系?；诙嗷蛐?yīng)的假設(shè)，這里假設(shè)u服從以0為均值，以

為方差的多元正態(tài)分布，即u～N(0,) 。e₁、e₂為誤差項且服從零均值正態(tài)分布。該方法將基因表達數(shù)據(jù)和疾病或性狀的GWAS數(shù)據(jù)進行聯(lián)合建模和估計，考慮了估算基因表達的不確定性，顯示出更高的統(tǒng)計效率。這類型方法也可以推廣到其他層面的組學數(shù)據(jù)（如蛋白組層面）和GWAS的集成分析研究中，具有極高的應(yīng)用價值。

3.2 遺傳變異的多效性

遺傳變異的多效性的存在決定人們可以集成多個GWAS數(shù)據(jù)進行分析。一方面，可以探討復(fù)雜疾病在遺傳結(jié)構(gòu)上的相關(guān)關(guān)系，這種相關(guān)關(guān)系在某種程度上可以反映復(fù)雜疾病之間的聯(lián)系；另一方面，利用這種多效性設(shè)計有效的統(tǒng)計模型，可以幫助人們找到更多與多種復(fù)雜疾病相關(guān)的多效性位點，同時，通過引入功能型數(shù)據(jù)進行分析，能進一步提高人們對疾病的遺傳機理的認識。

復(fù)雜疾病在遺傳結(jié)構(gòu)上的相關(guān)關(guān)系可以通過遺傳相關(guān)系數(shù)這一參數(shù)來度量，該參數(shù)反映兩個復(fù)雜性狀中由遺傳效應(yīng)引起的總體相關(guān)性程度。常見疾病（如心臟病、糖尿病、脂類代謝異常以及高血壓）兩兩之間存在較強的遺傳相關(guān)系數(shù)（0.27~0.43），反映了遺傳變異在這一類型的疾病中具有比較一致的效應(yīng)。基于線性混合模型，可以先了解遺傳相關(guān)系數(shù)的定義。假設(shè)已經(jīng)獲取兩個獨立樣本的關(guān)于兩個表征的觀測數(shù)據(jù)以及個體層面的基因測序數(shù)據(jù)（標準化后），表示為：y₁∈R^n1×1， G₁∈R^n1×m，y₂∈R^n2×1，G₂∈R^n2×1，對應(yīng)的協(xié)變量信息用X₁∈R^n1×p1、X₂∈R^n2×p2表示，這里的下標對應(yīng)兩個不同表征型，它們之間的關(guān)系由以下線性混合模型來建立：

基于多基因效應(yīng)，假設(shè)：

其中，

分別

為u₁和u₂的方差，ρ為u₁和u₂的相關(guān)系數(shù)，即遺傳相關(guān)系數(shù)，其取值范圍為[-1,1]。極大似然估計或限制性極大似然估計等方法常用于二元混合模型的參數(shù)估計。此外，基于概括性數(shù)據(jù)的LDScore回歸方法，可以針對任意兩個性狀進行分析，并且能有效地處理樣本重合的影響。在該方法的基礎(chǔ)上，研究者們提出了GNOVA，該方法可以計算不同性狀在功能型位點的分區(qū)遺傳相關(guān)程度，能夠幫助人們更具體地了解疾病的遺傳特征。更重要的是，集成多個GWAS數(shù)據(jù)的分析方法可以幫助人們進行多效性位點的檢測。研究者們通過對比不同性狀或疾病的GWAS研究的顯著位點，可以標注多效性基因位點。然而這種方法只考慮了顯著位點，因而缺乏統(tǒng)計效率。通過集成分析多個不同性狀的全基因組的GWAS數(shù)據(jù)（不只考慮顯著位點），建立有效的統(tǒng)計模型，可以檢測更多的風險變異點，并發(fā)現(xiàn)更多的多效性位點，進而提高統(tǒng)計分析的效率。多變量的線性混合模型也常被用于關(guān)聯(lián)變異位點的檢測，顯示出比一元線性混合模型更大的統(tǒng)計效率。目前常用的工具有GEMMA、mvLMM和BOLTREML等。該類型方法在實際應(yīng)用中往往受到很多限制，原因在于這類型方法需要用到GWAS個體層面的數(shù)據(jù)。受數(shù)據(jù)可獲得性和數(shù)據(jù)隱私保護等的限制，獲取GWAS個體層面的數(shù)據(jù)往往會比較困難?；诟爬ㄐ詳?shù)據(jù)的多個GWAS的集成分析，只需要GWAS概括性數(shù)據(jù)，而且不同性狀的GWAS數(shù)據(jù)可以來自不同的生物樣本集（考慮到基因的種群效應(yīng)，一般要求GWAS數(shù)據(jù)來自同種群），有利于大規(guī)模的集成分析。同時，該類型方法能夠更靈活地引入功能型數(shù)據(jù)，提高模型的統(tǒng)計效率。基于隱變量的統(tǒng)計模型在GWAS的集成分析中扮演著重要的角色。以基于基因多效性和功能型數(shù)據(jù)的遺傳分析（genetic analysis incorporating pleiotropy and annotation,GPA）方法為例，該方法基于概括性統(tǒng)計量（p值），集成分析多個GWA S數(shù)據(jù)和功能型數(shù)據(jù)，能有效地檢測出多效性位點。假設(shè)收集到K個GWAS的M個SNP的p值，存放在一個M×K維的矩陣中，用P表示，其中p_jk代表第j個SNP在第k個GWAS的p值。類似于單個GWAS危險變異的檢測，這里也假設(shè)存在一個M×K維的隱變量矩陣Z，Z_jk取0代表第j個SNP和第k個GWAS無關(guān)，Z_jk取1代表第j個SNP和第k個GWAS有關(guān)。變量之間對應(yīng)的關(guān)系可以通過下列模型表示：

以k=2為例，定義：

這里π₀₀、π₀₁、π₁₀、π₁₁分別表示Z _j1、Z_j2取不同值的比例，并且π₀₀+π₀₁+π₁₀+π₁₁=1，進一步地，GPA假設(shè)在給定GWAS的關(guān)聯(lián)狀態(tài)下，SNP的功能性狀態(tài)之間獨立。

在GPA模型中，π₀₀,π₀₁,π₁₀,π₁₁,q_00k,q_01k,q_10k,q_11k以及a₁,a₂為模型的未知參數(shù)，EM類型的算法可以有效地解決該模型的參數(shù)估計問題。模型通過計算局部錯誤發(fā)現(xiàn)率（false positive rate， FDR）Pr(Z_j1,Z_j2|P,A)、Pr(Z_j1|P,A)和和P和Pr(Z_j2|P,A)來檢測多效性位點以及單個GWAS的風險變異點；對 q_00k,q_01k,q_10k,q_11k,進行顯著性檢驗可以檢測與疾病相關(guān)聯(lián)的位點在第k類功能型上的富集情況。

進一步的研究顯示，通過Logistic回歸和Probit回歸實現(xiàn)隱變量和功能型數(shù)據(jù)的建模，可以解決GPA模型只能引入離散型、功能型數(shù)據(jù)的局限。同時，在多個GWAS數(shù)據(jù)的集成分析中，低秩性和稀疏性的疊加模型以及復(fù)合極大似然方法的使用，完美地避免了模型參數(shù)指數(shù)增長而無法估計的問題。

3.3 基于孟德爾隨機化法的復(fù)雜性狀因果推斷

基因的多效性普遍存在于人類基因組中，對復(fù)雜疾病的研究具有重要的價值?；虻闹薪槎嘈钥梢杂脕硖剿鲝?fù)雜疾病或性狀之間的因果關(guān)系，這類型的方法被稱作孟德爾隨機化法（Mendelian randomization，MR）。該方法可追溯至1986年，Katan M B提出載脂蛋白E （APOE）基因的變異（包括E2/E3/E4）能夠影響體內(nèi)膽固醇的水平，其中攜帶E2基因的人血漿膽固醇水平較低。由于親代到子代等位基因的隨機分配，個體APOE基因的攜帶情況不受其他混雜因素的影響。因此，可通過癌癥病人和非癌癥病人的E2基因攜帶情況是否存在差異來判斷低水平血漿膽固醇是否會增加癌癥風險。孟德爾隨機化法這一術(shù)語在1991年第一次被提出，逐漸成為因果推斷中的主流方法。該方法以基因變異（如APOE）為工具變量來研究暴露因素（如血漿膽固醇水平）和結(jié)局變量（如癌癥）之間的因果關(guān)系，可以突破傳統(tǒng)的觀察性流行病學在研究因果關(guān)系方面的多種限制，如反向因果、混雜因素等。近幾年來，基于GWAS概括性數(shù)據(jù)的MR分析方法的發(fā)展更能體現(xiàn)出其方法上的優(yōu)越性。首先，大量GWAS概括性數(shù)據(jù)的可公開獲得為研究者們提供了豐富的數(shù)據(jù)資源。MR檢驗可以在任意的性狀之間進行，大大擴展了因果關(guān)系的可研究范圍。相比之下，傳統(tǒng)的隨機控制實驗往往受到實驗倫理以及試驗設(shè)計的局限，對于罕見疾病的研究，往往因數(shù)據(jù)收集方面的困難而無法進行；其次，GWAS的研究基本上是基于大樣本的研究，目前還在不斷朝著更大樣本規(guī)模發(fā)展，MR因果分析也將受益于次。

孟德爾隨機化法以遺傳變異（G）為工具變量，在推斷暴露因素（X）對結(jié)局變量（Y）的因果關(guān)系時，需要滿足以下3條基本條件：①G和X相關(guān)；②假設(shè)存在混雜因素（U）同時影響暴露因素和結(jié)局變量，G和U無關(guān)；③G只能通過X影響Y，而不能存在其他的通道影響Y。只有在這3個條件同時滿足的情況下，MR方法才能正確地推斷X和Y的因果關(guān)系 （如圖6所示）。具體來說，條件①要求在做MR分析時，必須選擇與X顯著相關(guān)的位點，比如p<5×10-8。根據(jù)孟德爾遺傳定律，總是假設(shè)條件②成立，然而其他因素（如群體分層、樣本的重合）可能會導致這個條件不符合，因此在實際應(yīng)用中選擇的GWAS數(shù)據(jù)需來自同一種群且應(yīng)避免存在樣本重合；條件③要求G對Y不能有直接的影響，G對X和Y只能存在中介多效性?；蛑薪槎嘈砸渤１环Q作垂直多效性，對應(yīng)的基因的生物多效性也稱作水平多效性，用以描述基因與性狀之間的直接聯(lián)系。水平多效性的存在會影響MR分析結(jié)果的可靠性。

圖6 孟德爾隨機化法基本原理

MR分析中各變量間（包括G、X、Y、U）的關(guān)系可通過下面的線性結(jié)構(gòu)方程來描述：

其中，

為M個獨立的變異位點，γ_j為G _j對X的效應(yīng)大?。╡ffect size），α_j為G_j對Y的效應(yīng)大小。若該變異位點滿足MR分析的假設(shè)，則α_j=0；若存在水平多效性，則α_j≠0。η_x和η_j表示混雜因素對X和Y的影響大小，E_x和E_y為對應(yīng)的殘差向量，β表示X對Y的影響效應(yīng)的大小，若β不為0，則表示X和Y存在因果關(guān)系。MR分析的目的是準確地估計β的大小并進行因果關(guān)系的檢驗。上述模型要求GWAS個體層面的數(shù)據(jù)，而基于兩樣本的GWAS概括性統(tǒng)計量的MR分析方法不受數(shù)據(jù)的隱私保護等限制，逐漸成為熱點。這時，觀測到的數(shù)據(jù)為變異位點對暴露因素（X）和結(jié)局變量（Y）的回歸效應(yīng)（effect）的估計值及其標準差，記為

，它們的關(guān)系可以表示為：

由于兩組數(shù)據(jù)來自無樣本重合的GWAS研究，γ?j與 Γ?j相互獨立（給定的情況下）。若變異位點均滿足MR假設(shè)條件，不存在水平多效性，則α_j=0，j=1,?n；若存在某幾個α_j不為0，由于會影響β的估計，則會被當作異常點。若所有位點均存在水平多效性但平均水平為0，則稱為平衡的水平多效性或系統(tǒng)性多效性，反之則稱為非平衡的水平多效性或有向多效性。

目前關(guān)于MR的大量研究集中在處理變異位點的水平多效性方面，研究方法可大致分為3種類型。一類是矯正，如MREgger通過在傳統(tǒng)的異方差加權(quán)模型（inverse variance weighted，IVW）中引入常數(shù)項，修正由于非平衡的水平多效性（αj的均值不為0）對模型帶來的影響。該方法往往會引入較大的估計方差，降低統(tǒng)計效率。第二類是異常值剔除，這類型方法（如MR-PRESSO、GSMR等）先對可能存在水平效應(yīng)的異常位點進行檢測，剔除檢測到的異常位點，然后對模型參數(shù)重新估計。這類方法在存在少量的異常點時是有效的，但當存在大量的異常點或平衡的水平效應(yīng)時，則需謹慎使用。第三類方法是對違背MR假設(shè)的變異點進行具體的建模，如MR.raps和BWMR等。MR.raps和BWMR均假設(shè)存在系統(tǒng)性的水平效應(yīng)，而且αj服從零均值的正態(tài)分布，MR.raps在使用似然方法進行參數(shù)估計時，用穩(wěn)健的損失函數(shù)Tukey損失函數(shù)和Huber損失函數(shù)取締了均方損失函數(shù)，以得到穩(wěn)健的估計。BWMR則引入貝葉斯重加權(quán)（Bayesian reweighting）的思想，通過賦予異常點較低的權(quán)重來降低異常點的影響。

基于概括型統(tǒng)計量的MR分析，集成來自不同樣本的GWAS研究數(shù)據(jù)，可以推斷出不同性狀之間的因果關(guān)系，對了解復(fù)雜疾病的發(fā)病機制、疾病的干預(yù)治療以及制藥等方面具有重要的推進作用。

4 應(yīng)用案例

4.1 基因多效性分析：運用于神經(jīng)質(zhì)和重度抑郁癥GWAS數(shù)據(jù)分析

心理學上描述人的五大性格特質(zhì)包括：開放性、責任性、外傾性、宜人性和神經(jīng)質(zhì)，其中神經(jīng)質(zhì)主要反映平衡焦慮、敵對、壓抑、自我意識、沖動、脆弱等情緒的能力，即維持情緒穩(wěn)定性的能力。已有的研究表明，神經(jīng)質(zhì)和重度抑郁癥（major depression disorder，MDD）存在較強的關(guān)聯(lián)，在遺傳上存在較強的相關(guān)性。本文基于GPA方法對已有的神經(jīng)質(zhì)和重度抑郁癥的GWAS概括型數(shù)據(jù)和功能型數(shù)據(jù)（這里使用eQTL數(shù)據(jù)）進行基因的多效性分析，以探索這兩種性狀在遺傳上的關(guān)系。集成分析結(jié)果見表1，曼哈頓圖如圖7所示。集成兩組GWAS數(shù)據(jù)進行分析顯示出更高的統(tǒng)計效率，原因在于這兩組疾病具有高度相似的遺傳學機理。根據(jù)參數(shù)估計結(jié)果，π?00=0.793，π?01=0.014，π?10=0，π?11=0.194，其中，π?11=0.194（似然比檢驗p=0）表明神經(jīng)質(zhì)和重度抑郁癥具有高度相似的遺傳機理，也就是說兩個性狀共享信息，共同分析兩個GWAS數(shù)據(jù)能極大地提高效率。此外，表1的結(jié)果顯示通過集成功能型數(shù)據(jù)，能發(fā)現(xiàn)更多的顯著位點，進一步提高統(tǒng)計效率。

圖7 神經(jīng)質(zhì)和重度抑郁癥的 曼哈頓圖（藍線和紅線分別對應(yīng)局部錯誤率為0.1和0.05）

4.2 血脂代謝物與人類復(fù)雜疾病的因果關(guān)系

本節(jié)基于101種分子水平的血脂代謝物和57種人類復(fù)雜疾病的GWAS概括型數(shù)據(jù)，探討集成分析在實際中的應(yīng)用。復(fù)雜性狀或疾病的GWAS數(shù)據(jù)來源見表2，血脂代謝物包含87種血脂蛋白攜帶的脂類物質(zhì)以及14種游離的脂肪酸或相關(guān)度量。脂類物質(zhì)包括總膽固醇（total cholesterol，C）、膽固醇酯（cholesterol esters，CE）、總脂質(zhì)（total lipids， L）、磷脂（phospholipids，P）、游離膽固醇（free cholesterol，F(xiàn)C）、甘油三酯（triglycerides，TG）。血脂蛋白根據(jù)密度可分為低密度脂蛋白（low-density lipoprotein，LDL）、高密度脂蛋白（highdensity lipoprotein，HDL）、極低密度脂蛋白（very-low lipoprotein，VLDL）、中密度脂蛋白（intermediate-density lipoprotein，IDL）。各類血脂蛋白根據(jù)大小可以進行細分，如LDL可分為L.LDL、M.LDL、S.LDL等。57種人類復(fù)雜疾病包括與人體測量相關(guān)的性狀（如身高、體脂指數(shù)（BMI）等）、與心血管相關(guān)性狀（如冠心病（CAD）、收縮壓（SBP）、舒張壓（DBP）等）、代謝類疾?。ㄈ纰蛐吞悄虿?、脂類代謝異常等）、中樞神經(jīng)系統(tǒng)疾?。ㄈ绨柎暮ＤY、帕金森綜合征等）、精神疾?。ㄈ缰囟纫钟舭Y、多動癥等）以及自身免疫疾?。ㄈ纰裥吞悄虿　⑾到y(tǒng)性紅斑狼瘡等）。具體的分析分兩方面：一方面，基于基因的多效性系統(tǒng)地探討血脂代謝物和復(fù)雜疾病在遺傳上的相關(guān)關(guān)系，對了解復(fù)雜疾病的遺傳結(jié)構(gòu)和性質(zhì)具有重要的作用，本文基于GNOVA方法，對血脂代謝物和與人類復(fù)雜疾病的遺傳相關(guān)系數(shù)進行估計；另一方面，血脂代謝物和復(fù)雜疾病的因果關(guān)系一直是科學家們關(guān)注的重點，對于疾病的預(yù)防和診斷具有重要的指導作用。以心血管疾病為例，目前大量的研究表明，血液中的低密度膽固醇是造成心血管疾病的主要原因。本文應(yīng)用BWMR方法，對血脂代謝物和人類復(fù)雜疾病的因果關(guān)系和反向因果關(guān)系進行分析，從而深化對血脂代謝物和復(fù)雜疾病的認識。

圖8展示的是通過Bonferroni校正之后顯著的遺傳相關(guān)系數(shù)（p<0.05/(101×57)），其中640對具備顯著的遺傳相關(guān)關(guān)系。圖8中紅色為正相關(guān)，綠色為負相關(guān)，顏色深度和方格面積表示相關(guān)系數(shù)的大小，顏色越深或方格面積越大表示遺傳相關(guān)系數(shù)越大。遺傳相關(guān)系數(shù)的估計需調(diào)用GNOVA方法的軟件包。結(jié)果顯示，大部分復(fù)雜疾病與大部分血脂代謝物水平均有顯著的遺傳上的相關(guān)關(guān)系，如心血管疾病、高血壓疾病、血壓相關(guān)度量、體重相關(guān)度量等。此外，阿爾茨海默癥、神經(jīng)性厭食癥與多種類型的血脂代謝物有密切的關(guān)系。值得注意的是，雙相情感障礙癥與多種游離脂肪酸有顯著的遺傳上的相關(guān)關(guān)系。因果分析結(jié)果顯示，體質(zhì)指數(shù)、Ⅱ型糖尿病會影響體內(nèi)血脂代謝物的水平（結(jié)果如圖9所示）。

圖8 血脂代謝物（101種）與復(fù)雜疾病或性狀（57種）的相關(guān)系數(shù)

圖9 血脂代謝物和復(fù)雜性狀的因果分析結(jié)果

5 結(jié)束語

本文對GWA S研究中的集成分析進行了綜述，主要應(yīng)用于以下3個方面：檢測風險位點及其功能分析、基因多效性的分析、基于孟德爾隨機化的因果推斷。實際應(yīng)用顯示，集成分析在GWAS研究中顯示出重要的作用，有助于挖掘重要的信息。同時，多種類型的集成分析統(tǒng)計方法具備極大的拓展空間，將對未來的研究發(fā)揮更大的價值。未來GWAS的集成分析將更多地應(yīng)用于集成多組學數(shù)據(jù)的分析，對探索從遺傳變異到疾病發(fā)生的整個因果鏈條起著重要作用。隨著越來越多的高質(zhì)量數(shù)據(jù)的不斷產(chǎn)生，全方位生物醫(yī)療大數(shù)據(jù)（包括基因大數(shù)據(jù)、醫(yī)療圖像數(shù)據(jù)、電子病歷等）的集成分析將使得實現(xiàn)個性化的精準醫(yī)療成為可能。