約束標(biāo)準(zhǔn)化線性回歸法估計(jì)合成品種動物基因組品種構(gòu)成

何俊, 李智,2, 吳曉林,2

約束標(biāo)準(zhǔn)化線性回歸法估計(jì)合成品種動物基因組品種構(gòu)成

何俊1, 李智1,2, 吳曉林1,2

（1湖南農(nóng)業(yè)大學(xué)動物科技學(xué)院，中國長沙 410128；2美國紐勤公司生物信息與生物統(tǒng)計(jì)部，美國林肯市 68504）

合成品種是由至少兩種純種（祖先）培育的新品種，旨在兼顧祖先品種的有利遺傳特征，并且可以長期保持后代的雜種優(yōu)勢而不需要每個世代都雜交。合成品種的遺傳穩(wěn)定，不同于雜交群體，因而可以像純種一樣繁育。實(shí)踐中，估計(jì)合成品種的祖先品種對每個動物個體基因組的遺傳貢獻(xiàn)比例，即基因組品種構(gòu)成（genomic breeding composition, GBC），在畜禽品種登記、品種培育歷史和品種構(gòu)成分析、品種保護(hù)和雜交優(yōu)勢預(yù)測等方面有著非常重要的意義。利用基因組SNP基因型數(shù)據(jù)，采用合適的數(shù)學(xué)模型和統(tǒng)計(jì)方法，可以鑒定現(xiàn)有純種品種的動物個體或純種品種在雜交個體基因組的遺傳貢獻(xiàn)比例，而估計(jì)合成品種GBC的方法和研究都較少?！尽烤€性回歸是估計(jì)GBC的常用方法之一，但也存在諸多的問題。本研究旨在提出和評估一種約束的標(biāo)準(zhǔn)化線性回歸方法（restricted standardized linear regression, RSLR），作為傳統(tǒng)線性回歸方法的改進(jìn)方法，應(yīng)用于估計(jì)合成品種動物個體的GBC。采用肉牛王牛（Beefmaster）及其3個祖先品種（婆羅門牛、海福特牛和短角牛）的GGP 50K SNP芯片所測定的基因型數(shù)據(jù)，通過計(jì)算其基因頻率和歐氏距離，利用層次聚類分析方法解析了4個動物群體的遺傳關(guān)系，然后提出了RSLR方法，估計(jì)合成品種動物個體GBC的原理和方法。為了檢驗(yàn)該方法的估計(jì)效果，從基因型數(shù)據(jù)中選擇了均勻分布的分別包含1 000、5 000、10 000、20 000、30 000、40 000個SNP以及3個祖先品種共有的47 900個SNP的7個子集，分別采用RSLR和傳統(tǒng)線性回歸（linear regression, LR）兩種方法估計(jì)了4 323頭肉牛王牛的GBC，并比較了兩種方法的計(jì)算結(jié)果。聚類分析的結(jié)果與4個品種間的遺傳關(guān)系相吻合，表明肉牛王牛與婆羅門牛的遺傳關(guān)系最近，遺傳距離小于其與海福特牛和短角牛的遺傳距離。LR方法估計(jì)的GBC會低估婆羅門牛（0.459—0.462）和短角牛（0.208—0.212）對于肉牛王牛的基因組貢獻(xiàn)，同時高估海福特牛（0.326—0.333）的基因組貢獻(xiàn)。但RSLR方法估計(jì)的肉牛王牛GBC的平均值與3個祖先品種預(yù)期的基因組貢獻(xiàn)比例比較吻合：婆羅門牛為0.497—0.503，海福特牛為0.262—0.274，短角牛為0.229—0.231。此外，LR方法估計(jì)GBC的標(biāo)準(zhǔn)差和變異系數(shù)明顯大于用RSLR估計(jì)的結(jié)果。當(dāng)SNP子集數(shù)量在20 000以上時，LR方法估計(jì)牛肉王牛的3個祖先品種婆羅門牛、海福特牛和短角牛基因組貢獻(xiàn)的標(biāo)準(zhǔn)差分別為0.048、0.032和0.051—0.052，變異系數(shù)分別為10.46%—10.50%、9.61%—9.76%和23.94%—25.00%，而RSLR方法估計(jì)的標(biāo)準(zhǔn)差，3個祖先品種對應(yīng)為0.021、0.021—0.022和0.024—0.025，變異系數(shù)分別為4.18%—4.20%、7.89%—8.33%以及10.26%—10.68%。用RSLR方法估計(jì)的合成品種肉牛王牛動物個體的GBC，比LR方法的估計(jì)結(jié)果更加準(zhǔn)確，估計(jì)的結(jié)果比LR方法估計(jì)的結(jié)果更穩(wěn)定，且估計(jì)的一致性也更好，可以作為線性回歸方法的改進(jìn)，應(yīng)用于估計(jì)合成品種動物個體GBC。

SNP芯片；線性回歸；合成品種；基因組品種構(gòu)成

0 引言

【研究背景】合成品種是綜合了兩個或更多純種品種的性狀特征而培育的新品種，例如肉牛王牛、布蘭格斯牛等。合成品種不同于一般的簡單雜交群體，合成品種遺傳穩(wěn)定，可以像純種品種一樣進(jìn)行本品種內(nèi)繁育（包括一定程度的近交繁育）。事實(shí)上，現(xiàn)在許多的純種品種，如果追溯到足夠久遠(yuǎn)的年代，也都是合成品種。合成品種兼顧了其祖先品種性狀的優(yōu)勢，同時又可以避免這些品種的一些劣勢，而不需要繼續(xù)雜交繁育，因此具有較高的經(jīng)濟(jì)價值。通過合成雜交繁育而培育的肉牛、奶牛、綿羊、豬和家禽，已經(jīng)成為動物商品生產(chǎn)的一個重要方式，為畜禽商品生產(chǎn)提供了優(yōu)質(zhì)的種源[1]。合成品種動物個體的基因組品種構(gòu)成（genomic breed composition, GBC），是指祖先品種對于該個體的基因組遺傳率，也可以簡單理解為合成品種動物個體的基因組與祖先品種基因組的相似性的百分率。例如，布蘭格斯牛是安格斯牛和婆羅門牛的雜交后裔。從群體平均而言，布蘭格斯牛在遺傳上有3/8的婆羅門牛和5/8的安格斯牛血統(tǒng)[2]。肉牛王牛是20世紀(jì)30年代用海福特母牛和短角牛母牛與婆羅門公牛雜交培育而成的肉牛品種，平均含有25%的海福特牛、25%的短角牛以及50%的婆羅門牛血統(tǒng)[2]?！狙芯恳饬x】在動物遺傳育種中，估計(jì)動物個體的GBC具有廣泛的應(yīng)用價值，如了解和評估某動物品種的育成歷史和品種純度、地方品種保種、雜種優(yōu)勢預(yù)測以及設(shè)計(jì)雜交計(jì)劃和制定雜交育種方案等[3-4]。【前人研究進(jìn)展】動物個體的GBC通?？梢杂孟底V資料或基因組分子標(biāo)記來估計(jì)。從理論上講，后者比前者估計(jì)的GBC更準(zhǔn)確，因?yàn)橛没蚪M分子標(biāo)記估計(jì)GBC，不僅不受系譜錯誤的影響，還可以反映出實(shí)際的遺傳抽樣導(dǎo)致GBC的偏差。并且，用系譜計(jì)算的是平均（期望）GBC，沒有反映出孟德爾抽樣所導(dǎo)致的親本遺傳貢獻(xiàn)率上的偏差[5]，而用基因組標(biāo)記估計(jì)的是真實(shí)GBC。因而利用一套全基因組的SNP基因分型數(shù)據(jù)，采用合適的數(shù)學(xué)模型和統(tǒng)計(jì)方法，可以鑒定現(xiàn)有純種品種的動物個體，或是估計(jì)純種品種在雜交個體基因組的遺傳貢獻(xiàn)比例[6-9]。采用SNP標(biāo)記估計(jì)動物個體GBC的統(tǒng)計(jì)方法很多[10-13]，例如主成分分析方法[14-16]、混合分布方法[3,8,17-18]、線性回歸分析方法[10-11]。此外，Dodds等[19]將基因組預(yù)測模型（基因組BLUP）方法應(yīng)用于動物個體基因組品種構(gòu)成的估計(jì)。在這些方法中，線性回歸模型的方法比較簡便。該方法以參考（祖先）品種的等位基因頻率作為自變量，待測動物的基因型為依變量，計(jì)算參考群體的基因頻率對于每個個體等位基因計(jì)數(shù)的回歸系數(shù)。該方法目前已用于估計(jì)豬和牛的品種遺傳構(gòu)成[3,11,20-21]?！颈狙芯壳腥朦c(diǎn)】線性回歸方法是估計(jì)GBC的最常用方法之一。但是，線性回歸模型估計(jì)的GBC實(shí)際上是各祖先的基因頻率對于個體動物基因型的回歸系數(shù)。對于動物個體而言，其多個祖先品種的回歸系數(shù)之和并不一定等于1，因?yàn)榛貧w系數(shù)是有理實(shí)數(shù)，其數(shù)值可以超過1，也可以為負(fù)數(shù)。因此，用傳統(tǒng)線性回歸模型估計(jì)的GBC需要校正，使其和為1[3,6,8]?！緮M解決的關(guān)鍵問題】本研究一是利用約束條件下的標(biāo)準(zhǔn)化變量的線性轉(zhuǎn)換，提出了一個改進(jìn)的線性回歸方法來估計(jì)GBC，稱為約束的標(biāo)準(zhǔn)化變量線性回歸分析（restricted standardized linear regression，RSLR）方法。該方法不需要對所估計(jì)的祖先品種的回歸系數(shù)近似校正，而是直接估計(jì)出動物個體的GBC；二是以合成品種肉牛王牛為例，比較了RSLR方法和傳統(tǒng)的線性回歸方法（linear regression，LR）估計(jì)GBC的實(shí)際效果，為估計(jì)合成品種動物個體的GBC提供更為合適的估計(jì)方法。

1 材料與方法

1.1 試驗(yàn)材料

收集了4 323頭肉牛王牛，68頭婆羅門牛，1 232頭短角牛和2 423頭海福特牛的GGP 50K SNP基因型數(shù)據(jù)?；蛐蛿?shù)據(jù)由美國紐勤GeneSeek公司提供，每個個體有49 463個SNP位點(diǎn)的基因型。缺失基因型通過FImpute軟件來填充[22]。在基因型數(shù)據(jù)中，刪掉Y染色體和線粒體上的SNP基因型，保留芯片共有的47 900個SNP用于后續(xù)分析。首先計(jì)算了4個品種所有SNP的等位基因頻率，利用品種間基因頻率的歐氏距離[3,8,23]進(jìn)行Ward.D2層次聚類分析[24-26]，解析了4個品種的群體結(jié)構(gòu)；所有計(jì)算和分析過程均采用R及自編的R程序包完成。4個群體的動物數(shù)量及GGP 50K SNP的小等位基因頻率（minor allele frequency，MAF）的群體均值和標(biāo)準(zhǔn)差列于表1。

表1 4個牛群體的動物個體數(shù)目以及GGP 50K SNP芯片的基本信息

1.2 約束的標(biāo)準(zhǔn)化變量線性回歸分析方法

通過約束條件下標(biāo)準(zhǔn)化變量的線性轉(zhuǎn)換，改進(jìn)了傳統(tǒng)的線性回歸模型估計(jì)GBC的方法。該方法提供的標(biāo)準(zhǔn)化線性回歸系數(shù)可以作為其基因組品種構(gòu)成的直接估計(jì)，因而不再需要對線性回歸系數(shù)進(jìn)行校正。該方法的具體過程介紹如下。

設(shè)為一個個體所有個SNP的基因型向量（×1），其中SNP基因型分別用0 (AA)、1 (AB)、2 (BB) 表示。設(shè)=（1…T）是一個×的參考群體（祖先品種）基因頻率的向量，其中f為一個×1的向量，包含第個參考群體中所有SNP座位的等位基因（例如B等位基因）的頻率=（1,…,）。因此，GBC可以采用下列的線性回歸模型估計(jì)：

式中，是總體均數(shù)，=(1…T)′ 是×1的品種回歸系數(shù)向量，是誤差向量。以上即為LR模型估計(jì)GBC的方法。理論上，每個動物個體的GBC之和應(yīng)該為1。但在LR模型中，對每一個動物個體而言，個品種的回歸系數(shù)之和并不一定等于1，因此，用LR方法估計(jì)GBC需要對回歸系數(shù)近似校正，使其和為1[3,6,8]。

如果對上述線性回歸模型中的線性變量先標(biāo)準(zhǔn)化，然后約束標(biāo)準(zhǔn)化的（祖先品種）線性回歸系數(shù)為1，就可以避免對回歸系數(shù)進(jìn)行近似校正。因此標(biāo)準(zhǔn)化的（祖先品種）線性回歸系數(shù)可以直接作為GBC的估計(jì)值。

首先計(jì)算式（1）的平均值：

因?yàn)?)=，()=0。然后將線性變量標(biāo)準(zhǔn)化，即將式（1）兩邊同時減去式（2）兩邊的相應(yīng)的平均值，然后再除以g的標(biāo)準(zhǔn)差（σ），這樣可得到：

式中，σ為第個參考群體（祖先品種）的M個SNP的等位基因(例如B)頻率的標(biāo)準(zhǔn)差。令，，，，則式（3）可進(jìn)一步簡化為：

式中，為標(biāo)準(zhǔn)化的基因型向量（×1），x為第個參考群體（祖先品種）的標(biāo)準(zhǔn)化等位基因頻率向量（×1），p為標(biāo)準(zhǔn)化的回歸系數(shù)，即通徑系數(shù)[27]。

從祖先品種對動物個體基因組遺傳貢獻(xiàn)的角度看，每個個體的GBC總和應(yīng)該等于1。因此在的約束條件下做回歸變量的線性轉(zhuǎn)換。令，則式（4）可變?yōu)椋?/p>

又令=-x，z= x-x，c= p，因此式（5）可以表示如下：

式中，c為第個參考品種（群體）的GBC，= 1, …, T-1。最后一個參考群體（祖先品種）的GBC（c）可通過c=1-（1+…+c-1）進(jìn)行計(jì)算。

1.3 SNP子集的選擇

使用了7個SNP子集估計(jì)肉牛王牛動物個體的GBC，其中6個為從GGP 50K SNP中選擇的均勻分布的SNP子集，SNP的數(shù)目分別為1 000、5 000、10 000、20 000、30 000和40 000。還有一個SNP子集為包括了數(shù)據(jù)清理后的全部共有的47 900 SNP。

2 結(jié)果

2.1 4個品種的群體聚類分析

為了了解四個品種的群體結(jié)構(gòu)和遺傳背景，對4個群體的聚類分析表明（圖1），肉牛王牛和婆羅門牛先聚成一類，然后和聚成一類的海福特牛和短角牛再聚在一起，這與肉牛王牛的3個祖先品種的血緣構(gòu)成比例是相符合的，婆羅門牛占血緣構(gòu)成的50%，所以和肉牛王牛遺傳距離最近，其他兩個祖先品種各占25%，相對于肉牛王牛而言，它們距離相似，因而聚成一類。

2.2 7個子集的SNP在染色體上的分布

選擇了6個均勻分布的SNP子集以及數(shù)據(jù)清理后的全部47 900 SNP，每個SNP子集中的SNP數(shù)目從1 000到47 900不等，每個子集中的SNP在每條染色體上的分布數(shù)量見表2。

表2 選擇的7個子集中的SNP數(shù)量在染色體上的分布

0號染色體表示該SNP所在的染色體信息未知 Chromosome 0 indicates that the information of the chromosome where the SNP is located is unknown

2.3 肉牛王牛祖先品種的GBC估計(jì)

分別用LR和RSLR方法，估計(jì)了4 323頭肉牛王牛的GBC（表2）。用LR方法估計(jì)的3個祖先品種對于肉牛王牛的GBC分別為：0.457—0.463（婆羅門牛），0.322—0.338（海福特牛）以及0.208—0.216（短角牛）標(biāo)準(zhǔn)差依次分別為0.048—0.060、0.032—0.054、0.051—0.073。采用RSLR方法估計(jì)3個祖先品種對于肉牛王牛的GBC分別為：0.497—0.503（婆羅門牛）、0.262—0.274（海福特牛）和0.229—0.235（短角牛），3個品種的標(biāo)準(zhǔn)差依次分別為：0.021—0.029、0.021— 0.036、0.024—0.038?？梢?，用RSLR方法估計(jì)的肉牛王牛的GBC比LR方法所估計(jì)的GBC更加接近于所期望的群體均值。從所估計(jì)的GBC中位數(shù)看也是如此（表3）。相比之下，LR方法估計(jì)的GBC與期望的GBC偏差較大，特別是低估了肉牛王牛與婆羅門牛的基因組相似性，同時高估了肉牛王牛與海福特牛的基因組相似性。

圖1 四個品種的群體結(jié)構(gòu)分析

表3 兩種回歸分析方法和7個SNP集分別估計(jì)的肉牛王牛祖先品種的GBC

比較了LR和RSLR兩個方法用7個不同SNP子集計(jì)算的GBC的變異系數(shù)（表3）?？梢钥闯觯旱谝?，LR方法估計(jì)GBC的變異系數(shù)（10.46%—34.43%）明顯大于用RSLR方法計(jì)算的GBC變異系數(shù)（4.18%—16.59%），表明用RSLR方法估計(jì)GBC的個體間差異要遠(yuǎn)小于LR方法。第二，兩個方法估計(jì)的GBC的變異系數(shù)都隨著子集SNP數(shù)增加而降低，但是，RSLR估計(jì)的GBC的變化趨勢也要遠(yuǎn)小于LR估計(jì)的GBC的變化趨勢。例如，當(dāng)SNP數(shù)由1 000逐步增加到47 900時，用LR估計(jì)3個祖先品種的遺傳貢獻(xiàn)比例分別由12.99%降到10.46%（婆羅門牛），16.56%降到9.61%（海福特牛），34.43%降到25.00%（短角牛）。與此相比，RSLR方法在7個SNP子集中，除了1 000 SNP時估計(jì)的變異系數(shù)稍高，隨著SNP數(shù)增加，3個祖先品種的變異系數(shù)都比較小，而且取值范圍都比較接近，分別為4.19%—4.59%（婆羅門牛），7.89%—9.36%（海福特牛）和10.26%—11.69%（短角牛）。兩個方法都表明隨著SNP數(shù)的增加，GBC估值在個體間的變異呈現(xiàn)降低的趨勢?？傮w而言，用回歸模型的方法，GBC估值的變異系數(shù)在5 000 SNP以上基本都趨于穩(wěn)定。

作為初步的研究結(jié)果，本研究參考群體（祖先品種）中婆羅門牛的樣本數(shù)目偏少，因此有必要將來用更大的參考群體樣本進(jìn)行驗(yàn)證。從本研究的結(jié)果看，所估計(jì)的GBC與預(yù)期的GBC基本吻合，表明估計(jì)的基因型頻率大體上是比較準(zhǔn)確的。小樣本數(shù)據(jù)中主要對MAF很低的SNP的基因頻率估計(jì)偏差比較大（如稀有小等位基因頻率位點(diǎn)），但這些SNP等位基因頻率的偏差，對估計(jì)GBC的影響非常有限。

從動物個體看，估計(jì)肉牛王牛個體的3個祖先品種的GBC有一定的變化幅度，這是由于在品種繁育過程中實(shí)際遺傳抽樣的結(jié)果。如以RSLR方法用全部47 900 SNP估計(jì)的結(jié)果看，GBC的范圍為：[0.401，0.575]（婆羅門牛）、[0.116，0.338]（海福特牛）、[0.167，0.393]（短角牛）；GBC的95%的置信區(qū)間為：[0.454，0.541]（婆羅門牛）、[0.223，0.308]（海福特牛）、[0.197，0.302]（短角牛）。RSLR方法用全部47 900 SNP估計(jì)肉牛王牛個體的3個祖先品種GBC的分布見圖2。

圖2 采用RSLR方法用全部47900 SNP估計(jì)3個祖先品種GBC的分布

3 討論

3.1 用線性回歸的方法估計(jì)GBC

用線性回歸方法估計(jì)動物GBC，方法簡單實(shí)用，是一個非常值得推廣的方法。但傳統(tǒng)的LR方法估計(jì)的GBC實(shí)際上是動物個體基因型對于參考群體（祖先品種）相應(yīng)等位基因頻率的回歸系數(shù)，就數(shù)值而言，回歸系數(shù)可以取任何一個實(shí)數(shù)數(shù)值。因此每個個體的所有祖先品種的回歸系數(shù)之和不一定等于1。VanRaden等[6]提出一個校正品種回歸系數(shù)的方法，用校正后的回歸系數(shù)的相對值作為GBC的估計(jì)，但是該校正方法在計(jì)算上比較繁瑣。作者等曾提出了一個簡化方法，即將所有負(fù)回歸系數(shù)設(shè)為零，然后計(jì)算每個個體的參考群體回歸系數(shù)的相對值作為GBC的估計(jì)值[3,8]。這兩個方法在結(jié)果上接近，然而這些校正方法是經(jīng)驗(yàn)式的，沒有任何的理論依據(jù)。

本研究采用標(biāo)準(zhǔn)化線性變量的約束條件作為LR的改進(jìn)方法，約束條件是標(biāo)準(zhǔn)化的回歸系數(shù)之和為1。這樣就可以避免對于傳統(tǒng)回歸系數(shù)的校正。當(dāng)祖先品種間完全沒有遺傳親緣關(guān)系的時候，這個約束條件是合理的，否則就是近似的。標(biāo)準(zhǔn)化的回歸系數(shù)，即通徑系數(shù)。從通徑分析的理論看，決定兩個變量（個體或群體）間相似性（相關(guān)系數(shù)）的因素包括它們二者之間的直接通徑關(guān)系和通過第三個變量（個體或群體）的間接通徑關(guān)系。當(dāng)間接通徑關(guān)系忽略不計(jì)的時候，兩個變量（個體或群體）間的相關(guān)系數(shù)，就等于二者之間的直接通徑系數(shù)[27]。因此可以合理假設(shè)，如果祖先品種間沒有遺傳親緣關(guān)系，用改進(jìn)線性回歸模型估計(jì)的祖先品種的標(biāo)準(zhǔn)化回歸系數(shù)（通徑系數(shù)）可以作為每個祖先品種和合成品種動物個體基因組貢獻(xiàn)率（或基因組相似程度）的估計(jì)。從品種馴化的歷史過程看，每個畜禽品種在起源上都可能是相關(guān)聯(lián)的，但在祖先品種間的遺傳親緣關(guān)系比較久遠(yuǎn)的情況下，這個假設(shè)是近似成立的。此外，需要說明的是，本研究中約束條件是標(biāo)準(zhǔn)化的回歸系數(shù)（通徑系數(shù)）之和為1，這不完全等同于通徑分析。就后者而言，所有因素直接通徑的決定系數(shù)和間接通徑的決定系數(shù)之和為1，因此，從通徑分析的角度，RSLR仍然是一個近似的方法。

3.2 SNP的選擇

SNP的選擇對GBC的估計(jì)結(jié)果有一定影響。并且不同的方法對于SNP選擇的要求也不盡相同。例如，混合模型方法要求選擇信息量高的SNP，這包括群體特有或是群體間差異大的SNP。Hulsegge等[12]比較了3個統(tǒng)計(jì)指標(biāo)用以衡量標(biāo)記信息量的效果，這3個統(tǒng)計(jì)指標(biāo)分別是Delta、Wright的FST以及Weir和Cockerham的FST。筆者通過最大化SNP基因頻率的平均歐式距離來篩選SNPs[8]。除此而外，信息熵[28-29]和主成分分析[15-16]中的加載系數(shù)[30]也是衡量SNP信息量的指標(biāo)[31]。但值得說明的是，回歸模型中選擇變量（SNP）可能導(dǎo)致選擇偏性，特別是對于線性回歸的方法。因篇幅所限，本文沒有詳細(xì)討論這個問題。本研究中沒有選擇信息含量高的SNP，而是選擇均勻分布的SNP。另一方面，線性回歸模型一般都需要比較多的SNP數(shù)目。在此情形下，使用均勻分布的SNP，可以較好的覆蓋整個基因組，使結(jié)果更具有代表性[8]。

降低SNP之間的連鎖不平衡也是一個需要考慮的因素。特別是對于混合分布模型，其似然函數(shù)的假設(shè)前提是SNP之間沒有關(guān)聯(lián)。尤其用高密度SNP估計(jì)GBC，需要盡量減少或刪除處于高度連鎖不平衡的SNP。Hulsegge等[12]采用LD的2＞0.30作為刪除SNP的標(biāo)準(zhǔn)，結(jié)果表明在保持相同準(zhǔn)確性的前提下，使用這種方法篩選SNP，可以明顯降低所需SNP標(biāo)記的數(shù)目。SNP間LD的程度對于線性回歸模型而言，沒有混合分布模型那樣重要。本研究沒有選擇信息含量高的SNP，也沒有作降低SNP之間LD的處理，而是選擇均勻分布的SNP。結(jié)果表明，對于中、低密度的SNP（50K SNP以內(nèi)），在不考慮SNP間LD的情形下，所估計(jì)的肉牛王牛的GBC與期望的群體均值也是基本上吻合的。此外，值得一提的是，本研究中當(dāng)SNP子集為5 000以上時，估計(jì)的結(jié)果已趨于穩(wěn)定，在20 000以上時結(jié)果已經(jīng)穩(wěn)定，說明在不增加實(shí)驗(yàn)室檢測成本的情況下，利用現(xiàn)有SNP芯片數(shù)據(jù)篩選可應(yīng)用于GBC估計(jì)的SNP子集是完全可行的，因而當(dāng)前使用的中低密度芯片數(shù)據(jù)完全可以滿足品種GBC的分析，這是對現(xiàn)有SNP芯片功能的深入開發(fā)與拓展，也是對芯片數(shù)據(jù)的分析和應(yīng)用的進(jìn)一步挖掘。

3.3 肉牛王牛的基因組品種構(gòu)成

肉牛王牛于1954年首次被美國農(nóng)業(yè)部認(rèn)定為新品種。最初的目的是培育出能夠適應(yīng)德克薩斯州南部環(huán)境的一個牛品種。目前的肉牛王牛是一個多用途品種，可用于牛奶和牛肉生產(chǎn)。根據(jù)官方數(shù)據(jù)，肉牛王牛平均包含50%的婆羅門牛、25%的海福特牛以及25%的短角牛的血統(tǒng)。本研究中，RSLR方法估計(jì)4 323頭肉牛王牛GBC的結(jié)果，估計(jì)的3個祖先品種的GBC的群體均值分別為：0.501（婆羅門牛）、0.265（海福特牛）和0.234（短角牛），基本與官方數(shù)據(jù)相符。肉牛王牛與海福特牛的基因組相似性稍高于25%，而與短角牛則稍低于25%，這個差異可能是由于該品種合成過程中因?yàn)檫x擇而產(chǎn)生的偏差。當(dāng)然，統(tǒng)計(jì)方法在估計(jì)上的偏差也不能完全排除。對于肉牛王牛的3個祖先品種而言，婆羅門牛是從印度進(jìn)口的牛品種中繁殖而來的，該品種與海福特牛和短角牛在遺傳關(guān)系上比較遠(yuǎn)。相比之下，海福特牛和短角牛都屬于原產(chǎn)于英國的牛品種，它們之間可能存在一定的遺傳關(guān)系。這可能也是導(dǎo)致肉牛王牛與海福特牛和短角牛的基因組相似性產(chǎn)生偏離的原因之一。

用基因組標(biāo)記估計(jì)動物個體GBC，可以反應(yīng)出個體水平上的遺傳抽樣，是實(shí)現(xiàn)了的個體基因組品種構(gòu)成的估計(jì)值。因此所估計(jì)的動物個體GBC在群體中存在一定的變異。本研究用RSLR方法估計(jì)肉牛王牛3個祖先品種的基因組貢獻(xiàn)率。實(shí)踐中，GBC的95%的置信區(qū)間可以作為肉牛王牛品種登記的分子標(biāo)記依據(jù)，從而可避免由于系譜資料缺失或誤差所導(dǎo)致的錯誤。

4 結(jié)論

本研究利用基因組SNP數(shù)據(jù)，對傳統(tǒng)的LR方法進(jìn)行了改進(jìn)，提出了RSLR的估計(jì)方法估計(jì)動物個體GBC。在對合成品種肉牛王牛個體的GBC估計(jì)中，與LR方法比較，RSLR方法的估計(jì)結(jié)果的準(zhǔn)確度和一致性更好，可將RSLR方法作為一種估計(jì)合成品種GBC的合適方法。若對方法做進(jìn)一步改進(jìn)，將需考慮親本品種間的遺傳相關(guān)，采用完全的通徑分析方法來估計(jì)GBC。

[1] 劉文忠.家畜合成群體保留雜種優(yōu)勢的預(yù)測與培育效果評價. 遺傳, 2009, 31(8):791-798.

Liu W Z. Prediction of retained heterosis and evaluation on breeding effects of composite livestock populations., 2009, 31(8):791-798.(in Chinese)

[2] Marshall B H, Briggs D M.. 4th ed. New York: MacMillian Company, 1980.

[3] 何俊, 錢長嵩, Richard G Tait Jr, Stewart Bauck, 吳曉林. SNP芯片數(shù)據(jù)估計(jì)動物個體基因組品種構(gòu)成的方法及應(yīng)用. 遺傳, 2018, 40(4):305-314.

He J, Qian C S, Tait Jr R G, Bauck S, Wu X L. Estimating genomic breed composition of individual animals using selected SNPs., 2018, 40(4):305-314. (in Chinese)

[4] Wu X L, Liu R Z, Shi Q S, Liu X C, Li X, Wu M S. Marker-assisted mating applied in in-situ conservation of indigenous animals in small populations: (1) Choosing mating schemes for maximum heterozygosity., 2000, 13(4): 431-434.

[5] 楊子博, 王安邦, 冷蘇鳳, 顧正中, 周羊梅. 小麥新品種淮麥33的遺傳構(gòu)成分析. 中國農(nóng)業(yè)科學(xué), 2018, 51 (17):3237-3248.

YANG Z B, WANG A B, LENG S F, GU Z Z, ZHOU Y M. Genetic analysis of the novel high-yielding wheat cultivar Huaimai33., 2018, 51(17): 3237-3248. ( in Chinese)

[6] VanRaden P M, Cooper T A. Genomic evaluations and breed composition for crossbred U.S. dairy cattle.Orlando, Florida, 2015.

[7] Pritchard J K, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data., 2000, 155(2): 945-959.

[8] He J, Guo Y G, Xu J, Li H, Fuller A, Tait R G, Wu X L, Bauck S. Comparing SNP panels and statistical methods for estimating genomic breed composition of individual animals in ten cattle breeds.2018, 19: 56.

[9] Gobena M, Elzo M A, Mateescu R G. Population structure and genomic breed composition in an Angus-Brahman crossbred cattle population., 2018, 9: 90.

[10] Chiang C W K, Gajdos Z K Z, Korn J M, Kuruvilla F G, Butler J L, Hackett R, Guiducci C, Nguyen T T, Wilks R, Forrester T, Haiman C A, Henderson K D, Le Marchand L, Henderson B E, Palmert M R, McKenzie C A, Lyon H N, Cooper R S, Zhu X F, Hirschhorn J N. Rapid assessment of genetic ancestry in populations of unknown origin by genome-wide genotyping of pooled samples., 2010, 6(3): e1000866.

[11] Kuehn L A, Keele J W, Bennett G L, McDaneld T G, Smith T P L, Snelling W M, Sonstegard T S, Thallman R M. Predicting breed composition using breed frequencies of 50,000 markers from the US Meat Animal Research Center 2,000 Bull Project., 2011, 89(6): 1742-1750.

[12] Hulsegge B, Calus M P, Windig J J, Hoving-Bolink A H, Maurice-van Eijndhoven M H, Hiemstra S J. Selection of SNP from 50K and 777K arrays to predict breed of origin in cattle, 2013, 91:5128-5134.

[13] AKANNO E C, CHEN L, ABO-ISMAIL M K, CROWLEY J J, WANG Z, LI C, BASARAB J A, MACNEIL M D, PLASTOW G. Genomic prediction of breed composition and heterosis effects in Angus, Charolais, and Hereford crosses using 50K genotypes., 2017, 97(3):431-438.

[14] McVean G. A Genealogical Interpretation Of Principal Components Analysis.. 2009, 5(10): e1000686.

[15] Ma J, Amos C I. Principal components analysis of population admixture.,2012,7(7): e40115.

[16] LEWIS J, ABAS Z, DADOUSIS C, LYKIDIS D, PASCHOU P, DRINEAS P. Tracing cattle breeds with principal components analysis ancestry informative SNPs.. 2011, 6(4):e18007.

[17] Bansal V, Libiger O. Fast individual ancestry inference from DNA sequence data leveraging allele frequencies for multiple populations., 2015, 16: 4.

[18] ALEXANDER D H, LANGE K. Enhancements to the ADMIXTURE algorithm for individual ancestry estima-tion., 2011, 12: 246.

[19] Dodds K G, Auvray B, Newman N S A, McEwan C J. Genomic breed prediction in New Zealand sheep., 2014, 15:92

[20] Funkhouser S A, Bates R O, Ernst C W, Newcom D, Steibel J P. Estimation of genome-wide and locus-specific breed composition in pigs., 2017, 1(1):36-44.

[21] GOBENA M, ELZO M A, MATEESCU R G. Population structure and genomic breed composition in an Angus-Brahman crossbred cattle population.2018, 9:90.

[22] SARGOLZAEI M, CHESNAIS J P, SCHENKEL F S. A new approach for efficient genotype imputation using information from relatives., 2014,15:478.

[23] HE J, GUO YG, XU JQ, LI H, FULLER A, RICHARD G JR, WU XL, BAUCK S. Estimating genomic breed composition of individual animals in ten cattle breeds: Comparison of SNP panels and statistical methodology//. New Zealand: Auckland, 2018, 684- 687.

[24] MURTAGH F, LEGENDRE P. Ward's hierarchical agglomerative clustering method: which algorithms implement Ward's criterion?2014, 31(3): 274-295.

[25] LEGENDRE P, LEGENDRE L.. 3rd ed. Developments in environmental modelling. 2012, 24.

[26] 桑世飛, 王會, 梅德圣, 劉佳, 付麗, 王軍, 汪文祥, 胡瓊. 利用全基因組SNP芯片分析油菜遺傳距離與雜種優(yōu)勢的關(guān)系. 中國農(nóng)業(yè)科學(xué), 2015, 48(12): 2469-2478.

SANG S F, WANG H, MEI D S, LIU J, FU L, WANG J, WANG W X, HU Q. Correlation analysis between heterosis and genetic distance evaluated by genome-wide SNP chip in., 2015, 48(12): 2469-2478.

[27] Wright S. Correlation and causation., 1921, 20(7): 557-585.

[28] HAN T S, KOBAYASHI K.. Boston, MA, USA: American Mathematical Society, 2001.

[29] WU X L, XU J Q, FENG G F, WIGGANS G R, TAYLOR J F, HE J, QIAN C S, QIU J S, SIMPSON B, WALKER J, BAUCK S. Optimal design of low-density SNP arrays for genomic prediction: algorithm and applications., 2016, 11(9): e0161719.

[30] ABDI H, WILLIAMS L J. Principal component analysis., 2010, 2(4): 433-459.

[31] Wu X L, Xu J Q, Feng G F, Wiggans G R, Taylor J F, He J, Qian C S, Qiu J S, Simpson B, Walker J, Bauck S. Optimal design of low-density SNP arrays for genomic prediction: algorithm and applications., 2016, 11(9): e0161719.

Using Restricted Standardized Linear Regression Model to Estimate Genomic Breed Composition in Composite Breed Animals

HE Jun1, LI Zhi1,2, Wu XiaoLin1,2

(1College of Animal Science and Technology, Hunan Agricultural University, Changsha 410128, China;2Biostatistics and Bioinformatics,Neogen GeneSeek, Lincoln, NE 68504, USA)

【】A composite breed is made up of two or more purebreds (ancestries), designed to combine advantageous genetic characteristics from the ancestry breeds and to retain heterosis in future generations without crossbreeding. Unlike crossbred populations, composite variety can be maintained as a purebred. In practice, knowing the ratio of genomic contribution of an ancestry breed to individual composite animals, referred to as the genomic breed composition (GBC), is of importance in animal breed registration, tracing breeding history and population structure, breed conservation, and the prediction of heterosis. Using a set of genomic SNP genotype and an appropriate statistical model, GBC of a purebred or crossbred animal can be estimated. So far, studies on statistical methods devote to the estimation of GBC in composite breed are limited. Linear regression (LR) analysis was commonly used to estimated GBC of individual animals, but it had some limitations such as the coefficients of ancestral breeds does not add to 1.【】The purpose of the present study was to propose and evaluate the use of restricted standardized regression analysis, as an improved approach of linear regression analysis to estimate GBC in composite animals. 【】The dataset consisted of 4 323 Beefmaster cattle and purebred animals belonging to their ancestry breeds, namely Brahman, Hereford and Shorthorn. All these animals were genotyped by GeneSeek Genomic Profiling (GGP) bovine 50K SNP chips. Allelic frequencies of each SNP and the Euclidean distance between breeds were computed for the four animal populations, and their genetic relationships were revealed by Hierarchical Clustering based on Euclidean distance of SNP allele frequencies among the four populations. Genomic breed composition of the 4 323 Beefmaster cattle were estimated using RSLR and LR, respectively, based on 7 SNP panels(1K, 5K, 10K, 20K, 30K, 40K, and all the common 47 900 SNP). 【】The results of the clustering analysis agreed well with the genetic relationships of Beefmaster and the three ancestral breeds, showing that Beefmaster was more related to Brahman than Herdford and Shorhorn. Linear regression analysis underestimated the genomic contribution ratios of Brahman cattle (0.459-0.462) and shorthorn cattle (0.208-0.212) and at the same time overestimated that of Hereford cattle (0.326-0.333) to Beefmaster cattle. In contrast, estimated GBC of the 4 323 Beefmaster cattle obtained by using RSLR agreed well with expected genomic contribution ratios of the three ancestry breeds, which were 0.497-0.503 for Brahman, 0.262-0.274 for Hereford, and 0.229-0.231 for Shorthorn, respectively. Furthermore, the standard deviations (SD) and coefficients of variance (CV) of GBC obtained by using LR were larger than those obtained using RSLR. With 20K or more SNPs as the reference panels, the SD of GBC estimated by using LR were 0.048 (Brahman), 0.032 (Hereford) and 0.051-0.052 (Shorthorn), and the corresponding CV were 10.46%-10.50% (Brahman), 9.61%-9.76% (Hereford) and 23.94%-25.00% (Shorthorn), respectively. Using RSLR, on the other hand, the SD of GBC pertaining to each of the three ancestry breeds were 0.021 (Brahman), 0.021-0.022(Hereford) and 0.024-0.025 (Shorthorn), and the responding CV were 4.18%-4.20% (Brahman), 7.89%-8.33% (Hereford) and 10.26%-10.68% (Shorthorn), correspondingly. 【】The RSLR method provided more accurate and consistent estimates of GBC in the 4 323 Beefmaster cattle than the LR approach. It thus provided a new statistical method for the estimation of GBC in composite animals.

SNP chip; linear regression; composite breeds; genomic breed composition

10.3864/j.issn.0578-1752.2020.01.018

2019-03-01；

2019-05-30

湖南省科技計(jì)劃重點(diǎn)項(xiàng)目（2018NK2081）、長沙市科技計(jì)劃重點(diǎn)項(xiàng)目（kq1801014）、湖南省百人計(jì)劃項(xiàng)目和湖南省畜禽安全協(xié)同創(chuàng)新中心項(xiàng)目

何俊，Tel：0731-84618176；E-mail：hejun@hunau.edu.cn

（責(zé)任編輯 林鑒非）