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

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

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

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

1 材料與方法

1.1 樣本信息

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

1.2 質(zhì)量控制

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

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

1.3 位點(diǎn)篩選及評(píng)估

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

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

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

2 結(jié)果與分析

2.1 73個(gè)CpG位點(diǎn)的篩選

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

73 個(gè)CpG 位點(diǎn)在南北方漢族群體的甲基化水平表達(dá)特征如圖3 所示，藍(lán)色代表北方漢族群體，紅色代表南方漢族群體，圖中的白點(diǎn)是中位數(shù)，小提琴圖中的黑色粗條范圍是上下四分位數(shù)，外部形狀即為分布密度?？梢钥闯鲞@些甲基化位點(diǎn)在兩個(gè)群體中具有不同程度的差異。73 個(gè)CpG 位點(diǎn)的南北方漢族群體的主成分分析結(jié)果見圖4。圖4a可以看到，前兩個(gè)主成分解釋變異的23%，基本可以區(qū)分南北方漢族人群。圖4b 中不同漢族人群用不同顏色標(biāo)注，南北方漢族人群較各亞人群之間的差異更大，北方人群中的河南漢族相對(duì)居中。絕大部分的江西漢族和四川漢族與山西漢族和山東漢族可以很好地區(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)建與評(píng)估10 次重復(fù)的結(jié)果具體見表3，最終多元邏輯回歸30%測(cè)試集的預(yù)測(cè)準(zhǔn)確率為99.03%。Kappa系數(shù)均值為0.979 6；靈敏性是兩個(gè)相關(guān)個(gè)體正確歸類為相關(guān)的概率，特異性反之，分別為0.990 7 和0.989 5；陽性預(yù)測(cè)值是判為相關(guān)個(gè)體中有多少是真的相關(guān)，陰性預(yù)測(cè)值反之，結(jié)果分別為0.993 1 和0.986 2。10 次十折交叉驗(yàn)證的結(jié)果均在98%以上，最終平均準(zhǔn)確率為98.79%，其余各項(xiàng)模型預(yù)測(cè)性能指標(biāo)均超過0.95。

3 討 論

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

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

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

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