焦會永 ，孫亞男 ，景曉溪 ，劉 京 ，江 麗 ，李彩霞 ，葉 健 ，劉 凡 ，黃艷梅 ，趙雯婷

（1.新鄉(xiāng)醫(yī)學(xué)院法醫(yī)學(xué)系，河南 新鄉(xiāng) 453003；2.公安部物證鑒定中心 北京市現(xiàn)場物證檢驗工程技術(shù)研究中心 現(xiàn)場物證溯源技術(shù)國家工程實驗室 法醫(yī)遺傳學(xué)公安部重點實驗室，北京 100038；3.濟寧醫(yī)學(xué)院法醫(yī)學(xué)與醫(yī)學(xué)檢驗學(xué)院，山東 濟寧 272067；4.中國科學(xué)院北京基因組研究所，北京 100029）

通過DNA進行人類身高、外貌等表型特征的預(yù)測，是目前法醫(yī)DNA領(lǐng)域的研究熱點之一[1-3]。身高特征由于比較容易識別、可以量化、成年后相對穩(wěn)定，成為描述和識別個體特征的重要指標(biāo)，對于案件偵查具有重要的意義。研究[4]表明，身高是典型的多基因遺傳性狀，受遺傳、環(huán)境等因素影響。其中，遺傳因素對身高的影響較環(huán)境因素更加明顯，占80%左右[5]。目前國內(nèi)外研究報道了大量身高相關(guān)基因、SNP位點及稀有變異等遺傳信息[6]。

利用法醫(yī)人類學(xué)方法對身高進行預(yù)測的研究取得了一定的進展[7-8]，然而可以應(yīng)用于實際辦案的身高預(yù)測技術(shù)方法仍然較少。在法醫(yī)DNA領(lǐng)域，目前構(gòu)建的身高預(yù)測模型主要是針對歐洲極高身高個體的定性研究，使用的主要方法有通過利用分類與回歸樹（classification and regression tree，CART）法、支持向量機（support vector machine，SVM）法及加權(quán)等位基因求和（weighted allele sums，WAS）等模型來預(yù)測身高[9]。AULCHENKO等[10]首次報道了采用WAS模型針對54個SNP位點進行荷蘭人群中5%的極高身高個體的定性預(yù)測研究，其預(yù)測準(zhǔn)確性的曲線下面積（area under curve，AUC）值為 0.65，比隨機預(yù)測的 0.5稍高。LIU等[9]基于歐源樣本的研究建立了180個SNP位點的WAS預(yù)測模型，極高身高個體的定性預(yù)測準(zhǔn)確性提升到了0.75。2014年，WOOD等[11]以包含253 288個歐洲個體的79個研究項目為基礎(chǔ)，通過Meta分析得到了697個與身高顯著關(guān)聯(lián)的SNP位點，是目前國際上身高研究領(lǐng)域影響力較大的一項研究成果。目前，針對亞洲人群的身高模型研究還較少。本研究擬依據(jù)WOOD等[11]研究中相關(guān)SNP位點構(gòu)建身高預(yù)測模型，評估其在中國山東漢族男性群體中的預(yù)測能力。

1 對象與方法

1.1 研究對象

樣本采集前對志愿者進行問卷調(diào)查，主要調(diào)查包括三代以內(nèi)親屬的民族、年齡及是否患有影響身高的相關(guān)疾病等問題，篩選出符合標(biāo)準(zhǔn)的志愿者進行樣本采集。納入標(biāo)準(zhǔn)：（1）三代以內(nèi)均為漢族且無異族通婚史；（2）年齡為 20～50 歲的男性；（3）身高≥180cm（高身高）或≤165cm（低身高）。排除標(biāo)準(zhǔn)：（1）曾接受過激素治療；（2）有先天性發(fā)育異?；蚝筇炱髻|(zhì)性畸形；（3）患有甲狀腺疾病、垂體疾病及腫瘤等；（4）患有身高缺陷性、家族遺傳性、精神性、影響骨代謝及其他影響生長發(fā)育的急慢性疾病。本研究共采集了59例山東漢族男性無關(guān)個體（表1）的血液樣本，分別用于全基因組關(guān)聯(lián)分析和二代測序，所有志愿者均簽署知情同意書。

表1 研究對象一般資料

身高由經(jīng)過專業(yè)培訓(xùn)的人員進行測量，采用統(tǒng)一的測量工具及測量方法。測量方法：被測量者呈立正姿勢，赤腳站立于身高測量計的底板上，足跟并攏，足尖分開約呈60°，上肢自然下垂，肩甲、骶骨及腳跟緊靠身高測量計的立柱，雙眼平視前方；測量者站在被測量者一側(cè)，調(diào)整身高測量計平板至適宜位置，記錄身高數(shù)據(jù)，精確至1cm。

1.2 方法

1.2.1 DNA提取與定量

使用 QIAamp DNA Blood Midi試劑盒（德國QIAGEN公司）提取所有樣本的全基因組DNA，采用NanoDrop 2000分光光度計（美國Thermo Scientific公司）檢測DNA濃度。

1.2.2 SNP分型檢測

全基因組關(guān)聯(lián)分析（Genome-wide association study，GWAS）研究：采用 Affymetrix SNP Array 6.0芯片（美國Thermo Fisher Scientific公司）對42例山東漢族男性樣本進行基因分型。用Sty或NSP酶對250ng DNA樣本進行消化，對消化后的小片段DNA進行連接反應(yīng)，然后用通用引物進行PCR擴增。擴增反應(yīng)完成后，對產(chǎn)物進行片段化處理并用生物素對片段進行標(biāo)記，再與芯片進行雜交。經(jīng)過16～18h的雜交反應(yīng)，采用GeneChipTMScanner 3000 7G（美國AB公司）檢測雜交信號熒光強度并讀取數(shù)據(jù)。

二代測序：采用HiSeq 4000測序平臺（美國Illumina公司）對17例山東漢族男性樣本進行全基因組測序，經(jīng)過BWA軟件（美國Illumina公司）[12]將測序短序列與參考基因組進行比對，確定短序列在基因組上的位置，使用SAMtools軟件（美國Illumina公司）[13]將這些短序列按一定順序排列并進行數(shù)據(jù)轉(zhuǎn)換，通過Picard軟件（美國博德研究所）進一步去除測序過程中產(chǎn)生的冗余信息及噪聲，使用GATK軟件（美國博德研究所）尋找SNP位點，最后使用ANNOVAR軟件[14-15]對SNP位點進行功能注釋。17例樣本的DNA測序檢測委托北京市理化分析測試中心及中國科學(xué)院北京基因組研究所進行。

1.2.3 數(shù)據(jù)處理及統(tǒng)計學(xué)分析

針對芯片檢測獲得的SNP分型結(jié)果采用GCTA軟件[16]得到樣本之間的遺傳相關(guān)矩陣（genetic relationship matrix，GRM），未檢測到有強親緣關(guān)系的樣本。對GRM進行奇異值分解（singular value decomposition，SVD）得到主要的特征向量（eigenvector），對主要的特征向量進行k-means聚類，從而進行人群分層評估，并將主要的特征向量在GWAS中作為協(xié)變量以校正可能存在的人群亞分層混淆因素。GWAS分析中用PLINK v1.9軟件（美國博德研究所）[17]的Logistic回歸檢測基因型-表型相關(guān)性，因變量為logit（P），其中P為低身高的概率，低身高設(shè)為1，高身高設(shè)為0；SNP分型及PLINK v1.9軟件[17]分析得到的三個主成分（C1、C2及C3）為自變量，自變量SNP基因型以效應(yīng)頻率等位基因數(shù)目賦值，即非效應(yīng)等位基因純合型賦值為0、雜合型賦值為1，效應(yīng)等位基因純合型賦值為2，挖掘與表型顯著相關(guān)的SNP位點（P≤0.05認(rèn)為與身高潛在關(guān)聯(lián)，P≤5×10-8認(rèn)為與身高達到全基因組水平顯著關(guān)聯(lián)[11]）。使用Mathematical軟件[18]對GWAS數(shù)據(jù)進行基因填補[5，19-20]，從得到的相應(yīng)SNP分型數(shù)據(jù)中篩選出與參考文獻[11]共有的SNP位點信息。

17例樣本的測序數(shù)據(jù)，使用PLINK v1.9軟件[17]進行主成分分析以去除人群分層，提取SNP位點分型信息，與上述GWAS分型結(jié)果進行合并分析。由于樣本量較少，這17例樣本未進行GWAS分析。

采用R v3.2.2軟件繪制GWAS分析Q-Q plot圖及曼哈頓圖，通過Q-Q plot圖直觀顯示觀測值與預(yù)測值之間的差異，判斷圖中的SNP位點分布是否合理，以排除潛在的偏差SNP位點，如分型錯誤；通過曼哈頓圖直觀顯示SNP位點與身高的關(guān)聯(lián)性。

1.3 身高預(yù)測模型建立

首先根據(jù)參考文獻[11]中報道的697個SNP位點效應(yīng)等位基因的Beta值（其值大小和正負與該等位基因?qū)ι砀叩淖饔贸收嚓P(guān)），獲得每個樣本的WAS值，即：

其中，bm為第m個SNP位點的效應(yīng)等位基因的Beta值，xm為每個樣本第m個SNP位點的效應(yīng)等位基因數(shù)目。繼而依照以二進制身高為因變量（y）的Logistic回歸模型：

計算每個樣本為高身高或低身高的概率，即：

然后參考LIU等[9]的研究方法，使用SPSS 19.0軟件根據(jù)樣本數(shù)據(jù)繪制受試者工作特征（receiver operating characteristic，ROC）曲線并計算其曲線下面積（area under curve，AUC），采用 AUC 值評價預(yù)測模型的準(zhǔn)確性，AUC值分布為0.5～1.0，0.5表示沒有預(yù)測能力，數(shù)值越接近1.0表示模型的預(yù)測能力越準(zhǔn)確[21]。

將芯片篩選結(jié)果與二代測序提取結(jié)果合并，篩選出共有的SNP位點，用于身高預(yù)測模型準(zhǔn)確性評價。針對上述樣本數(shù)據(jù)，用Excel軟件（美國Microsoft公司）繪制散點圖，使用SPSS 19.0軟件（美國IBM公司）對模型所用SNP位點的Beta值在本研究與參考文獻[9]中的差異進行卡方檢驗，并使用R v3.2.2軟件繪制高身高、低身高分布密度曲線。

2 結(jié) 果

2.1 GWAS結(jié)果

采用Affymetrix SNP Array 6.0芯片進行檢測的42例樣本共檢測了909622個SNP位點，使用PLINK v1.9軟件對檢測的全部SNP位點進行質(zhì)量控制：（1）剔除分型成功率小于90%的樣本2個；（2）剔除分型成功率小于90%的SNP位點90 039個；（3）剔除效應(yīng)等位基因最小等位基因頻率（minor allele frequency，MAF）小于0.01的SNP位點161595個。最終共得到40例樣本的657988個SNP位點用于身高關(guān)聯(lián)分析。

二代測序17例樣本的質(zhì)量控制：（1）剔除分型成功率小于90%的染色體位點5489036個；（2）剔除較小等位基因頻率小于0.05的染色體位點1047 604個。最終共獲得17例樣本的3123513個SNP位點。

通過PLINK v1.9軟件中的主成分分析（principal component analysis，PCA）模塊將個體的高維度（維度等于該個體檢測的SNP位點數(shù)目）全基因組SNP位點分型數(shù)據(jù)降到低維度（維度不大于個體數(shù)目）。分析得到的前三個主成分（C1、C2及C3）累積貢獻率〉85%，故以這三個主成分作為協(xié)變量加入到關(guān)聯(lián)分析模型中，進一步校正可能的人群分層[22]，17例樣本和40例樣本的分析均未發(fā)現(xiàn)人群分層，Q-Q plot圖可見樣本存在較高的偏度和峰度（圖1A）。通過GWAS樣本分析，繪制曼哈頓圖顯示無達到全基因組水平顯著關(guān)聯(lián) SNP位點（P≤5×10-8），未發(fā)現(xiàn)與身高顯著相關(guān)的SNP 位點（圖 1B）。

2.2 用WAS法構(gòu)建中國漢族身高預(yù)測模型

40例樣本的分型結(jié)果與參考文獻[11]中697個SNP位點的共有位點為186個，經(jīng)基因填補[23]后達到了612個。對二代測序的17例樣本分型結(jié)果進行提取，得到595個SNP位點。兩者統(tǒng)一整合后篩選出547個共有的SNP位點，并將效應(yīng)等位基因（頻率較小的等位基因）與參考文獻[11]進行比對，根據(jù)效應(yīng)等位基因是否一致調(diào)整Beta值方向，進而獲得適用于構(gòu)建本研究身高預(yù)測模型的Beta值，構(gòu)建WAS預(yù)測模型。利用身高預(yù)測模型繪制ROC曲線并計算其AUC值為0.67（95%置信區(qū)間為0.53～0.90）。根據(jù)1.3節(jié)的公式獲得WAS值繪制密度曲線（圖2），提示此模型對高身高（30例）及低身高（27例）的預(yù)測有一定區(qū)分度。

圖1 身高關(guān)聯(lián)研究Logistic回歸分析結(jié)果

圖2 57例樣本身高的WAS密度曲線

身高相關(guān)SNP位點對人類身高的影響具有人群差異。Beta值方向分布散點圖（以兩組Beta值交點為原點，處在第一、三象限的Beta值方向相同，第二、四象限的Beta值方向相反）顯示，本研究與WOOD等[11]研究人群之間Beta值方向分布有較大差異（圖3），提示身高相關(guān)SNP位點存在種族特異性。針對歐洲及中國漢族547個SNP位點對應(yīng)Beta值在高身高及低身高水平的分布（表3）進行卡方檢驗，χ2=6.934，雙側(cè)精確P=0.008，單側(cè)精確P=0.040，兩者Beta值之間差異具有統(tǒng)計學(xué)意義，進一步提示種族特異性的存在。

圖3 547個SNP位點對應(yīng)Beta值在本研究及參考文獻[11]中的分布

表3 547個SNP位點分布情況 （個）

3 討 論

身高是由遺傳和環(huán)境因素共同影響的復(fù)雜多基因性狀，隨著新一代測序技術(shù)及GWAS的推廣，大量身高相關(guān)SNP位點被發(fā)現(xiàn)。近幾年，國外不少研究開始嘗試基于身高相關(guān)SNP位點來預(yù)測人類身高[9]。對人類身高的預(yù)測研究在臨床診斷、法庭科學(xué)等多個領(lǐng)域都有潛在的應(yīng)用價值。在臨床方面，如果能從DNA水平預(yù)測未來身高發(fā)育方面的異常，有助于醫(yī)生盡早診斷先天性巨人癥及生長限制類疾病并及時進行干預(yù)[9]。在法醫(yī)學(xué)實踐方面，通過DNA數(shù)據(jù)預(yù)測個體的身高信息，可以為案件偵查提供線索。

由于國內(nèi)尚未發(fā)現(xiàn)足夠數(shù)量的中國漢族人群身高顯著相關(guān)SNP位點，因此無法直接構(gòu)建用于預(yù)測中國漢族人群身高的預(yù)測模型。本研究選取了當(dāng)前國際上在身高預(yù)測方面效果較好的歐洲高加索人群身高顯著相關(guān)SNP位點，構(gòu)建預(yù)測模型，用于評估其預(yù)測中國漢族人群身高的準(zhǔn)確性。評價標(biāo)準(zhǔn)即為ROC曲線的AUC。ROC曲線源于信號探測理論，最初用于描述信號與噪聲之間的關(guān)系，以評價雷達性能[24]，隨著ROC分析應(yīng)用于醫(yī)學(xué)診斷，日益受到重視，到目前已經(jīng)發(fā)展成為一種分析方法。ROC曲線采用二分法進行分析，以真陽性率（靈敏度）為縱坐標(biāo)、假陽性率（特異性）為橫坐標(biāo)，AUC越大，判斷準(zhǔn)確性就越大[25]。目前，國內(nèi)外人類身高預(yù)測模型較少，WAS即為其中之一且預(yù)測準(zhǔn)確性較高[9]。AULCHENKO等[10]首次將ROC曲線應(yīng)用于身高預(yù)測，根據(jù)每例樣本的WAS值繪制ROC曲線，構(gòu)建身高預(yù)測模型，在僅有54個身高相關(guān)SNP位點的情況下，AUC值達到了0.65。隨后，LIU等[9]根據(jù)180個身高相關(guān)SNP位點計算每例樣本W(wǎng)AS值并進行了ROC分析，其AUC值達到了0.75（95%置信區(qū)間為 0.72～0.79）。

目前，中國人群身高相關(guān)遺傳位點研究還處于起步階段，需要進行人群遺傳多態(tài)性分析、人群身高模型構(gòu)建、極端身高個體定性、個體身高預(yù)測等多個階段的研究，實現(xiàn)對實際案件生物檢材來源人身高的精細推測還有很長的路要走。本研究初步構(gòu)建了可用于中國漢族人群身高預(yù)測的計算模型并進行效果評估，借鑒AULCHENKO等[10]研究方法來評估此模型預(yù)測中國漢族人群身高的準(zhǔn)確性，為后續(xù)研究提供了數(shù)據(jù)支撐。本研究首先應(yīng)用Affymetrix Array SNP 6.0芯片對40例中國漢族男性樣本進行GWAS研究，僅發(fā)現(xiàn)了與身高潛在相關(guān)的SNP位點，但在全基因組水平未發(fā)現(xiàn)與身高顯著相關(guān)的SNP位點。可能由以下原因?qū)е拢海?）樣本量過少；（2）大部分與身高相關(guān)的SNP位點效應(yīng)較小，很難被檢測到。本研究采用WAS分析法獲得的AUC值為0.67（95%置信區(qū)間為0.53～0.90）。與AULCHENKO等[10]研究相比，本研究在使用的SNP位點數(shù)量大大提升的條件下，預(yù)測準(zhǔn)確性稍高，但仍低于LIU等[9]在歐洲人群中的研究。另外，本研究身高預(yù)測模型所用的547個身高相關(guān)SNP位點是基于歐洲人群研究得出的，而卡方檢驗結(jié)果提示身高相關(guān)SNP位點存在種族差異性[26]，推測本研究預(yù)測模型AUC值比LIU等[9]低的主要原因，是因為在歐洲人群發(fā)現(xiàn)的身高相關(guān)SNP對中國漢族人群并不完全適用。

在后續(xù)的研究工作中，需要關(guān)注身高遺傳學(xué)研究的最新進展，增加用于篩選驗證的相關(guān)SNP位點；加大人群驗證樣本量，不斷挖掘適用于中國漢族人群的身高相關(guān)SNP位點；改善計算方法體系，繼續(xù)優(yōu)化模型算法。

[1]KAYSER M.Forensic DNA Phenotyping：Predicting human appearance from crime scene material for investigative purposes[J].Forensic Sci Int Genet，2015，18：33-48.

[2]TURKHEIMER E.Genetic Prediction[J].Hastings Cent Rep，2015，45（S5）：32-38.

[3]WALSH S， LIU F， WOLLSTEIN A， et al.The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA[J].Forensic Sci Int Genet，2013，7（1）：98-115.

[4]楊應(yīng)忠，王亞平，馬蘭，等.中國漢族高原肺水腫易感基因的全基因組關(guān)聯(lián)研究[J].遺傳，2013，35（11）：1291-1299.

[5]DU M，AUER P L，JIAO S，et al.Whole-exome imputation of sequence variants identified two novel alleles associated with adult body height in African Americans[J].Hum Mol Genet，2014，23（24）：6607-6615.

[6]CHO Y S， GO M J， KIM Y J， et al.A largescale genome-wide association study of Asian populations uncovers genetic factors influencing eight quantitative traits[J].Nat Genet，2009，41（5）：527-534.

[7]THODBERG H H，JENNI O G，CAFLISCH J， et al.Prediction of adult height based on automated determination of bone age[J].J Clin Endocrinol Metab，2009，94（12）：4868-4874.

[8]MARTIN D D， SCHITTENHELM J， THODBERG H H.Validation of adult height prediction based on automated bone age determination in the Paris Longitudinal Study of healthy children[J].Pediatr Radiol，2016，46（2）：263-269.

[9]LIU F， HENDRIKS A E，RALF A， et al.Common DNA variants predict tall stature in Europeans[J].Hum Genet，2014，133（5）：587-597.

[10]AULCHENKO Y S，STRUCHALIN M V， BELONOGOVA N M，et al.Predicting human height by Victorian and genomic methods[J].Eur J Hum Genet，2009，17（8）：1070-1075.

[11]WOOD A R， ESKO T， YANG J， et al.Defining the role of common variation in the genomic and biological architecture of adult human height[J].Nat Genet，2014，46（11）：1173-1186.

[12]LI H，DURBIN R.Fast and accurate short read alignment with Burrows-Wheeler transform[J].Bioinformatics，2009，25（14）：1754-1760.

[13]LI H.A statistical framework for SNP calling，mutation discovery，association mapping and population genetical parameter estimation from sequencing data[J].Bioinformatics，2011，27（21）：2987-2993.

[14]WANG K，LI M，HAKONARSON H.ANNOVAR：functional annotation of genetic variants from high-throughput sequencing data[J].Nucleic Acids Res，2010，38（16）：e164.

[15]YANG H，WANG K.Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR[J].Nat Protoc，2015，10（10）：1556-1566.

[16]YANG J， LEE S H，GODDARD M E，et al.GCTA：a tool for genome-wide complex trait analysis[J].Am J Hum Genet，2011，88（1）：76-82.

[17]CHANG C C，CHOW C C，TELLIER L C，et al.Second-generation PLINK：rising to the challenge of larger and richer datasets[J].Gigascience，2015，4：7.

[18]LIU E Y，LI M，WANG W，et al.MaCH-admix：genotype imputation for admixed populations[J].Genet Epidemiol，2013，37（1）：25-37.

[19]ESTRADA K，KRAWCZAK M，SCHREIBER S，et al.A genome-wide association study of northwestern Europeans involves the C-type natriuretic peptide signaling pathway in the etiology of human height variation[J].Hum Mol Genet，2009，18（18）：3516-3524.

[20]LI Y，WILLER C，SANNA S，et al.Genotype imputation[J].Annu Rev Genomics Hum Genet，2009，10：387-406.

[21]陳衛(wèi)中，倪宗瓚，潘曉平，等.用ROC曲線確定最佳臨界點和可疑值范圍[J].現(xiàn)代預(yù)防醫(yī)學(xué)，2005，32（7）：729-731.

[22]PRICE A L，PATTERSON N J，PLENGE R M，et al.Principal components analysis corrects for stratification in genome-wide association studies[J].Nat Genet，2006，38（8）：904-909.

[23]何桑，丁向東，張勤.基因型填充方法介紹及比較[J].中國畜牧雜志，2013，49（23）：95-100.

[24]SWETS J A.Measuring the accuracy of diagnostic systems[J].Science，1988，240（4857）：1285-1293.

[25]王敬瀚.ROC曲線在臨床醫(yī)學(xué)診斷實驗中的應(yīng)用[J].中華高血壓雜志，2008，16（2）：175-177.

[26]LEI S F， TAN L J， LIU X G， et al.Genome-wide association study identifies two novel loci containing FLNB and SBF2 genes underlying stature variation[J].Hum Mol Genet，2009，18（9）：1661-1669.