B2016年優(yōu)秀論文BK0209

1、- 1 - 參賽密碼 （由組委會填寫） 全全全全 “華為華為華為華為杯杯杯杯”第十三屆全國研究生第十三屆全國研究生第十三屆全國研究生第十三屆全國研究生 數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽 學(xué)學(xué)校校華東師范大學(xué)華東師范大學(xué) 參賽隊號參賽隊號K0209 隊員姓名隊員姓名 1.賈柯賈柯 2.崔軒崔軒 3.陳陳嘉嘉駿駿 - 2 - 參賽密碼 （由組委會填寫） “華為杯華為杯華為杯華為杯”第十三屆全國研究生第十三屆全國研究生第十三屆全國研究生第十三屆全國研究生 數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽數(shù)學(xué)建模競賽 題題 目目具有遺傳性疾病和性狀的遺傳位點分析具有遺傳性疾病和性狀的遺傳位點分析

2、摘要： 大量研究表明，人體的許多表型性狀差異以及對藥物和疾病的易感性等都 可能與某些位點相關(guān)聯(lián)，或和包含有多個位點的基因相關(guān)聯(lián)。因此，定位與性 狀或疾病相關(guān)聯(lián)的位點在染色體或基因中的位置，能幫助研究人員了解性狀和 一些疾病的遺傳機(jī)理，也能使人們對致病位點加以干預(yù)，防止一些遺傳病的發(fā) 生。 對于問題一對于問題一，根據(jù)位點中堿基對的特征，基于生物基因的加性效應(yīng)，位點 中的三種組合分別可編碼為 0、1,、2，其中 1 代表雜合子基因，0 和 2 分別代 表純合子基因中的主要等位基因(major allele)與次要等位基因(minor allele)。 對于問題二對于問題二，需通過一定方法計算出每個

3、位點與疾病之間的關(guān)聯(lián)程度，本 文首先通過卡方檢驗進(jìn)行建模，并且分別用 Benjamini BONF0.0006024;BONF 校正后校正后 P P 值為值為 0.00060240.0006024） 。 除此之外， 還采用置換檢驗?zāi)Ｐ秃?貝葉斯模型進(jìn)行檢驗，檢測出來最顯著的致病位點與卡方檢驗相一致，因此最 終得出與疾病相關(guān)的位點有一個，位點名稱為：rs2273298rs2273298（置換檢驗?zāi)Ｐ托Ｕ脫Q檢驗?zāi)Ｐ托Ｕ?后后 p p 值為值為 0.0094450.009445，貝葉斯因子，貝葉斯因子的對數(shù)取值的對數(shù)取值為為 4.512384.51238） 。 對于問題三對于問題三，根據(jù)基因可以表

4、示為位點的集合這一特征，本文采用 Set-based test 和 VEGAS 模型對一個基因內(nèi)連鎖不平衡的 SNP 位點進(jìn)行建模， 并且都采用置換算法進(jìn)行模型求解，閾值為經(jīng)過 BH 校正后 P 值小于 0.05，最 - 3 - 終兩個模型得到了一致且較好的效果。與疾病相關(guān)聯(lián)的基因有三個，基因所屬 序列為：gene_gene_5555、gene_gene_102102、gene_gene_217217（Set-basedSet-based testtest 模型模型 BHBH 校正校正 P P 值后值后 分 別 為分 別 為 0.149985,0.04,0.1499850.149985,0.0

5、4,0.149985 ； VEGASVEGAS 模 型模 型 BHBH 校 正 后校 正 后 P P 值 分 別 為值 分 別 為 0.00165,0.0184,0.00090.00165,0.0184,0.0009） 。 對于問題四對于問題四，多個性狀往往表現(xiàn)為一個整體來進(jìn)行衡量，本文分別采用 mv-plink 模型和 MultiPhen 模型對多個表型之間的關(guān)聯(lián)進(jìn)行建模，并且找出這 些關(guān)聯(lián)表型的致病位點。閾值選取為經(jīng)過 BH 校正后 P 值小于 0.05。最終兩個 模型都得出與樣本中十個性狀有關(guān)聯(lián)的位點有一個，位點名稱為：rs12746773rs12746773 （mv-plinkmv-p

6、link 模型模型 BHBH 校正后校正后 P P 值為值為 2.8684472.868447 20 10； MultiPhenMultiPhen 模型模型 BHBH 校正校正后后 P P 值為值為 8.3061988.306198 21 10） 。 關(guān)鍵詞：關(guān)鍵詞：遺傳統(tǒng)計學(xué)， 全基因組關(guān)聯(lián)性分析(GWAS)，位點(SNPs) ，卡方 檢驗 - 4 - 目 錄 一、 問題描述.- 5 - 二、二、 合理假設(shè)與符號說明. - 7 - 2.1 合理假設(shè).- 7 - 2.2 符號說明.- 7 - 三、 問題分析.- 8 - 3.1 問題一.- 8 - 3.2 問題二.- 8 - 3.3 問題三.-

7、 8 - 3.4 問題四.- 8 - 四、 模型特點介紹.- 9 - 4.1 問題二的建模.- 9 - 4.1.1 卡方檢驗?zāi)Ｐ?- 9 - 4.1.1.2 列聯(lián)表的獨立性檢驗?zāi)Ｐ? - 9 - 4.1.2 基于貝葉斯的 GWAS 模型.- 12 - 4.1.3 置換檢驗?zāi)Ｐ? - 14 - 4.2 問題三的建模.- 15 - 4.2.1 基于集合的基因檢驗?zāi)Ｐ停⊿et-based test）.- 15 - 4.2.2 全面基于基因關(guān)聯(lián)分析模型(VEGAS).- 16 - 4.3 問題四的建模.- 17 - 4.3.1 基于典型關(guān)聯(lián)分析的多表型模型（MV-Plink）.- 17 - 4.3.2

8、 MultiPhen 模型. - 18 - 五、 問題求解.- 19 - 5.1 問題一求解.- 19 - 5.2 問題二求解.- 19 - 5.3 問題三求解.- 22 - 5.4 問題四求解.- 23 - 六、 模型評價.- 25 - 參考文獻(xiàn).- 26 - 附件.- 27 - - 5 - 一、一、問題描述問題描述 人體的每條染色體攜帶一個 DNA 分子， 人的遺傳密碼由人體中的 DNA 攜帶。 DNA 是由分別帶有 A,T,C,G 四種堿基的脫氧核苷酸鏈接組成的雙螺旋長鏈分子。 在這條雙螺旋的長鏈中，共有約 30 億個堿基對，而基因則是 DNA 長鏈中有遺傳 效應(yīng)的一些片段。在組成 DN

9、A 的數(shù)量浩瀚的堿基對（或?qū)?yīng)的脫氧核苷酸）中， 有一些特定位置的單個核苷酸經(jīng)常發(fā)生變異引起 DNA 的多態(tài)性，我們稱之為位 點。染色體、基因和位點的結(jié)構(gòu)關(guān)系見圖 1-1。 在 DNA 長鏈中，位點個數(shù)約為 堿基對個數(shù)的 1/1000。由于位點在 DNA 長鏈中出現(xiàn)頻繁，多態(tài)性豐富，近年來 成為人們研究 DNA 遺傳信息的重要載體，被稱為人類研究遺傳學(xué)的第三類遺傳 標(biāo)記。 圖 1-1 染色體、基因和位點的結(jié)構(gòu)關(guān)系 大量研究表明，人體的許多表型性狀差異以及對藥物和疾病的易感性等都 可能與某些位點相關(guān)聯(lián)，或和包含有多個位點的基因相關(guān)聯(lián)。因此，定位與性 狀或疾病相關(guān)聯(lián)的位點在染色體或基因中的位置，能

10、幫助研究人員了解性狀和 一些疾病的遺傳機(jī)理，也能使人們對致病位點加以干預(yù)，防止一些遺傳病的發(fā) 生。近年來，研究人員大都采用全基因組的方法來確定致病位點或致病基因， 具體做法是：招募大量志愿者（樣本） ，包括具有某種遺傳病的人和健康的人， 通常用 1 表示病人，0 表示健康者。 對每個樣本，采用堿基(A,T,C,G)的編碼 方式來獲取每個位點的信息(因為染色體具有雙螺旋結(jié)構(gòu)， 所以用兩個堿基的組 合表示一個位點的信息） 。 本文提出的問題包括： 一、請用適當(dāng)?shù)姆椒ǎ?把 genotype.dat 中每個位點的堿基(A,T,C,G) 編碼 方式轉(zhuǎn)化成數(shù)值編碼方式，便于進(jìn)行數(shù)據(jù)分析。 二、根據(jù)附錄中

11、 1000 個樣本在某條有可能致病的染色體片段上的 9445 個 位點的編碼信息(見 genotype.dat)和樣本患有遺傳疾病 A 的信息（見 phenotype.txt 文件） 。設(shè)計或采用一個方法，找出某種疾病最有可能的一個或 幾個致病位點，并給出相關(guān)的理論依據(jù)。 三、同上題中的樣本患有遺傳疾病 A 的信息（phenotype.txt 文件） ?，F(xiàn)有 300 個基因，每個基因所包含的位點名稱見文件夾 gene_info 中的 300 個 dat 文件，每個 dat 文件列出了對應(yīng)基因所包含的位點(位點信息見文件 genotype.dat)。由于可以把基因理解為若干個位點組成的集合，遺傳

12、疾病與基 - 6 - 因的關(guān)聯(lián)性可以由基因中包含的位點的全集或其子集合表現(xiàn)出來請找出與疾病 最有可能相關(guān)的一個或幾個基因，并說明理由。 四、在問題二中，已知 9445 個位點，其編碼信息見 genotype.dat 文件。 在實際的研究中，科研人員往往把相關(guān)的性狀或疾病看成一個整體，然后來探 尋與它們相關(guān)的位點或基因。試根據(jù) multi_phenos.txt 文件給出的 1000 個樣 本的 10 個相關(guān)聯(lián)性狀的信息及其 9445 個位點的編碼信息(見 genotype.dat)， 找出與 multi_phenos.txt 中 10 個性狀有關(guān)聯(lián)的位點。 - 7 - 二、二、合理假設(shè)與符號說明

13、合理假設(shè)與符號說明 2.12.1 合理假設(shè)合理假設(shè) 根據(jù)題意，可以進(jìn)行如下假設(shè)： 1) 假設(shè)附件的樣本均來自同一人種； 2) 性別對該疾病影響不大，即不考慮性別對疾病的影響； 3) 假設(shè)題目中所給數(shù)據(jù)都是有效數(shù)據(jù)； 2.22.2 符號說明符號說明 符號定義 GWAS全基因組關(guān)聯(lián)分析 SNP單核苷酸的多態(tài)性 0 H 原假設(shè) 1 H 備擇假設(shè) 2 卡方統(tǒng)計量 COV協(xié)方差矩陣 T置換檢驗次數(shù) - 8 - 三、三、問題分析問題分析 3.13.1 問題一問題一 問題一要求附件 genotype.dat 中的位點的堿基轉(zhuǎn)化成數(shù)值編碼。 通過查閱 相關(guān)資料可得，生物基因具有加性效應(yīng)，并且同一個位點均表現(xiàn)為

14、三種不同的 堿基對，因此采用 0,1,2 的編碼方式對位點堿基進(jìn)行數(shù)值編碼。 3.23.2 問題二問題二 問題二要求從附錄中 1000 個樣本在某條有可能致病的染色體片段上的 9445 個位點的編碼信息和樣本患有遺傳疾病 A 的信息中， 找出相關(guān)的致病位點。 針對這個問題，我們可能要采用檢驗獨立性的模型去檢驗疾病與致病位點的關(guān) 系。由于涉及到的假設(shè)檢驗次數(shù)很多，這就會形成多重假設(shè)檢驗問題。這就需 要我們尋找 P 值校正的方法，從而找出可信度較高的致病位點。 3.33.3 問題三問題三 問題三要求找出疾病 A 可能的致病基因。在進(jìn)行單個 SNP 與疾病關(guān)聯(lián)分析 的時候，有許多的 SNP 與疾病的

15、關(guān)系可能處于接近關(guān)聯(lián)的程度，而這些 SNP 往 往很容易被傳統(tǒng)關(guān)聯(lián)分析方法當(dāng)作隨機(jī)噪聲而被舍棄掉。針對這樣的情景，要 根據(jù)一個基因內(nèi) SNP 間連鎖不平衡的關(guān)系，把這些信號相對較低的 SNP 結(jié)合到 基因的層級上，從而增強(qiáng)信號，檢測出與疾病相關(guān)的基因。 3.43.4 問題四問題四 問題四要求找出 10 個相關(guān)聯(lián)性狀的可能致病位點。在實際的研究過程中， 往往把相關(guān)的性狀或疾病看成一個整體。由于要關(guān)注性轉(zhuǎn)間的關(guān)聯(lián)，因此問題 四的解決可能會用到回歸類的模型以及典型顯著分析模型，因為這兩種模型都 能較好地對變量間的關(guān)系進(jìn)行建模。 - 9 - 四、四、模型特點介紹模型特點介紹 4.14.1 問題二的建模

16、問題二的建模 4 4.1.1.1.1 卡方檢驗?zāi)Ｐ涂ǚ綑z驗?zāi)Ｐ?4.1.1.24.1.1.2 列聯(lián)表的獨立性檢驗?zāi)Ｐ土新?lián)表的獨立性檢驗?zāi)Ｐ?我們把按兩個或多個特征分類的頻數(shù)數(shù)據(jù)稱為交叉分類數(shù)據(jù)，這些數(shù)據(jù)一 般以表格的形式給出，稱為列聯(lián)表。一般，若總體中的個體可按兩個屬性 A 與 B 分類，A 有 個類， 1 A， r A，B 有個類 1 B， c B，從總體中抽取大小 為 n 的樣本，設(shè)其中有 ij n個個體既屬于類 i A又屬于 j B， ij n稱為頻數(shù)，將cr個 排列為一個 行 列的二維列聯(lián)表，簡稱為cr列聯(lián)表。若所考慮的屬性多于 兩個，也可以按類似的方式作出列聯(lián)表，稱為多維列聯(lián)表。 列

17、聯(lián)表分析的基本問題是，考慮各屬性之間有無關(guān)聯(lián)，即判別兩屬性是否 獨立。在cr列聯(lián)表中，若以 . i p, j p.和 ij p分別表示總體中的個體僅屬于 i A， 僅屬于 j B和同時屬于 i A和 j B的概率，可得到一個二維離散分布表，則“A、B 兩屬性獨立”的假設(shè)可以表述為公式（4.1） ： ., 1, 1,: .0 cjripppH jiij （4.1） 這里諸 ij p共有rc個參數(shù)，在原假設(shè) 0 H成立時，這rc個參數(shù) ij p由cr 個參 數(shù) . 1 p, . r p和 1 . p, c p.決定，在這后cr 個參數(shù)中存在兩個約束條 件: r i i p 1 . 1, c j j

18、 p 1 . 所以， 此時 ij p實際上由2cr個獨立參數(shù)所確定的。 據(jù)此， 檢驗統(tǒng)計量為公式（4.2） r i c j ij ijij pn pnn 11 2 (4.2) 在原假設(shè) 0 H成立時上式近似服從自由度為12 crrc=11cr 的分布，其中諸 ij p 是在 0 H成立下得到的 ij p的最大似然估計，其表達(dá)式(4.3) - 10 - 為 n n n n ppp j i jiij . . . (4.3) 對給定的顯著性水平(01)，檢驗的拒絕域為公式（4.4） ： 111 22 crW（4.4） 在針對疾病相關(guān)位點檢測的時候，我們這樣使用卡方檢驗：假設(shè)其中一個 SNP 位點基因

19、型的分布情況為 GG,GT,TT。它們在患病樣本和非患病樣本的分布 如表 4-1： 表 4-1 患病和非患病樣本分布表 GGGTTT總數(shù) 患病R 非患病S 總數(shù)N 在進(jìn)行卡方檢驗的時候，需要把表轉(zhuǎn)化為表 4-2： 表 4-2 卡方檢驗轉(zhuǎn)化表 GT總數(shù) 患病2R 非患病2S 總數(shù)2N 4.1.1.24.1.1.2 BenjaminiBenjamini ( d ds m sd d d dsmf （4.21） 其中m為樣本均值， 為樣本方差，d為自由度，如果保持那么可 得公式（4.22）和（4.23）： , )( m)-1)(d - )|(log 22 T l mdsd MPd （4.22） - 1

20、4 - T l mds mdsd d MPd 22 2 2 2 )( )()(1( 0 01 )|(log （4.23） 4.1.34.1.3 置換檢驗?zāi)Ｐ椭脫Q檢驗?zāi)Ｐ?模擬運(yùn)算法的基本原理是：根據(jù)所研究的問題構(gòu)造一個檢驗統(tǒng)計量 , 并利用樣本 ,按排列組合的原理 ,導(dǎo)出檢驗統(tǒng)計量的理論抽樣分布 ;若難以導(dǎo) 出確切的理論分布 ,則采用抽樣模擬的方法做估計其近似分布。然后求出從該 分布中獲得樣本及更極端樣本的概率 ( P 值) ,并界定此概率值作出推論。若 檢驗統(tǒng)計量的抽樣分布是基于樣本的所有可能的排列(或組合) 條件下的分布 , 則稱之 為“Exact Permutation Test ( E

21、PT) ”,可譯為“確切排列 (組合) 檢 驗”,其思路類似于秩和檢驗。 對實際問題來說 ,往往得不到檢驗統(tǒng)計量的確 切抽樣分布 ,可通過基于樣本的大量重復(fù)的隨機(jī) 排列(或組合) 估計其近似的 抽樣分布 ,則稱之為“Randomized Permutation test (RPT) ”,可譯為 “隨 機(jī)排列(或組合) 檢驗”。 在全基因組關(guān)聯(lián)分析中，采用標(biāo)簽置換和基因拋除是采用模擬運(yùn)算法的兩 個基本做法。一般采用自適應(yīng)的模擬運(yùn)算法來估計理論抽樣分布。在自適應(yīng)算 法中，對于沒有太大重要性的位點會比看起來有一些重要性的位點更快地被拋 棄。也就是說，如果在十次對于某一個位點的模擬運(yùn)算過程中有 9 次

22、統(tǒng)計值都 大于觀測到的抽樣統(tǒng)計，那么再繼續(xù)對這個位點模擬運(yùn)算也沒有什么必要了， 因為這個位點已經(jīng)很不可能會顯現(xiàn)出比別的位點更加顯著的結(jié)果。這個步驟可 以大大加快模擬運(yùn)算的過程，因為在剛開始的階段大量沒有意義的位點都會被 拋棄掉，這樣可以使得模擬運(yùn)算到有意義的位點上的概率大大增加。因此很自 然地，那些對疾病有意義的位點統(tǒng)計得到的 P 值會比那些對于疾病沒有太大意 義的位點得到的 P 值更加準(zhǔn)確，但是這并不對我們的結(jié)果造成影響。 自適應(yīng)標(biāo)簽置換檢驗?zāi)Ｐ椭幸还灿?6 個參數(shù)，參數(shù)含義如表 4-3： 表 4-3 參數(shù)含義 參數(shù)表示默認(rèn)值參數(shù)涵義 5每一個位點的最小 置換次數(shù) 1000000每一個位點的

23、最大 置換次數(shù) 0級別閾值 0.0001經(jīng)驗 P 值的置信區(qū)間 1位點調(diào)整區(qū)間 0.001位點調(diào)整區(qū)間減慢 速度 - 15 - 上面參數(shù)表示：對于每一個位點，至少要經(jīng)過次置換才能進(jìn)行位點的調(diào) 整，并且最多不能超過次置換。表示經(jīng)過次之后進(jìn)行位點調(diào)整，并且 下一次調(diào)整間隔次數(shù)增加為R001. 0次，R 是當(dāng)前進(jìn)行過的置換次數(shù)。在 每一次的置換過程中，對每一個期望 P 值的置信區(qū)間為公式（4.24）： % 2 1 100 T CI （4.24） 其中 T 為位點個數(shù)。 相應(yīng)流程圖如圖所示： 圖 4-1 置換檢驗?zāi)Ｐ土鞒虉D 4.24.2 問題三的建模問題三的建模 4.2.14.2.1 基于集合的基因檢

24、驗?zāi)Ｐ停ɑ诩系幕驒z驗?zāi)Ｐ停⊿et-basedSet-based testtest） 兩個 SNP 之間的連鎖不平衡是用來度量兩個 SNP 位點之間的相關(guān)程度。其 中 D 和r平方量化兩個 SNP 位點之間連鎖不平衡程度的指標(biāo)。在這里，基因座 L1 位點為 A 的概率為 A p,基因座 L1 位點為 a 的概率為 a q;基因座 L2 位點為 B 的概率為 B p，基因座 L2 位點為 b 的概率為 b q。因此我們有下表 4-4： - 16 - 表 4-4 L2 BL2 b L1 A DppP BAAB DqpP bAAb L1 a DPqP BaaB DqqP baAB 由上表我們可以

25、得到公式（4.25）和（4.26） ： aBAbabAB PPPPD（4.25） bBaA qpqp D r 2 2 （4.26） 該模型的過程如下： 步驟一：對于每一個基因，我們檢測每個 SNP 與其他 SNP 的是否連鎖不平 衡，檢驗得到的 r 若超過閾值 R 的時候， 則認(rèn)為兩個 SNP 之間是連鎖不平衡的。 步驟二：對于每一個 SNP，我們進(jìn)行單個 SNP 和疾病之間的關(guān)聯(lián)分析，在 這里我們選用卡方檢驗來進(jìn)行單個 SNP 與疾病之間的關(guān)聯(lián)分析。 步驟三：對于每個基因，我們選取 N 個獨立的 SNP（來源于步驟一的連鎖 不平衡檢驗） ，P 值小于步驟一中的閾值 P。首先選取 P 值最小的

26、那個 SNP，去 除與選定 SNP 連鎖不平衡的其他的 SNP，根據(jù)統(tǒng)計量的顯著程度進(jìn)行排列。 步驟四：對于這些 SNP 的子集，每個基因的統(tǒng)計量由基因中的 SNP 的統(tǒng)計 量的均值表示。 步驟五：對數(shù)據(jù)集進(jìn)行大量的置換，保證 SNP 之間的連鎖不平衡狀態(tài)不變。 步驟六：對于每個被置換的數(shù)據(jù)集，重復(fù)步驟二和步驟四。 步驟七：基因的經(jīng)驗 P 值計算，由置換后基因的統(tǒng)計量超過未經(jīng)置換的基 因的統(tǒng)計量的次數(shù)除以所有的置換次數(shù)。 4.2.24.2.2 全面基于基因關(guān)聯(lián)分析模型全面基于基因關(guān)聯(lián)分析模型(VEGAS)(VEGAS) 基于基因的關(guān)聯(lián)分析模型在進(jìn)行關(guān)聯(lián)分析的時候，通?？紤]的是一個基因 內(nèi)的所有

27、 SNP，而不是像傳統(tǒng)的關(guān)聯(lián)分析那樣只考慮一個 SNP。根據(jù)基因組得而 組織結(jié)構(gòu)來看，基于基因的全基因組關(guān)聯(lián)分析會比傳統(tǒng)基于單個 SNP 的關(guān)聯(lián)分 析效果更加好。例如，當(dāng)一個基因包含多個致病的 SNP，而這幾個在相同基因 里的 SNP，在進(jìn)行單個 SNP 的關(guān)聯(lián)分析的時候，可能只是處在接近顯著的程度， 而 SNP 和可能會被當(dāng)成隨機(jī)的噪音被舍棄掉。 當(dāng)把這些同一個基因中的致病 SNP 結(jié)合到一起進(jìn)行相關(guān)的統(tǒng)計檢驗和連鎖不平衡方面的校正，基于基因的方法可 能就會檢測出這些接近顯著的 SNP 相關(guān)的基因。 下面是 VEGAS 模型的過程： 步驟一：對于一個基因中的 n 個 SNP，產(chǎn)生一個 n 維

28、向量 d，該向量服從正 態(tài)分布，該正態(tài)分布的均值為 0，方差矩陣為。其中 為nn矩陣，該矩陣 的元素為成對連鎖不平衡中的 r 值 - 17 - 步驟二： 對矩陣進(jìn)行 Cholesky 矩陣分解， 我們得到nn下三角矩陣 C， 其中 CC。 步驟三：然后用向量 d 與下三角矩陣進(jìn)行相乘，得到 n 維向量 original Q 步驟四：產(chǎn)生一個新的 n 維向量 n zzZ, 1 ,Z 服從正態(tài)分布 , 0 n N, 然后把向量轉(zhuǎn)化為自由量為 1 的相關(guān)卡方變量， 2 1 , iinnew zqqqQ。 步驟五：統(tǒng)計量為 n 維向量 Q 中元素的和。 步驟六：對數(shù)據(jù)進(jìn)行置換，重復(fù)步驟四和步驟五 步驟

29、七：基因的經(jīng)驗 P 值計算，由置換后基因的統(tǒng)計量超過未經(jīng)置換的基 因的統(tǒng)計量的次數(shù)除以所有的置換次數(shù)。 4.34.3 問題四的建模問題四的建模 4.3.14.3.1 基于典型關(guān)聯(lián)分析的多表型模型（基于典型關(guān)聯(lián)分析的多表型模型（MV-PlinkMV-Plink） 典型關(guān)聯(lián)分析是利用綜合變量對之間的相關(guān)關(guān)系來反應(yīng)兩組指標(biāo)之間 的整體相關(guān)性的多元統(tǒng)計分析方法。它的基本原理是：為了把握兩組指標(biāo)之間 的相關(guān)關(guān)系，分別在兩組變量中提取有代表性的兩個綜合變量，利用兩個綜合 變量之間的相關(guān)關(guān)系來反映兩組指標(biāo)之間的整體相關(guān)性。 典型關(guān)聯(lián)分析在這里的應(yīng)用主要是為了尋找表型 Y 和基因型 X 之間的相關(guān) 關(guān)系最大化

30、。這里，XYcorrp, 成為典型相關(guān)系數(shù)估計值。為了得到 ，矩 陣 Y 和矩陣 X 的協(xié)方差矩陣會被分塊為如公式（4.27） ： XXXY YXYY Y X COV（4.27） 其中YY為矩陣 Y 的KK 方差矩陣，XX為矩陣 X 的KK 方差矩陣， YX 和XY為矩陣 X 和矩陣 Y 的1K或者K1協(xié)方差矩陣。 通過上述的矩陣分塊， 可以計算出矩陣 11 YYYXXXXY B, p 為 矩陣 B 最大的特征值開平方根， 2 1 YYXX XY aa a p其中向量a為矩陣 B 最大 的特征值鎖對應(yīng)的特征向量。此處基于典型關(guān)聯(lián)分析多表型模型采用 Wiki lambda 2 1p和 Wiki

31、lambda 相關(guān)的 F 檢驗來獲取 P 值。其中 F 統(tǒng)計量為公式 - 18 - （4.28） 11 2 2 KnpK p F（4.28） 4.3.24.3.2 MultiPhenMultiPhen 模型模型 針對第四問，我們試圖從單個表型信號接近顯著的 SNP 中結(jié)合多個表型， 然后找出信號顯著的 SNP。在這里我們使用 MultiPhen 模型，從多個相關(guān)的表 型中提煉出具有強(qiáng)相關(guān)聯(lián)的 SNP。 在標(biāo)準(zhǔn)的多元 GWAS 方法分析中， 通常使用線性回歸對表型 Y 和基因型 X 進(jìn) 行相關(guān)的建模，其中我們令,其中 i 為第 i 個人，k 表示第 k 個 表型，令,其中，i 表示第 i 個人，

32、g 表示第 g 個 基因型。通過回歸的方法來檢測第 g 個基因型 SNP 和第 k 個表型之間的關(guān)系， 我們有如下公式（4.29） ： （4.29） 其中表示殘差， 而且假設(shè)服從正態(tài)分布。 而在 MultiPhen 模型中所采用 的方法是逆回歸的方法，具體來說就是把基因型 X 作為因變量，表型 Y 作為自 變量?；蛐蛿?shù)據(jù) X 是等位基因的數(shù)量，需要用到序數(shù)回歸來進(jìn)行相關(guān)的建模， 這里使用的是比例邏輯回歸（proportional odds logistic regression） ，這 個模型的類別概率定義如公式(4.30)： (4.30) 其中對于每一個 SNP，g=1,.,G 使用似然比

33、率檢驗（likelihood ratio test）去檢驗原假設(shè)。這個檢驗不用假設(shè) Hardy-Weinberg 平衡。 - 19 - 五、五、問題求解問題求解 5.15.1 問題一求解問題一求解 根據(jù)本題數(shù)據(jù)特征，每個位點均表示為兩個不同堿基的組合形式，并且 由于染色體的雙螺旋結(jié)構(gòu)，每種位點只有三種組合方式，如在位點 rs100015 位 置，不同樣本的編碼都是 T 和 C 的組合，有三種不同編碼方式 TT,TC 和 CC。因 此基于生物基因的加性效應(yīng) (Additive genetic effects6)，位點中的三種 組合分別可編碼為 0、1,、2，其中 1 代表雜合子基因，0 和 2

34、分別代表另外純 合子基因中的主要等位基因(major allele)與次要等位基因(minor allele)。 主要等位基因和次要等位基因主要由頻率決定，頻率較低的為次要等位基因。 比如，在題中給出的位點樣本中，前十列位點對應(yīng)的基因編碼方式表 5-1： 表 5-1 問題一基因編碼方式表 位點名稱次要/主要 等位基因 次要等位 基因（純合子） 雜合子主要等位 基因（純合子） rs3094315C/TCC = 2CT/TC = 1TT = 0 rs3131972T/CTT = 2CT/TC = 1CC = 0 rs3131969T/CTT = 2CT/TC = 1CC = 0 rs1048488

35、C/TCC = 2CT/TC = 1TT = 0 rs12562034A/GAA = 2AG/GA = 2GG = 0 rs12124819G/AGG = 2GA/AG = 2AA = 0 rs4040617G/AGG = 2GA/AG = 2AA = 0 rs2980300T/CTT = 2CT/TC = 1CC = 0 rs4970383T/GTT = 2TG/GT = 1GG = 0 rs4475691T/CTT = 2CT/TC = 1CC = 0 5.25.2 問題二求解問題二求解 分析題目可得本題主要過程是建立計算樣本中位點與樣本整體之間的關(guān)聯(lián) 程度，即獨立性程度，因此本題采用卡方

36、檢驗?zāi)Ｐ?。在進(jìn)行卡方檢驗之后，我 們得到了 9445 次檢驗的 P 值，然后作出相關(guān)的 Q-Q 圖。如圖 5-1，可以發(fā)現(xiàn)部 分 P 值的實際值與觀測值存在較大的偏差。在 GWAS 研究中，則認(rèn)為這個 SNP 位 點的觀測值的偏離是由這個 SNP 突變所產(chǎn)生的。這些檢驗說明可能存在致病性 的 SNP 位點。下文中我們選用適合的閾值，對 SNP 致病性位點作進(jìn)一步的篩選。 - 20 - 01234 0246 卡方檢驗卡方檢驗Q-Q圖圖 Expectedlog10p Observedlog10p 圖 5-1 問題二 卡方檢驗 Q-Q 圖 由于問題二中涉及到的位點有 9445 個，也就是說要進(jìn)行 9

37、445 次的卡方檢 驗，這就涉及到了多重檢驗的問題，因此需要對 P 值的閾值進(jìn)行調(diào)整或者對 P 值本身進(jìn)行校正，要不然會產(chǎn)生假陽性。我們在對 P 值進(jìn)行閾值調(diào)整的時候選 用了 T 05. 0 新的閾值，在對 P 值本身進(jìn)行校正方面，還選用了 BH 校正以及 BONF 校正。 從圖 5-2 中， 我們可以看出 BONF 校正相對于 BH 校正來說，更加的保守。 除了在圖 5-2 的結(jié)果以外，也有其他論文報道過 BH 進(jìn)行校正的時候，能夠盡可 能在減少假陽性的同時，避免假陰性增加太多，因此下文所有的 P 值校正均采 用 BH 校正。分別使用置換檢驗?zāi)Ｐ团c貝葉斯模型進(jìn)行結(jié)果的統(tǒng)計分析與檢驗， 由表

38、5-2，表 5-4，表 5-5 的結(jié)果對比表明這兩種模型得到的實驗結(jié)果很好地印 證了卡方檢驗的結(jié)果，因此最終選取與某種疾病最相關(guān)的位點為：rs2273298。 通過卡方檢驗?zāi)Ｐ偷玫降穆D圖如圖 5-2 所示： - 21 - 圖 5-2 問題二 卡方檢驗的曼哈頓圖 卡方檢驗?zāi)Ｐ偷玫降那?10 個與某種疾病相關(guān)的位點關(guān)聯(lián)程度如表 5-2 所 示： 表 5-2 問題二 位點關(guān)聯(lián)表 P 值 位點 卡方檢驗?zāi)Ｐ?調(diào)整前BH 調(diào)整BONF 調(diào)整 rs22732986.378e-80.00060240.0006024 rs9323728.028e-50.30130.7582 rs120362169.571

39、e-50.30130.904 rs28073450.00017380.34851 rs43916360.00024140.34851 rs75223440.000250.34851 rs94263060.00025830.34851 rs121339560.00060220.71091 rs115802180.0009640.9871 rs5903680.0015760.9871 置換檢驗?zāi)Ｐ筒捎玫膮?shù)如表 5-3 所示： 表 5-3 問題二 參數(shù)列表 參數(shù)表示參數(shù)值 10 1000000 0.0001 0.01 5 0.001 - 22 - 相關(guān)參數(shù)含義參考模型介紹。最終置換檢驗?zāi)Ｐ偷玫降膱D

40、 5-2。 置換檢驗?zāi)Ｐ偷玫降那?10 個與某種疾病相關(guān)的位點關(guān)聯(lián)程度如表 5-4 所 示： 表 5-4 問題二 置換檢驗的位點關(guān)聯(lián)表 期望 P 值 位點 置換檢驗?zāi)Ｐ?調(diào)整前調(diào)整后 rs22732981e-60.009445 rs9323727.6e-50.311685 rs120362169.9e-50.311685 rs43916360.00019780.4298824 rs94263060.00029820.4298824 rs75223440.00031340.4298824 rs28073450.00031860.4298824 rs115802180.00077810.86935

41、98 rs121339560.00082840.8693598 rs7074720.0011861 貝葉斯模型得到的前 10 個與某種疾病相關(guān)的位點關(guān)聯(lián)程度如表 5-5 所示： 表 5-5 問題二 貝葉斯模型的位點關(guān)聯(lián)表 log10_bf 位點 貝葉斯模型 rs22732984.51238 rs9323722.13662 rs43916362.12278 rs120362162.04464 rs75223442.03471 rs94263062.01492 rs28073451.9341 rs121339651.65489 rs115802181.57159 rs7074721.4864 5.

42、35.3 問題三求解問題三求解 根據(jù)問題三種數(shù)據(jù)的特征，基因與疾病的關(guān)聯(lián)性可由基因相關(guān)的位點的關(guān) 聯(lián)性表示，因此本題采用基于集合的基因檢驗?zāi)Ｐ团c VEGAS 模型來對樣本進(jìn)行 建模求解，最終得到與疾病關(guān)聯(lián)最有可能的基因為：gene_102、gene_55、 gene_217。 采用基于集合的基因檢驗?zāi)Ｐ团c VEGAS 模型對基因樣本進(jìn)行關(guān)聯(lián)性分析結(jié) 果如圖 5-3 所示： - 23 - Set-based test VEGAS 1 2 3 1 2 3 0100200300 基因基因 -log10(p) 圖 5-3 問題三 關(guān)聯(lián)性分析結(jié)果圖 其中兩種模型前 10 個（兩個模型前 10 排序相同）

43、與疾病相關(guān)的基因 BH 校 正 p 值如表 5-6 所示： 表 5-6 問題三 兩個模型的 BH 校正 P 值對比表 校正 P 值 基因 VEGAS 模 型 Set-based 檢 驗?zāi)Ｐ?gene_2170.00090.0149985 gene_550.001650.0149985 gene_1020.01840.04 gene_1690.07590.15 gene_1490.55740.53394 gene_220.79928570.7885714 gene_2930.79928750.7885714 gene_2350.81562500.822 gene_1500.87733330.82

44、2 gene_1330.92340.822 5.45.4 問題四求解問題四求解 針對問題中多個性狀的特征， 采用 mv-plink 模型和 MultiPhen 模型進(jìn)行問 題求解。閾值選取為 BH 校正過后，P 值小于 0.05。兩個模型均在樣本中保持較 好且較一致的效果， 最終得到與樣本中 10 個性狀有關(guān)聯(lián)的位點為： rs12746773。 兩個模型對應(yīng)的 BH 校正曼哈頓圖如圖 5-4 所示： - 24 - MultiPhen MV-plink 0 5 10 15 20 0 5 10 15 20 0250050007500 BP區(qū)間區(qū)間 -log10(p) 圖 5-4 問題四 BH 校正

45、曼哈頓 其中兩種模型前 10 個（兩個模型前 10 排序相同）與疾病相關(guān)的基因 BH 校 正 P 值如表 5-7 所示： 表 5-7 問題四 兩個模型的 BH 校正 P 值對比表 校正 P 值 位點 mv-plink校正P值 位點 MultiPhen rs127467732.868447e-2 0 rs127467738.306198e-2 1 rs24832755.408352e-1rs109167555.772251e-1 rs21856395.408352e-1rs101578355.772251e-1 rs109154235.408352e-1rs49088775.772251e-1

46、rs49088775.408352e-1rs121443755.772251e-1 rs121443755.408352e-1rs109167555.772251e-1 rs49088035.408352e-1rs112492485.772251e-1 rs6295245.408352e-1rs107991285.772251e-1 rs75495995.408352e-1rs115772625.772251e-1 rs109167555.408352e-1rs118114976.199957e-1 - 25 - 六、六、模型評價模型評價 本文針對問題中給出的數(shù)據(jù)的結(jié)構(gòu)特點，建立不同的模型進(jìn)行

47、數(shù)據(jù)分類， 并綜合運(yùn)用各種方法得到較好效果。同時為了滿足不同數(shù)據(jù)的需求，不斷對模 型進(jìn)行改進(jìn)提升，主要工作如下： 1. 針對本題并沒有獨立的數(shù)據(jù)集進(jìn)行驗證，問題二、問題三、問題四，我們 都給出了至少兩種的模型，并且通過這些不同模型的結(jié)果進(jìn)行相互的驗證。 2.針對單個 SNP 與疾病進(jìn)行關(guān)聯(lián)分析的時候，接近顯著的 SNP 往往會被忽略， 這里我們利用同一基因中 SNP 可能的連鎖不平衡關(guān)系，找出與疾病相關(guān)的 基因。并且對 VEGAS 模型中的基因列表以及基因與 SNP 的關(guān)系進(jìn)行精簡， 使得模型更加適合題目的應(yīng)用場景。 - 26 - 參考文獻(xiàn)參考文獻(xiàn) 1Benjamin Y,Hochberg Y

48、. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple TestingJ. Statist Society Series, 1995,57(2):289-300. 2Scott, L., Mohlke, K., Bonnycastle, L., Willer, C., Li, Y., Duren, W., Erdos, M., Stringham, H., Chines, P., Jackson, A., Prokunina-Olsson, L., Ding, C., Swift

49、, A., Narisu, N., Hu, T., Pruim, R., Xiao, R., Li, X., Conneely, K., Riebow, N., Sprau, A., Tong, M., White, P., Hetrick, K., Barnhart, M., Bark, C., Goldstein, J., Watkins, L., Xiang, F., Saramies, J., Buchanan, T., Watanabe, R., Valle, T., Kinnunen, L., Abecasis, G., Pugh, E., Doheny, K., Bergman, R., Tuomilehto, J., Collins, F., and Boehnke, M. A genome-wide association study of type 2 diabetes in finns detects multiple susceptibility variants. Science 316(5829), 13415 (2007). 3Guan,