九色国产,午夜在线视频,新黄色网址,九九色综合,天天做夜夜做久久做狠狠,天天躁夜夜躁狠狠躁2021a,久久不卡一区二区三区

打開APP
userphoto
未登錄

開通VIP，暢享免費(fèi)電子書等14項(xiàng)超值服

開通VIP

首頁

好書

留言交流

下載APP

聯(lián)系客服
全基因組選擇介紹及實(shí)踐-2：構(gòu)建H矩陣

1, 編者自語

H矩陣作為一步法的入門技術(shù), 是需要掌握的, 本文以一篇文獻(xiàn)為例, 介紹如何從頭構(gòu)建H矩陣. 文章包括H矩陣推導(dǎo)過程和代碼實(shí)現(xiàn).

2, H矩陣定義

基因組選擇中, GBLUP的一個(gè)挑戰(zhàn)是, 在參考群構(gòu)建時(shí), 需要兩步, 第一步根據(jù)系譜和表型數(shù)據(jù), 計(jì)算出偽數(shù)據(jù)(pseudo-data)(比如, 根據(jù)系譜計(jì)算公牛的女兒產(chǎn)奶偏差作為表型值, 因?yàn)楣]有產(chǎn)奶數(shù)據(jù)), 然后用基因組信息進(jìn)行評(píng)估建模, 這就造成信息的損耗和偏離. 解決的方法是, 可以通過一種手段, 將系譜關(guān)系A(chǔ)矩陣和基因組信息構(gòu)建的親緣關(guān)系G矩陣合并為H矩陣, 這樣就成了一步法(Single-setp).

As not all animals can be genotyped, a 2- or 3-step procedure has to be followed; first, a regular genetic evaluation is run; then, corrected phenotypes or pseudo-data are used in the second step, where the marker-assisted selection model is effectively applied (Guillaume et al., 2008; VanRaden et al., 2009). These phenotypes are daughter yield deviations (DYD) and yield deviations (YD) for dairy cattle.

3, H矩陣推導(dǎo)過程—三組進(jìn)行推導(dǎo)

假設(shè)A矩陣, 包括三部分:

A is the numerator relationship matrix based on pedigree. Consider three types of animals in u: 1) ungenotyped ancestors with breeding values u1; 2) genotyped animals, with breeding values u2 (no ancestor is genotyped and phantom parents can be generated if necessary); and 3) ungenotyped animals with breeding values u3, which might descend from either one of the three types of animals.

  • 編號(hào)1 是沒有測序的個(gè)體, 祖先

  • 編號(hào)2 是測序的個(gè)體(當(dāng)代和后代)

  • 編號(hào)3 是沒有測序個(gè)體(當(dāng)代和后代)

如果所有個(gè)體A矩陣按照世代排, 那么可以將其分為如下部分:

其中是測序個(gè)體構(gòu)成的矩陣, 相應(yīng)個(gè)體對(duì)應(yīng)的是G矩陣,

是所有個(gè)體的親緣關(guān)系矩陣, 如果在所有個(gè)體矩陣中剖分出A矩陣, 那么:

為3組育種值, 它的構(gòu)成為:


那么3組的方差, 3組和1組, 3組和2組的協(xié)方差為:


那么矩陣變?yōu)?

剖出A矩陣

4, H矩陣推導(dǎo)過程—二組進(jìn)行推導(dǎo)

那么A矩陣和A逆矩陣可以剖分為:

  • 編號(hào)1為非測序個(gè)體

  • 編號(hào)2為測序個(gè)體

測序個(gè)體和非測序個(gè)體的方差以及協(xié)方差:

假定H矩陣為所有個(gè)體的矩陣:
那么H矩陣可以分為:

5, H逆矩陣推導(dǎo)過程

H矩陣還可以寫為:

對(duì)其進(jìn)行求逆:

這里$A22^{-1}$為測序個(gè)體的親緣關(guān)系逆矩陣

6, 系譜數(shù)據(jù)構(gòu)建A矩陣

R語言生成系譜

ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) ped_full

結(jié)果:

> ped_full  ID Sire Dam 1  9    1   2 2 10    3   4 3 11    5   6 4 12    7   8 5 13    9  10 6 14   11  12 7 15   11   4 8 16   13  15 9 17   13  14

利用nadiv包計(jì)算A矩陣
因?yàn)閚adiv中的makeA函數(shù), 計(jì)算的A矩陣行名沒有安裝1~17編號(hào), 需要生成一個(gè)id_r的編號(hào)對(duì)矩陣進(jìn)行排序, 排序后的矩陣按照1 ~ 17編號(hào).

ped = nadiv::prepPed(ped_full) A_mat = as.matrix(nadiv::makeA(ped)) id = row.names(A_mat) id_r = match(1:17,id) A_mat_sort = A_mat[id_r,id_r] Matrix::Matrix(A_mat_sort, sparse=TRUE) # 查看稀疏矩陣

結(jié)果:

> Matrix::Matrix(A_mat_sort, sparse=TRUE) 17 x 17 sparse Matrix of class "dsCMatrix"   [[ suppressing 17 column names '1’, '2’, '3’ ... ]] 1  1.000 .     .     .     .     .     .     .     0.50 .     .    .    0.2500 .     .      0.12500 0.12500 2  .     1.000 .     .     .     .     .     .     0.50 .     .    .    0.2500 .     .      0.12500 0.12500 3  .     .     1.000 .     .     .     .     .     .    0.500 .    .    0.2500 .     .      0.12500 0.12500 4  .     .     .     1.000 .     .     .     .     .    0.500 .    .    0.2500 .     0.5000 0.37500 0.12500 5  .     .     .     .     1.000 .     .     .     .    .     0.50 .    .      0.250 0.2500 0.12500 0.12500 6  .     .     .     .     .     1.000 .     .     .    .     0.50 .    .      0.250 0.2500 0.12500 0.12500 7  .     .     .     .     .     .     1.000 .     .    .     .    0.50 .      0.250 .      .       0.12500 8  .     .     .     .     .     .     .     1.000 .    .     .    0.50 .      0.250 .      .       0.12500 9  0.500 0.500 .     .     .     .     .     .     1.00 .     .    .    0.5000 .     .      0.25000 0.25000 10 .     .     0.500 0.500 .     .     .     .     .    1.000 .    .    0.5000 .     0.2500 0.37500 0.25000 11 .     .     .     .     0.500 0.500 .     .     .    .     1.00 .    .      0.500 0.5000 0.25000 0.25000 12 .     .     .     .     .     .     0.500 0.500 .    .     .    1.00 .      0.500 .      .       0.25000 13 0.250 0.250 0.250 0.250 .     .     .     .     0.50 0.500 .    .    1.0000 .     0.1250 0.56250 0.50000 14 .     .     .     .     0.250 0.250 0.250 0.250 .    .     0.50 0.50 .      1.000 0.2500 0.12500 0.50000 15 .     .     .     0.500 0.250 0.250 .     .     .    0.250 0.50 .    0.1250 0.250 1.0000 0.56250 0.18750 16 0.125 0.125 0.125 0.375 0.125 0.125 .     .     0.25 0.375 0.25 .    0.5625 0.125 0.5625 1.06250 0.34375 17 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.25 0.250 0.25 0.25 0.5000 0.500 0.1875 0.34375 1.00000

和文章結(jié)果一致:

7, 基因組信息G矩陣

令個(gè)體9~12為測序個(gè)體, G矩陣為對(duì)角線為1, 非對(duì)角線為0.7的矩陣, 有行名和列名.

代碼

G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 G

結(jié)果

> G     9  10  11  12 9  1.0 0.7 0.7 0.7 10 0.7 1.0 0.7 0.7 11 0.7 0.7 1.0 0.7 12 0.7 0.7 0.7 1.0

8, 根據(jù)公式計(jì)算H矩陣

需要計(jì)算的元素:

這里非測序個(gè)體為1:8, 13:17, 測序個(gè)體為9:12

# 提取子集 id_g = 9:12 id_ng = c(1:8,13:17) A22 = A[id_g,id_g] A12 = A[id_ng,id_g] A21 = A[id_g,id_ng] A22 = A[id_g,id_g] iA22 = solve(A22)

根據(jù)公式, 計(jì)算H11, H12, H21, H22

# 構(gòu)建H矩陣 H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21 H12 = A12 %*% iA22 %*%G H21 = G %*% iA22 %*% A21 H22 = G H = cbind(rbind(H11,H21),rbind(H12,H22)) round(H,4)

結(jié)果

> round(H,2)      1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 1  1.00 0.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 2  0.00 1.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 3  0.18 0.18 1.00 0.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.18 0.30 0.39 4  0.18 0.18 0.00 1.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.68 0.55 0.39 5  0.18 0.18 0.18 0.18 1.00 0.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 6  0.18 0.18 0.18 0.18 0.00 1.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 7  0.18 0.18 0.18 0.18 0.18 0.18 1.00 0.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 8  0.18 0.18 0.18 0.18 0.18 0.18 0.00 1.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 9  0.50 0.50 0.35 0.35 0.35 0.35 0.35 0.35 1.00 0.70 0.70 0.70 0.85 0.70 0.52 0.69 0.77 10 0.35 0.35 0.50 0.50 0.35 0.35 0.35 0.35 0.70 1.00 0.70 0.70 0.85 0.70 0.60 0.73 0.77 11 0.35 0.35 0.35 0.35 0.50 0.50 0.35 0.35 0.70 0.70 1.00 0.70 0.70 0.85 0.68 0.69 0.77 12 0.35 0.35 0.35 0.35 0.35 0.35 0.50 0.50 0.70 0.70 0.70 1.00 0.70 0.85 0.52 0.61 0.77 13 0.43 0.43 0.43 0.43 0.35 0.35 0.35 0.35 0.85 0.85 0.70 0.70 1.35 0.70 0.56 0.96 1.02 14 0.35 0.35 0.35 0.35 0.43 0.43 0.43 0.43 0.70 0.70 0.85 0.85 0.70 1.35 0.60 0.65 1.02 15 0.26 0.26 0.18 0.68 0.34 0.34 0.26 0.26 0.52 0.60 0.68 0.52 0.56 0.60 1.18 0.87 0.58 16 0.34 0.34 0.30 0.55 0.34 0.34 0.31 0.31 0.69 0.73 0.69 0.61 0.96 0.65 0.87 1.41 0.80 17 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.77 0.77 0.77 0.77 1.02 1.02 0.58 0.80 1.52

文章結(jié)果:

9, 代碼匯總

# 構(gòu)建系譜 ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) ped_full # 計(jì)算A矩陣 ped = nadiv::prepPed(ped_full) A_mat = as.matrix(nadiv::makeA(ped)) id = row.names(A_mat) id_r = match(1:17,id) A_mat_sort = A_mat[id_r,id_r] A = A_mat_sort Matrix::Matrix(A_mat_sort, sparse=TRUE) # 構(gòu)建G矩陣 G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 G # 提取子集 id_g = 9:12 id_ng = c(1:8,13:17) A11 = A[id_ng,id_ng] A22 = A[id_g,id_g] A12 = A[id_ng,id_g] A21 = A[id_g,id_ng] A22 = A[id_g,id_g] iA22 = solve(A22) # 構(gòu)建H矩陣 H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21 H12 = A12 %*% iA22 %*%G H21 = G %*% iA22 %*% A21 H22 = G H = cbind(rbind(H11,H21),rbind(H12,H22)) id = rownames(H) id_r = match(1:17,id) H = H[id_r,id_r] round(H,2)

10, 編寫一個(gè)函數(shù)

這里編寫一個(gè)H_matrix函數(shù), 參數(shù)有:

  • G矩陣, 包括行名和列名

  • ped, 系譜文件

  • id_g, 是基因型id

  • id_full, 是所有個(gè)體的id

H_matrix = function(G = G, ped = ped,id_g, id_full){  ped = nadiv::prepPed(ped_full)  A_mat = as.matrix(nadiv::makeA(ped))  A = A_mat[id_full,id_full]  id_ng = setdiff(id_full,id_g)  A11 = A[id_ng,id_ng]  A22 = A[id_g,id_g]  A12 = A[id_ng,id_g]  A21 = A[id_g,id_ng]  A22 = A[id_g,id_g]  iA22 = solve(A22)  # 構(gòu)建H矩陣  H11 = A11 + A12 %*% iA22 %*% (G - A22) %*% iA22 %*% A21  H12 = A12 %*% iA22 %*%G  H21 = G %*% iA22 %*% A21  H22 = G  H = cbind(rbind(H11,H21),rbind(H12,H22))  id = rownames(H)  id_r = match(1:17,id)  H = H[id_r,id_r]  return(H) }

測試

# pedigree ped_full <- data.frame(ID=9:17,Sire=c(1,3,5,7,9,11,11,13,13),Dam=c(2,4,6,8,10,12,4,15,14)) # genotype G <- matrix(0.7,4,4) diag(G) <- 1 rownames(G) <- colnames(G) <- 9:12 # id_g id_g = 9:12 # id_full id_full = 1:17 H = H_matrix(G,ped_full,id_g,id_full) round(H,2)

結(jié)果

> round(H,2)      1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 1  1.00 0.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 2  0.00 1.00 0.18 0.18 0.18 0.18 0.18 0.18 0.50 0.35 0.35 0.35 0.43 0.35 0.26 0.34 0.39 3  0.18 0.18 1.00 0.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.18 0.30 0.39 4  0.18 0.18 0.00 1.00 0.18 0.18 0.18 0.18 0.35 0.50 0.35 0.35 0.43 0.35 0.68 0.55 0.39 5  0.18 0.18 0.18 0.18 1.00 0.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 6  0.18 0.18 0.18 0.18 0.00 1.00 0.18 0.18 0.35 0.35 0.50 0.35 0.35 0.43 0.34 0.34 0.39 7  0.18 0.18 0.18 0.18 0.18 0.18 1.00 0.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 8  0.18 0.18 0.18 0.18 0.18 0.18 0.00 1.00 0.35 0.35 0.35 0.50 0.35 0.43 0.26 0.31 0.39 9  0.50 0.50 0.35 0.35 0.35 0.35 0.35 0.35 1.00 0.70 0.70 0.70 0.85 0.70 0.52 0.69 0.77 10 0.35 0.35 0.50 0.50 0.35 0.35 0.35 0.35 0.70 1.00 0.70 0.70 0.85 0.70 0.60 0.73 0.77 11 0.35 0.35 0.35 0.35 0.50 0.50 0.35 0.35 0.70 0.70 1.00 0.70 0.70 0.85 0.68 0.69 0.77 12 0.35 0.35 0.35 0.35 0.35 0.35 0.50 0.50 0.70 0.70 0.70 1.00 0.70 0.85 0.52 0.61 0.77 13 0.43 0.43 0.43 0.43 0.35 0.35 0.35 0.35 0.85 0.85 0.70 0.70 1.35 0.70 0.56 0.96 1.02 14 0.35 0.35 0.35 0.35 0.43 0.43 0.43 0.43 0.70 0.70 0.85 0.85 0.70 1.35 0.60 0.65 1.02 15 0.26 0.26 0.18 0.68 0.34 0.34 0.26 0.26 0.52 0.60 0.68 0.52 0.56 0.60 1.18 0.87 0.58 16 0.34 0.34 0.30 0.55 0.34 0.34 0.31 0.31 0.69 0.73 0.69 0.61 0.96 0.65 0.87 1.41 0.80 17 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.77 0.77 0.77 0.77 1.02 1.02 0.58 0.80 1.52

搞定.

公眾號(hào)ID：R-breeding

11, 數(shù)據(jù)及結(jié)果來源

paper:Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information.[J]. Journal of Dairy Science, 2009, 92(9):4656-63.

本站僅提供存儲(chǔ)服務(wù)，所有內(nèi)容均由用戶發(fā)布，如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊舉報(bào)。
打開APP，閱讀全文并永久保存 查看更多類似文章
猜你喜歡
類似文章
論numpy中matrix 和 array的區(qū)別
OpenCV學(xué)習(xí)筆記（四十）——再談OpenCV數(shù)據(jù)結(jié)構(gòu)Mat詳解
單細(xì)胞數(shù)據(jù)挖掘?qū)崙?zhàn)：文獻(xiàn)復(fù)現(xiàn)（一）批量讀取數(shù)據(jù)
機(jī)器人如何從各軸角度算出當(dāng)前XYZ
多維數(shù)組
OpenCV學(xué)習(xí)筆記（三十六）——Kalman濾波做運(yùn)動(dòng)目標(biāo)跟蹤
更多類似文章 >>
生活服務(wù)
熱點(diǎn)新聞
分享 收藏 導(dǎo)長圖 關(guān)注 下載文章
綁定賬號(hào)成功
后續(xù)可登錄賬號(hào)暢享VIP特權(quán)！
如果VIP功能使用有故障，
可點(diǎn)擊這里聯(lián)系客服！

聯(lián)系客服