中國南北方漢族人群DNA甲基化表觀遺傳差異研究*

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

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

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

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

1 材料與方法

1.1 樣本信息

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

1.2 質量控制

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

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

1.3 位點篩選及評估

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

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

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

2 結果與分析

2.1 73個CpG位點的篩選

483 例漢族男性樣本數據的EWAS 分析結果見圖1。挑選355 個差異性顯著的CpG 位點（P ＜1×10-6）。使用Lasso回歸進一步篩選位點剔除冗余變量。Lasso 回歸的特點是在擬合廣義線性模型的同時進行變量篩選和復雜度調整。變量篩選是指有選擇的把變量放入模型從而得到更好的性能參數。復雜度調整的程度由參數λ 來控制，λ 越大對變量較多的線性模型的懲罰力度就越大，從而精簡變量，結果見圖2a、b。根據圖2a曲線最低點確定懲罰值λ，在圖2b的相應懲罰值的位置確定出模型最后所納入的變量，最終篩選出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 所示，藍色代表北方漢族群體，紅色代表南方漢族群體，圖中的白點是中位數，小提琴圖中的黑色粗條范圍是上下四分位數，外部形狀即為分布密度。可以看出這些甲基化位點在兩個群體中具有不同程度的差異。73 個CpG 位點的南北方漢族群體的主成分分析結果見圖4。圖4a可以看到，前兩個主成分解釋變異的23%，基本可以區(qū)分南北方漢族人群。圖4b 中不同漢族人群用不同顏色標注，南北方漢族人群較各亞人群之間的差異更大，北方人群中的河南漢族相對居中。絕大部分的江西漢族和四川漢族與山西漢族和山東漢族可以很好地區(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.

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

3 討 論

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

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

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

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

附件 PIBB_20210091_Table_S1.xlsx 見本文網絡版（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