The pedirgree file contain 1037 individuals from 3 generations \(males=518\\ females=519\)
ped.msuprp<-as.data.frame(MSUPRP_meat$pedigree)
rownames(ped.msuprp)<-ped.msuprp$ID
head(ped.msuprp)
## ID Par1 Par2 gener sex
## 6070 6070 0 0 0 0
## 6071 6071 0 0 0 0
## 6086 6086 0 0 0 0
## 6088 6088 0 0 0 0
## 6092 6092 0 0 0 0
## 6323 6323 0 0 0 1
table(ped.msuprp$sex)
##
## 0 1
## 518 519
table(ped.msuprp$gener)
##
## 0 1 2
## 20 57 960
dim(G)
## [1] 1015 1015
length(intersect(rownames(G), rownames(pheno.msuprp)))
## [1] 940
length(intersect(ped.msuprp$ID, rownames(pheno.msuprp)))
## [1] 960
length(intersect(rownames(G), ped.msuprp$ID))
## [1] 1015
length(intersect(rownames(G), rownames(ped.msuprp)))
## [1] 1015
# Individual genotyped (Generations)
idped.gen<-ped.msuprp[rownames(G),]
table(idped.gen$gener)
##
## 0 1 2
## 20 55 940
colnames(idped.gen)[1]<-"id"
colnames(idped.gen)[5]<-"Sex"
phenofile<-pheno.msuprp%>%left_join(idped.gen)
## Joining, by = "id"
table(phenofile$gener)
##
## 2
## 940
It were selected 12 meat quality and carcass traits recorded in 960 individuals, which were genotyping for 42428 SNP. For the statistical analysis was fitted a linear mixed model described in Cariró 2017, in which included fixed effects of \(sex\), \(slaughter group\) and \(carcass\ weight\), with exception for the number of ribs trait, that only had sex as fixed effect, the genomic relationship matrix (G) was computed as VanRadem 2008.
\(\mathbf{Y= X \beta + Zu +e}\)
# 2. Gblup analysis meat traits
cont_struc2<-rbind(
c("last_lum.1",y~sex + car_wt + slgdt_cd),
c("dress_ptg.1",y~sex + car_wt + slgdt_cd),
c("car_length.1",y~sex + car_wt + slgdt_cd),
c("num_ribs.1",y~sex),
c("belly.1",y~sex + car_wt + slgdt_cd),
c("juiciness.1",y~sex + car_wt + slgdt_cd),
c("WBS.1",y~sex + car_wt + slgdt_cd),
c("ph_24h.1",y~sex + car_wt + slgdt_cd),
c("driploss.1",y~sex + car_wt + slgdt_cd),
c("protein.1",y~sex + car_wt + slgdt_cd),
c("cook_yield.1",y~sex + car_wt + slgdt_cd),
c("bf10_22wk",y~sex + car_wt + slgdt_cd)
)
gb.meatraits<-list()
for (i in 1:nrow(cont_struc2)) {
gb.meatraits[[i]]<-gblup(cont_struc2[i,][[1]],pheno.msuprp,
c(cont_struc2[i,][[2]]),
G,pos=c(T,T))
names(gb.meatraits)[i]<-cont_struc2[i]
}
## Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid
#save(gb.meatraits, file = "gbMeatratisMsuPRP.Rdata")
lapply(gb.meatraits, function(x)varcomp(x))
## $last_lum.1
## Estimate StdError prop.var se
## G 11.94571 2.109141 0.4079896 0.05732552
## In 17.33374 1.178952 0.5920104 0.06185818
##
## $dress_ptg.1
## Estimate StdError prop.var se
## G 0.8023418 0.1777845 0.2824343 0.05422069
## In 2.0384664 0.1276088 0.7175657 0.06414282
##
## $car_length.1
## Estimate StdError prop.var se
## G 1.914497 0.3014001 0.4848341 0.05739534
## In 2.034271 0.1461378 0.5151659 0.05816738
##
## $num_ribs.1
## Estimate StdError prop.var se
## G 0.2767628 0.06142882 0.3580153 0.06718035
## In 0.4962848 0.03999038 0.6419847 0.07587294
##
## $belly.1
## Estimate StdError prop.var se
## G 0.02476174 0.006784861 0.1982324 0.04957400
## In 0.10015092 0.005952166 0.8017676 0.06240203
##
## $juiciness.1
## Estimate StdError prop.var se
## G 0.02367782 0.01209957 0.07286584 0.03631702
## In 0.30127313 0.01644872 0.92713416 0.05148086
##
## $WBS.1
## Estimate StdError prop.var se
## G 0.1185178 0.02584294 0.2925219 0.05496538
## In 0.2866409 0.01820819 0.7074781 0.06459903
##
## $ph_24h.1
## Estimate StdError prop.var se
## G 0.002886079 0.0008285234 0.1868492 0.04932895
## In 0.012559954 0.0007510999 0.8131508 0.06273813
##
## $driploss.1
## Estimate StdError prop.var se
## G 0.3073669 0.06849498 0.2803005 0.05420036
## In 0.7891952 0.04941001 0.7196995 0.06422381
##
## $protein.1
## Estimate StdError prop.var se
## G 0.4460666 0.08317114 0.3783189 0.05738021
## In 0.7330089 0.04926695 0.6216811 0.06334544
##
## $cook_yield.1
## Estimate StdError prop.var se
## G 2.378149 0.4996554 0.3111316 0.05563172
## In 5.265398 0.3380695 0.6888684 0.06450203
##
## $bf10_22wk
## Estimate StdError prop.var se
## G 7.53506 1.388605 0.3808797 0.05692426
## In 12.24825 0.816722 0.6191203 0.06270037
lapply(gb.meatraits, function(x)dim(x$model$G))
## $last_lum.1
## [1] 929 929
##
## $dress_ptg.1
## [1] 931 931
##
## $car_length.1
## [1] 930 930
##
## $num_ribs.1
## [1] 655 655
##
## $belly.1
## [1] 930 930
##
## $juiciness.1
## [1] 925 925
##
## $WBS.1
## [1] 920 920
##
## $ph_24h.1
## [1] 910 910
##
## $driploss.1
## [1] 929 929
##
## $protein.1
## [1] 918 918
##
## $cook_yield.1
## [1] 921 921
##
## $bf10_22wk
## [1] 931 931
n<-lapply(gb.meatraits, function(x)length(x$model$y))%>%map_df(as_tibble)
colnames(n)<-"n"
load("A_matrix_MSUPRP_Meatq.Rdata")
dim(A)
## [1] 1037 1037
A<-A[rownames(G),rownames(G)]
dim(A)
## [1] 1015 1015
dim(G)
## [1] 1015 1015
gbAmat.meatraits<-list()
for (i in 1:nrow(cont_struc2)) {
gbAmat.meatraits[[i]]<-gblup(cont_struc2[i,][[1]],pheno.msuprp,
c(cont_struc2[i,][[2]]),
G = A ,pos=c(T,T))
names(gbAmat.meatraits)[i]<-cont_struc2[i]
}
## Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid
#save(gbAmat.meatraits, file = "gbAmat_MsuPRPmeat.Rdata")
lapply(gbAmat.meatraits, function(x)varcomp(x)[1,c(3,4)])
## $last_lum.1
## prop.var se
## G 0.2564387 0.09229331
##
## $dress_ptg.1
## prop.var se
## G 0.1793848 0.07670021
##
## $car_length.1
## prop.var se
## G 0.3846533 0.1157862
##
## $num_ribs.1
## prop.var se
## G 0.4989639 0.147774
##
## $belly.1
## prop.var se
## G 0.1504623 0.070337
##
## $juiciness.1
## prop.var se
## G 0.06254925 0.04720592
##
## $WBS.1
## prop.var se
## G 0.2429087 0.09004674
##
## $ph_24h.1
## prop.var se
## G 0.1646012 0.07413536
##
## $driploss.1
## prop.var se
## G 0.4197666 0.1222509
##
## $protein.1
## prop.var se
## G 0.2020227 0.08197036
##
## $cook_yield.1
## prop.var se
## G 0.344088 0.1088444
##
## $bf10_22wk
## prop.var se
## G 0.1508952 0.07054918
lapply(gbAmat.meatraits, function(x)dim(x$model$G))
## $last_lum.1
## [1] 929 929
##
## $dress_ptg.1
## [1] 931 931
##
## $car_length.1
## [1] 930 930
##
## $num_ribs.1
## [1] 655 655
##
## $belly.1
## [1] 930 930
##
## $juiciness.1
## [1] 925 925
##
## $WBS.1
## [1] 920 920
##
## $ph_24h.1
## [1] 910 910
##
## $driploss.1
## [1] 929 929
##
## $protein.1
## [1] 918 918
##
## $cook_yield.1
## [1] 921 921
##
## $bf10_22wk
## [1] 931 931
hG<-lapply(gb.meatraits, function(x)varcomp(x)[1,c(3,4)])%>%
map_df(as_tibble)%>%mutate(trait=names(gb.meatraits))%>%select(trait, prop.var, se)
hG<-cbind(n,hG)%>%select(trait,n,prop.var,se)
hA<-lapply(gbAmat.meatraits, function(x)varcomp(x)[1,c(3,4)])%>%map_df(as_tibble)%>%
mutate(trait=names(gbAmat.meatraits))%>%select(trait, prop.var, se)
h2<-left_join(hA,hG, "trait")
h2<-h2%>%select(trait,n,prop.var.x:se.y)
colnames(h2)<-c("Trait","n","Heritability","se",
"Heritability", "se")
kable(h2)%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),full_width = F,
position = "center",font_size = 14) %>%
add_header_above(c(" ", " ","A matrix" = 2, "G matrix" = 2), italic = T,
bold = T)
| Trait | n | Heritability | se | Heritability | se |
|---|---|---|---|---|---|
| last_lum.1 | 929 | 0.2564387 | 0.0922933 | 0.4079896 | 0.0573255 |
| dress_ptg.1 | 931 | 0.1793848 | 0.0767002 | 0.2824343 | 0.0542207 |
| car_length.1 | 930 | 0.3846533 | 0.1157862 | 0.4848341 | 0.0573953 |
| num_ribs.1 | 655 | 0.4989639 | 0.1477740 | 0.3580153 | 0.0671804 |
| belly.1 | 930 | 0.1504623 | 0.0703370 | 0.1982324 | 0.0495740 |
| juiciness.1 | 925 | 0.0625493 | 0.0472059 | 0.0728658 | 0.0363170 |
| WBS.1 | 920 | 0.2429087 | 0.0900467 | 0.2925219 | 0.0549654 |
| ph_24h.1 | 910 | 0.1646012 | 0.0741354 | 0.1868492 | 0.0493289 |
| driploss.1 | 929 | 0.4197666 | 0.1222509 | 0.2803005 | 0.0542004 |
| protein.1 | 918 | 0.2020227 | 0.0819704 | 0.3783189 | 0.0573802 |
| cook_yield.1 | 921 | 0.3440880 | 0.1088444 | 0.3111316 | 0.0556317 |
| bf10_22wk | 931 | 0.1508952 | 0.0705492 | 0.3808797 | 0.0569243 |