孫昌春 許繼臣 江 麗 郭曉媛1， 趙雯婷 葉 健* 李彩霞*

（1）山西醫(yī)科大學(xué)法醫(yī)學(xué)院，太原 030001；2）公安部物證鑒定中心，法醫(yī)遺傳學(xué)公安部重點實驗室，現(xiàn)場物證溯源技術(shù)國家工程實驗室，北京 100038）

飲食、氣候等環(huán)境因素可能會影響人的表觀遺傳學(xué)結(jié)構(gòu)［1?3］。DNA甲基化是一種重要的表觀遺傳標(biāo)記，具有遺傳穩(wěn)定、含量豐富、隨齡變化等特點［4］。伴隨高通量的DNA甲基化數(shù)據(jù)的不斷涌現(xiàn)，大量與腫瘤等疾病發(fā)生相關(guān)的DNA 甲基化位點被發(fā)現(xiàn)［5］。在法醫(yī)學(xué)領(lǐng)域，DNA 甲基化已被用于年齡推斷［6?9］、組織屬性判別［10］、同卵雙胞胎的鑒別［11?12］等?；跉W美等國外人群的研究證明DNA甲基化在族群地域間存在差異［13?14］。

漢族是世界上人口最多的民族。研究表明，漢族人群具有混合特征［15］，呈現(xiàn)明顯的南北分化，在基于STR［16］、線粒體DNA（mtDNA）［17?18］、Y染色體SNP［15］、常染色體SNP［19］等遺傳標(biāo)記的研究均已表明漢族內(nèi)部存在明顯的南北遺傳差異。北方人群由于受到來自中亞和歐洲遺傳成分的影響，呈現(xiàn)東西走向的變化趨勢，南北方人群遺傳的差異以秦嶺淮河到長江為地理分界［20］。漢族人群與當(dāng)?shù)厣贁?shù)民族之間的遺傳差異小于南北方漢族之間的遺傳差異［21］。已有研究多局限在DNA遺傳標(biāo)記層面。中國不同地域的氣候、飲食等環(huán)境因素有很大差異，目前缺乏不同地域人群之間是否存在表觀遺傳學(xué)差異相關(guān)研究。本文應(yīng)用表觀基因組關(guān)聯(lián)分析（epigenome?wide association study，EWAS）技術(shù)研究了中國南北方漢族人群之間DNA甲基化的差異，并結(jié)合機器學(xué)習(xí)等算法構(gòu)建了DNA 甲基化南北方漢族人群推斷算法。

1 材料與方法

1.1 樣本信息

本研究使用的483 份漢族男性DNA 樣本來源于國家科技資源共享服務(wù)平臺計劃項目，使用Illumina 的甲基化芯片Infinium MethylationEPIC BeadChip（簡稱850K芯片）檢測（科技部備份號：*BF2020121803316），850K 芯片可以檢測整個基因組860 000 個CpG 位點，覆蓋CpG 島、啟動子區(qū)、編碼區(qū)、開放染色質(zhì)和增強子區(qū)域。樣本的數(shù)量和年齡、性別信息分別見表1及附件表S1。本研究通過公安部物證鑒定中心倫理委員會審查（編號：2017?001），所有參與者均簽署了書面知情同意書。

1.2 質(zhì)量控制

使用R 軟件的ChAMP 包對低質(zhì)量數(shù)據(jù)進行預(yù)處理，根據(jù)以下原則過濾探針：檢測到原始探針的P 值大于0.01 的甲基化位點；在大于等于5%的樣本的beads 數(shù)目小于3 的探針；非CpG 探針；是SNP 或探針覆蓋區(qū)域內(nèi)存在SNP 的探針［22］；被鑒定為 交 叉反 應(yīng) 的探 針（cross?reactive or multi?hit probes）［23］，以及性染色體上的探針。分析得到每個探針位點的甲基化β 值，然后應(yīng)用BMIQ（beta?mixture quantile）方法進行β 值的歸一化處理。使用奇異值分解（singular value decomposition，SVD）方法檢測是否存在批次效應(yīng)與甲基化水平的關(guān)系［24］。

本研究中樣本的DNA 提取自外周血，外周血中不同的細胞類型DNA 甲基化特征都不相同，細胞類型組成是EWAS 分析中的一個潛在混雜因素。使用ReFACTor 算法［25］計算不同細胞類型組成的主成分并在GLINT 軟件下游分析時將其添加為協(xié)變量，從而減少細胞類型對EWAS 分析結(jié)果的影響。參考Teschendorff 等［26］研究，設(shè)置參數(shù)K=7（7 種細胞類型，嗜酸性粒細胞、中性粒細胞、單核細胞、B 細胞、NK 細胞、CD8+ T 和CD4+ T細胞）。

1.3 位點篩選及評估

通過EWAS 結(jié)合機器學(xué)習(xí)Lasso 回歸的方法篩選南北方漢族人群的差異甲基化位點。GLINT 軟件［27?30］基于Python2.7 開發(fā)，可以實現(xiàn)組織異質(zhì)性校正和EWAS分析兩種功能?；谏鲜鲑|(zhì)量控制后的DNA 甲基化矩陣文件，首先對南方漢族和北方漢族人群進行二分類編碼，南方漢族編碼為1，包括江西漢族和四川漢族，北方漢族編碼為0，包括山東漢族、山西漢族和河南漢族，并以此編碼用作表型文件，通過GLINT 軟件的數(shù)據(jù)管理功能（data management）檢測并刪除異常值。然后使用ReFACTor 算法校正組織的異質(zhì)性（adjusting for tissue heterogeneity），最后對生成的數(shù)據(jù)文件進行EWAS分析，刪除染色體X和Y中的位點，刪除非特異性的位點［31］。年齡和細胞類型用作協(xié)變量。

使用Lasso 分析方法進一步篩選位點，使用R（version 4.0.2）軟件的glmnet 包建立Lasso?logistic回歸模型。該模型是通過構(gòu)造懲罰函數(shù)實現(xiàn)變量選擇和參數(shù)估計，通過將其回歸系數(shù)設(shè)置為0的方式剔除呈現(xiàn)共線性或者與因變量沒有相關(guān)性的冗余變量［32?33］，選擇對因變量影響較大的自變量并計算出相應(yīng)的回歸系數(shù)，最終得到一個預(yù)測模型。

使用R 軟件基于每個DNA 甲基化差異位點在中國南北方漢族群體的甲基化水平表達特征進行可視化展示。使用R軟件softmaxreg包構(gòu)建多元邏輯回歸模型，通過十折交叉驗證的方法評估模型的準(zhǔn)確性。其中多元邏輯回歸使用caret 包中的createDataPartition函數(shù)從訓(xùn)練集的每個標(biāo)記人群中隨機采樣70%個體構(gòu)建模型，其余30%個體測試模型準(zhǔn)確性，使用confusionMatrix函數(shù)對上述模型進行評價；十折交叉驗證則是使用caret 包中的createFolds函數(shù)對數(shù)據(jù)集進行劃分，將數(shù)據(jù)集分成10 份，輪流將其中9 份作為訓(xùn)練集構(gòu)建模型，1 份作為測試集來測試模型的準(zhǔn)確性。為了使結(jié)果更加精確，每種算法運行10 次，使用seed 函數(shù)設(shè)定隨機數(shù)的初始值，產(chǎn)生不同的樣本組合，10 次結(jié)果的均值作為對算法精度的估計。模型評價指標(biāo)包括Kappa 系 數(shù)、 靈 敏 性（sensitivity）、 特 異 性（specificity）、陽性預(yù)測值（PPV）和陰性預(yù)測值（NPV）。

2 結(jié)果與分析

2.1 73個CpG位點的篩選

483 例漢族男性樣本數(shù)據(jù)的EWAS 分析結(jié)果見圖1。挑選355 個差異性顯著的CpG 位點（P ＜1×10-6）。使用Lasso回歸進一步篩選位點剔除冗余變量。Lasso 回歸的特點是在擬合廣義線性模型的同時進行變量篩選和復(fù)雜度調(diào)整。變量篩選是指有選擇的把變量放入模型從而得到更好的性能參數(shù)。復(fù)雜度調(diào)整的程度由參數(shù)λ 來控制，λ 越大對變量較多的線性模型的懲罰力度就越大，從而精簡變量，結(jié)果見圖2a、b。根據(jù)圖2a曲線最低點確定懲罰值λ，在圖2b的相應(yīng)懲罰值的位置確定出模型最后所納入的變量，最終篩選出73 個CpG 位點，位點信息具體見表2。

Fig.1 Manhattan diagram of EWAS analysis results

Fig.2 Characteristic variable screening based on Lasso regressionThe figure shows the process of selecting the most appropriate value of parameter λ in the Lasso model by cross?validation. (a) Lasso regression cross?validation of the optimal parameter atlas.(b)Sites in the model of regression coefficient.

Continued to Table 2

Continued to Table 2

2.2 位點的評估

73 個CpG 位點在南北方漢族群體的甲基化水平表達特征如圖3 所示，藍色代表北方漢族群體，紅色代表南方漢族群體，圖中的白點是中位數(shù)，小提琴圖中的黑色粗條范圍是上下四分位數(shù)，外部形狀即為分布密度?？梢钥闯鲞@些甲基化位點在兩個群體中具有不同程度的差異。73 個CpG 位點的南北方漢族群體的主成分分析結(jié)果見圖4。圖4a可以看到，前兩個主成分解釋變異的23%，基本可以區(qū)分南北方漢族人群。圖4b 中不同漢族人群用不同顏色標(biāo)注，南北方漢族人群較各亞人群之間的差異更大，北方人群中的河南漢族相對居中。絕大部分的江西漢族和四川漢族與山西漢族和山東漢族可以很好地區(qū)分開來。

Fig.3 Violin diagram of methylation levels of 73 CpG sites in northern and southern Han populations

Fig.4 Principal component analysis of northern and southern Han populations using 73 CpG sites(a) Principal component analysis of northern and southern Han populations. (b) Principal component analysis of Han populations of different province.

模型構(gòu)建與評估10 次重復(fù)的結(jié)果具體見表3，最終多元邏輯回歸30%測試集的預(yù)測準(zhǔn)確率為99.03%。Kappa系數(shù)均值為0.979 6；靈敏性是兩個相關(guān)個體正確歸類為相關(guān)的概率，特異性反之，分別為0.990 7 和0.989 5；陽性預(yù)測值是判為相關(guān)個體中有多少是真的相關(guān)，陰性預(yù)測值反之，結(jié)果分別為0.993 1 和0.986 2。10 次十折交叉驗證的結(jié)果均在98%以上，最終平均準(zhǔn)確率為98.79%，其余各項模型預(yù)測性能指標(biāo)均超過0.95。

3 討 論

環(huán)境因素可通過表觀遺傳機制，如誘導(dǎo)DNA甲基化模式改變等，在不改變DNA 序列前提下，改變基因表達，從而引發(fā)表型變化。法醫(yī)學(xué)研究已經(jīng)發(fā)現(xiàn)大量不同組織之間、不同年齡之間存在差異的DNA 甲基化位點，并構(gòu)建預(yù)測模型［6?10］。不同地域的氣候、飲食等環(huán)境因素存在差異，環(huán)境通過影響DNA 甲基化水平進而影響基因表達，使得生物表型發(fā)生變化，進而造成不同地域人群間的差異。Fraser 等［13］研究北歐（CEU）和西非（YRI）多個家系的DNA 甲基化數(shù)據(jù)，發(fā)現(xiàn)在族群內(nèi)和族群間存在顯著的甲基化差異，這種差異可能來自等位基因突變、上位效應(yīng)，以及基因跟環(huán)境之間的相互作用，還發(fā)現(xiàn)在轉(zhuǎn)錄起始位點附近的DNA 甲基化有群體特異性。Yuan 等［14］研究了509 份胎盤的450K甲基化芯片數(shù)據(jù)，基于甲基化位點進行亞洲、非洲、及高加索人群的推斷，準(zhǔn)確性為0.938。東亞是全球人口最多的區(qū)域，占全球人口的22%，漢族是東亞主體民族，大量基因組學(xué)研究發(fā)現(xiàn)漢族內(nèi)部存在明顯的南北遺傳差異，然而目前缺乏針對不同地域漢族人群的表觀遺傳學(xué)差異研究。

本文研究了483份漢族個體的甲基化數(shù)據(jù)，雖然南方人群沒有使用廣東廣西樣本，而是相對靠北的江西和四川人群樣本，仍然研究發(fā)現(xiàn)了南北方漢族人群之間的甲基化差異，并最終篩選出73 個CpG位點。圖4a、b展示了不同地域人群間的甲基化差異，這些差異可能歸因于等位基因頻率的差異以及上位效應(yīng)或基因與環(huán)境的相互作用［13］。地理學(xué)第一和第二定律指出地物之間的距離越近，相關(guān)性越大；空間隔離又造成地物之間的空間異質(zhì)性［34?36］。從圖4b 可以看出河南漢族居中，可能是因為河南在地理位置上屬于華中地區(qū)，地理位置接近的人群之間的基因交流頻繁，遺傳距離接近，遺傳差異度就小。對于距離相近的人群往往需要更多的位點和更大的參考人群數(shù)據(jù)才可以實現(xiàn)精確區(qū)分。模型預(yù)測性能驗證采用多元邏輯回歸隨機抽取70%個體構(gòu)建模型，30%個體進行測試以及十折交叉驗證的方法。未來可增加樣本量進一步驗證這一組甲基化位點的人群區(qū)分準(zhǔn)確性。

生物體的一些性狀是由多基因共同控制的，其測量值可以用連續(xù)的數(shù)量進行表示，這些性狀統(tǒng)稱為數(shù)量性狀（quantitative trait，QT）。遺傳變異可以調(diào)節(jié)多種機體代謝功能，包括DNA 甲基化在表觀遺傳學(xué)中所表現(xiàn)出的基因表達調(diào)節(jié)作用。目前比較重要的一種研究方法稱為數(shù)量性狀定位分析（quantitative trait locus，QTL），即定位控制數(shù)量性狀的基因位點在基因組中的位置。DNA 甲基化數(shù)量性狀位點分析（DNA methylation quantitative trait locus，meQTLs）是以DNA甲基化作為數(shù)量性狀，探討DNA 甲基化水平與基因表達之間關(guān)系［37?38］。除了年齡、性別以及疾病等環(huán)境因素會影響個體的甲基化水平，DNA 序列也會對甲基化產(chǎn)生影響，尤其是特異性SNP 的改變。有研究發(fā)現(xiàn)非洲和歐洲祖先群體之間70%的差異甲基化位點至少與一個meQTL相關(guān)［39］，這表明很大一部分的DNA甲基化的群體差異可能主要是由DNA序列變異引起。而從表觀遺傳角度入手實現(xiàn)地域人群的區(qū)分關(guān)鍵在于尋找獨立于SNP 的甲基化位點，從而作為STR、SNP等基因組遺傳標(biāo)記的有效補充。

近距離人群的區(qū)分和推斷是法醫(yī)DNA 領(lǐng)域的難點之一。目前已報道的AISNP 體系初步實現(xiàn)了東亞南北方人群的區(qū)分，未來結(jié)合表觀遺傳標(biāo)記，有望實現(xiàn)更加精細的人群區(qū)分。本文探索研究了利用甲基化進行不同地域人群推斷的可行性，研究表明南北方漢族人群之間存在表觀遺傳差異，未來需進一步增加人群數(shù)據(jù)獲取更多地域相關(guān)甲基化位點，與族群相關(guān)SNP 位點配合使用實現(xiàn)東亞人群的精細區(qū)分。

附件 PIBB_20210091_Table_S1.xlsx 見本文網(wǎng)絡(luò)版（http://www.pibb.ac.cn或http://www.cnki.net）。

Table 1 Size of samples used in this study

Table 2 Information of the 73 CpGs

Table 3 The performance statistics of the prediction model