秦 婕，王嘉博，益西康珠，鐘金城

(1.青藏高原動物遺傳資源保護與利用教育部重點實驗室，四川成都 610041；2.青藏高原動物遺傳資源保護與利用四川省重點實驗室，四川成都 610041)

全基因組測序(Whole Genome Sequencing，WGS)是對細胞或生物體所有的基因組進行測序，獲得完整的基因組信息.全基因組測序數(shù)據(jù)可以找到基因與表型之間的聯(lián)系，在挖掘動植物重要經(jīng)濟性狀相關的功能基因、分析遺傳機制等方面有重要意義.

我國牦牛主要分布在青藏高原海拔2 500米以上的高寒草原及高山峽谷區(qū)域，有著耐高寒和耐低氧的優(yōu)良性狀，是牧民生活依賴的重要生產(chǎn)資料[1-3].由于傳統(tǒng)牦牛飼養(yǎng)方式造成的長期過度繁殖和近親繁殖，家養(yǎng)牦牛的繁殖能力、生長速度、成年體型和產(chǎn)奶量下降，給當?shù)匦竽翗I(yè)的發(fā)展造成了影響[4]. 對全基因組測序數(shù)據(jù)進行分析，可以找到與生長性狀相關的候選基因，但是牦牛的全基因組測序數(shù)據(jù)一般有較多缺失值，需要利用基因型填充的方法填充缺失值.

基因型填充方法的基本原理是根據(jù)參考群體提供的基因型信息，構建出參考群體和目標群體之間共享的單倍型信息，比對目標群體與參考群體之間共享的單倍型信息，將目標群體缺失的基因型信息填充完整，得到完整的基因型數(shù)據(jù)[5]. 基因型填充方法大致分為兩類，一類是計算密集型的，如IMPUTE、MACH和fastPHASE 等，還有一類是計算高效型，如PLINK、MINIMAC 和BEAGLE 等. 計算密集型的基因型填充軟件在填充的過程中充分考慮到所有已知的基因型信息，使得對未知的基因型信息的估算更加精確，但是填充耗時較長.而計算高效型的基因型填充軟件在填充過程中僅僅關注與特定的SNP 位點相鄰的一部分標記的基因型信息，填充耗時較短但是填充準確性相對要低一點. 在日本黑牛的研究中，利用BEAGLE軟件從缺失基因型的數(shù)量、參考群體的大小、參考群體和目標群體之間的遺傳關系評估了日本黑牛測序數(shù)據(jù)的填充準確性.對于不同數(shù)量的缺失基因型的影響，50 K、26 K 和20 K 的測序數(shù)據(jù)的相關性都比較高，但7 K 的相關性較低；對于參考群體的大小的影響，50 K、26 K 和20 K 的測序數(shù)據(jù)在參考群體超過400 時的相關性趨于平穩(wěn)，但當參考群體數(shù)量增加到400 以上時，7 K 測序數(shù)據(jù)的相關性略有增加；對于參考群體和試驗群體之間遺傳關系的影響，使用相關參考群體的相關性高于使用無關參考群體的相關性[6].在肉牛的研究中，使用 BEAGLE、FIMPUTE 和 IMPUTE2 軟件對多品種的肉牛群體進行了填充. FIMPUTE 軟件對純種群體填充的填充匹配率在94.20%到 97.93% 之 間， IMPUTE2 軟 件 為 95.35% 到98.31%，BEAGLE 軟件為90.02%到96.38%.雜交動物的填充匹配率為54.15% ～97.53% (FIMPUTE)，57.04% ～ 97.46% (IMPUTE2)， 以 及 54.35% ～95.64%(BEAGLE)[7].在牛的研究中，選擇三種參考基因組，對數(shù)據(jù)進行基因分型后，利用PEDIMPUTE、FINCHAP、FIMPUTE 和 BEAGLE 軟件進行填充，F(xiàn)IMPUTE 軟件的填充匹配率最高，約為95%，其次是BEAGLE 軟件，約為92%，另外兩種軟件的填充匹配率受基因分型的影響變化較大[8].基因型填充技術在人類[9]、雞[10]、羊[11]、豬[12]、植物[13]上都有廣泛的應用.

本文將利用StochasticImpute 函數(shù)、impute. knn 算法和BEAGLE 軟件三種基因型填充方法對牦牛的測序數(shù)據(jù)進行填充，探究三種填充方法在不同的缺失率條件下的填充效果，以填充匹配率、相關性和填充耗時為評價指標，期望得到填充效果更好的填充方法，用于后續(xù)的研究中.

1 材料與方法

1.1 數(shù)據(jù)處理

本文所用的牦?；蛐蛿?shù)據(jù)來源于已發(fā)表的文章(JIA C，Animal Genetics，2020)[14]，文章測量了青海大通牧場354 頭雌性阿什旦牦牛的4 個生長性狀，并對所有個體的全基因組測序數(shù)據(jù)進行了基因分型和質量控制，得到了三個類型的文件，分別是ped 文件、map 文件和 csv 文件.其中 ped 文件包含了 354 頭牦牛的98 688 個SNP 位點的信息，map 文件包含了354 頭牦牛的遺傳圖譜信息，csv 文件包含了354 頭牦牛體重、肩高、體長和胸圍4 個生長性狀的信息.本研究將ped 文件轉換成hapmap 格式基因組信息文件.首先，用 PLINK 軟件將 ped 文件轉換為 0、1、2 型的文件，其中 0、1、2 型分別對應主要(major)等位基因、雜合基因和次要等位基因(minor). 因為本文是對不同的基因型填充方法進行比較，所以需要先用GATK 軟件對原始的98 k 數(shù)據(jù)進行過濾. 從過濾后的98 k 數(shù)據(jù)中隨機選取10 k 數(shù)據(jù)，總共生成30 組不同的10 k數(shù)據(jù)作為本文的真實數(shù)據(jù)集. 將每組數(shù)據(jù)按照5%、10%、15%和20%的比例進行人為的缺失，缺失的基因型用“NA”表示，把含有缺失基因型的數(shù)據(jù)作為參考數(shù)據(jù)集.缺失掉的基因型數(shù)據(jù)保留在另一份獨立文件中作為預測數(shù)據(jù)集的真實基因型.利用三種基因型填充方法對參考數(shù)據(jù)集中的缺失基因型進行預測，得到預測數(shù)據(jù)集.最后比較真實基因型數(shù)據(jù)和預測基因型數(shù)據(jù)的匹配率和相關性.

1.2 基因型填充方法

StochasticImpute 函數(shù)的填充原理是利用等位基因的頻率進行缺失值的填充，即以所有個體為參考群體，計算每個SNP 位點中所有等位基因的頻率，用頻率最高的等位基因填充該SNP 位點中的缺失值.本研究根據(jù)數(shù)值型基因型文件，利用StochasticImpute 函數(shù)分別計算0、1、2 三種基因型頻率，用頻率最高的等位基因型作為填充的基因型.

impute. knn 算法需要調(diào)用 R 語言中的“impute”軟件包，這種算法利用基因型數(shù)據(jù)中特定數(shù)目近鄰基因型值來填充含有缺失值的個體的基因型[15]. 首先需要將候選鄰居進行分類，再使用距離公式計算含有缺失值的基因與候選鄰居之間的距離，其中用來計算距離的基因坐標應為基因中未缺失的元素.對于候選鄰居可能缺少用于計算距離的坐標的情況，需要計算非缺失元素的平均值. 找到一個基因的k 個近鄰后，根據(jù)這k 個近鄰中的大部分鄰居所屬的類別決定含有缺失值的基因型，然后通過對其相鄰非缺失元素求平均值來估算缺失元素[16].本文設置的k 值為3，選擇了歐幾里德計算鄰近數(shù)據(jù)之間的距離.

本文選擇BEAGLE 5.1[17]作為填充軟件，該軟件利用隱馬爾可夫模型(Hidden Markov Models，HMM)進行基因型填充.首先利用HMM 計算參考面板的單倍型中一個標記到下一個標記的概率，其中每個標記指的是一段有相同等位基因的單倍型集合[18]. 在每個標記處，用等位基因標記的概率之和作為該等位基因的估算概率.從第一個標記到最后一個標記的概率之和就是特定的單倍型概率. 然后，對目標樣本進行基因分型，得到單倍型集合. 再根據(jù)目標樣本與參考面板之間共有的基因序列建立模型，利用計算得到的參考面板中的等位基因標記的概率，預測目標樣本中同樣標記處的缺失值，進行填充[17].本文需要將0、1、2 型文件進行轉換，將“0”轉換為“AA”，“1”轉換為“AB”，“2”轉換為“BB”，“NA”轉換為“??”，得到新的基因型文件，再利用BEAGLE 軟件進行填充.

1.3 評價基因型填充方法的標準

為了比較三種基因型填充方法的填充效果，本文將把填充匹配率、相關性和填充耗時作為評價標準.其中填充匹配率指的是填充正確的基因型個數(shù)與需要填充的基因型個數(shù)的比值. 其中Nmatch是填充后準確預測基因型的數(shù)目，Ntotal是所有缺失基因型的數(shù)目.

相關性計算為真實的基因型與填充得到的數(shù)值型基因型之間的相關性. 其中Gimpute是填充后的基因型，Greal是真實的基因型.

填充耗時指的是從填充開始一直到填充結束所用的時間.用R 語言中的system. time 函數(shù)來標記開始和結束時間，最后將軟件運行返回的時間作為填充耗時.

2 結果與分析

2.1 填充匹配率的比較

對于基因型數(shù)據(jù)，隨機設置5%、10%、15% 和20%的缺失率，用三種基因型填充方法進行填充，重復30 次，得到填充匹配率，如下圖1.

圖1 三種填充方法在不同缺失率條件下的匹配率Fig.1 Matching rate of three imputation methods under different missing rate

由圖1 可以看出，在不同的缺失率條件下，BEAGLE 軟件的填充匹配率最高，填充匹配率分別為0.863 0、0.861 7、0.860 9 和 0.859 9，其次是 StochasticImpute 函數(shù)，填充匹配率分別為0.812 5、0.812 3、0.812 5 和0.812 5，填充匹配率最低的是impute. knn算法，填充匹配率分別為0.626 5、0.621 4、0.614 8 和0.606 3. 三種方法的匹配率并不隨著缺失率的增大而產(chǎn)生的顯著的變化.

2.2 相關性的比較

利用30 次的重復抽樣數(shù)據(jù)，計算真實基因型數(shù)據(jù)數(shù)值與預測基因型數(shù)值之間的相關性，以30 次結果的平均值作為考量(圖2).

圖2 三種填充方法在不同缺失率條件下的相關性Fig.2 Correlation of three imputation methods under different missing rate conditions

由圖2 可以看出，在不同的缺失率條件下，StochasticImpute 函數(shù)和impute. knn 算法的相關性相對較高，StochasticImpute 函數(shù)的相關性分別為0.221 6、0.219 0、0.219 3、0.219 3，而 impute. knn 算法的相關性分別為 0.221 8、0.220 4、0.218 5、0.219 7. BEAGLE 軟件的相關性分別為 0.220 4、0.202 6、0.183 1、0.163 7，隨著缺失率的增大，BEAGLE 軟件的相關性在逐漸降低，并且?guī)缀醵嫉陀诹硗鈨煞N填充方法的相關性. 而隨著缺失率的增大，StochasticImpute 函數(shù)和impute. knn 算法的相關性基本沒有變化.

2.3 填充耗時的比較

在利用三種方法進行填充的時候，記錄三種填充方法重復30 次得到的填充耗時，取平均值，用對數(shù)函數(shù)進行標準化，得到圖3.

圖3 三種填充方法在不同缺失率條件下的填充耗時Fig.3 Impute time of three imputation methods under different deletion rates

由圖3 可以看出，BEAGLE 軟件的填充耗時最長，填充耗時分別為 380.7 s、465.3 s、531.0 s 和604.5 s，其次是 impute. knn 算法，填充耗時分別為13.0 s、14.0 s、14.4 s 和 15.4 s，填充耗時最少的是StochasticImpute 函數(shù)，填充耗時為7.2 s、7.7 s、7.9 s和8.0 s.并且隨著缺失率的增大，三種填充方法的填充耗時都在增加，其中BEAGLE 軟件的填充耗時增加得最多，在缺失率為20% 的條件下，填充耗時為604.5 s.

3 討論

因為相關性只考慮數(shù)值(即數(shù)值型基因型值)之間的順序關系，在基因型數(shù)據(jù)用0、1、2 的數(shù)值表示時會造成了巨大的偏差.在Korku 等[19]的研究報道中，與本文同樣使用的是0、1、2 型的基因型文件，但是填充得到的相關性和填充匹配率都比較高，而且填充匹配率比較高的群體，對應的相關性也比較高. 本文的填充匹配率比較高，但是相關性卻比較低，可能是因為數(shù)據(jù)中的0、1 型的基因型比較多.雖然BEAGLE 軟件填充正確的基因型比較多，但是這些基因型大多是0、1 型，也可能對相關性造成影響，導致相關性比較低.

因為StochasticImpute 函數(shù)的原理是根據(jù)等位基因的頻率填充缺失的基因型，每個缺失基因型填充的時候只利用所有參考群體該位置的基因型，所以填充耗時比較短.對于參考群體數(shù)量比較多，突變位點比較少的數(shù)據(jù)，用StochasticImpute 函數(shù)的填充效果比較好，填充準確率比較高，填充耗時也比較短.但是StochasticImpute 函數(shù)沒有考慮到群體內(nèi)部具體的遺傳關系和分層結構，而群體內(nèi)部具體的遺傳關系和分層結構在牦牛數(shù)據(jù)中尤其重要，因此StochasticImpute 函數(shù)不適用于填充牦牛數(shù)據(jù).

impute. knn 算法的填充效果與k 值有關，選擇合適的k 值，得到的填充效果比較好. 在Di 等[20]的研究中，因為選擇了合適的k 值，所以利用impute. knn算法填充得到的填充效果比較好.因為impute. knn 算法的k 值不能通過計算確定，需要人為設置，對于樣本量比較小的數(shù)據(jù)，可以計算所有的k 值找到最合適的值，而本文的樣本量比較大，沒有辦法計算所有的k 值，選擇k 值為3，相對于354 個個體而言很少，所以填充效果比較差.impute. knn 算法適合樣本量比較小的數(shù)據(jù).

在這三種填充方法中，BEAGLE 軟件的填充效果比較好，填充匹配率最高，填充耗時也較短. 在Yang等[21]的報道中，利用BEAGLE 軟件對水牛的測序數(shù)據(jù)進行填充，得到的填充匹配率為83%左右. 而本文利用BEAGLE 軟件對牦牛的測序數(shù)據(jù)進行填充得到的填充匹配率為86%左右，填充匹配率更高，說明可以利用BEAGLE 軟件對牦牛的基因型數(shù)據(jù)進行填充.在鄧天宇等[22]的研究中，利用BEAGLE 軟件填充模擬數(shù)據(jù)得到的填充耗時和本文的結果比較接近，都相對較短.并且隨著缺失率的增加，BEAGLE 軟件的填充耗時也在增加，但是增長的趨勢較慢.

本文研究了三種不同的基因型填充方法對牦牛的測序數(shù)據(jù)進行填充的效果，但是只對比了在不同缺失率條件下的結果，對于參考群體大小、目標群體位點比例、SNP 位點數(shù)目等方面并未進行研究，希望能夠在之后的研究中針對不同的方面展開更多的研究.

因為本文使用的測序數(shù)據(jù)是由354 頭牦牛，10 k測序數(shù)據(jù)組成，所以在基因型數(shù)據(jù)相對較小的時候，用BEAGLE 軟件進行填充的匹配率比較高，填充耗時也比較少.當需要填充的基因型數(shù)據(jù)比較多的時候，用BEAGLE 軟件進行填充時，填充匹配率會提高，但是填充耗時也會增長.

牦牛群體一般處于半野生的放牧環(huán)境下，因此家系和亞群的背景在遺傳關系中具有非常強的結構關系.利用這些信息來進一步提升牦?；蛐蛿?shù)據(jù)填充的準確性，可能會是一個理想的研究方向.

4 結論

通過對三種填充方法的填充效果進行對比，BEAGLE 軟件的填充匹配率更高，impute. knn 算法的相關性更高，StochasticImpute 函數(shù)的填充耗時更短.在缺失率小于20%時，基因型填充方法的填充效果較好，建議使用20%的缺失率對測序數(shù)據(jù)進行過濾，過濾后的數(shù)據(jù)利用基因型填充技術填充可以得到可信度較高的基因型數(shù)據(jù).