First we load the data and perform a few data object checks
rm(list=ls())
setwd("C:\\Users\\steibelj\\OneDrive\\Documents\\kaitlin")
library(plyr)
#--------------------
load("lesions_wide_format.Rdata")
all_lesions<-lesions_wide_format
library(xlsx)## Loading required package: rJava## Loading required package: xlsxjarsall_lesions[is.na(all_lesions)]<-NA
#Create data
ls(all_lesions)## [1] "Finisher_Post_Front_ls" "Finisher_Post_Middle_ls"
## [3] "Finisher_Post_Rear_ls" "Finisher_Post_Total_ls"
## [5] "Finisher_Pre_Front_ls" "Finisher_Pre_Middle_ls"
## [7] "Finisher_Pre_Rear_ls" "Finisher_Pre_Total_ls"
## [9] "Finisher_Stable_Front_ls" "Finisher_Stable_Middle_ls"
## [11] "Finisher_Stable_Rear_ls" "Finisher_Stable_Total_ls"
## [13] "ids" "Nursery_Post_Front_ls"
## [15] "Nursery_Post_Middle_ls" "Nursery_Post_Rear_ls"
## [17] "Nursery_Post_Total_ls" "Nursery_Pre_Front_ls"
## [19] "Nursery_Pre_Middle_ls" "Nursery_Pre_Rear_ls"
## [21] "Nursery_Pre_Total_ls" "Nursery_Stable_Front_ls"
## [23] "Nursery_Stable_Middle_ls" "Nursery_Stable_Rear_ls"
## [25] "Nursery_Stable_Total_ls" "Sow_Post_Front_ls"
## [27] "Sow_Post_Middle_ls" "Sow_Post_Rear_ls"
## [29] "Sow_Post_Total_ls" "Sow_Pre_Front_ls"
## [31] "Sow_Pre_Middle_ls" "Sow_Pre_Rear_ls"
## [33] "Sow_Pre_Total_ls" "Sow_Stable_Front_ls"
## [35] "Sow_Stable_Middle_ls" "Sow_Stable_Rear_ls"
## [37] "Sow_Stable_Total_ls"lesion_index<-grep("_ls",colnames(all_lesions))
all_lesions[lesion_index]<-log(all_lesions[lesion_index]+1)
load("covar_wide_format.Rdata")
covar_w<-covar_wide_format
sum(all_lesions$ids!=covar_w$Animal)## [1] 0mp<-read.table("SNP_Map.txt",header=T,sep="\t")
map<-mp[,c(3,4)]
rownames(map)<-mp$Name
colnames(map)<-c("chr","pos")Create a new dataframe and prepare all columns
dts<-cbind(all_lesions,covar_w)
class(dts)## [1] "data.frame"sum(dts$ids!=dts$Animal)## [1] 0rownames(dts)<-dts$id
sum(rownames(dts)!=dts$Animal)## [1] 0pen_idx<-grep("Pen",colnames(dts))
penf<-as.data.frame(lapply(dts[,pen_idx],as.factor))
ob_idx<-grep("Obs",colnames(dts))
obsf<-as.data.frame(lapply(dts[,ob_idx],as.factor))
pastef<-function(v1,v2,sep="_"){
re<-paste(v1,v2,sep=sep)
re[is.na(v1)|is.na(v2)]<-NA
re<-as.factor(re)
return(re)
}
dts[,pen_idx]<-as.data.frame(lapply(penf,pastef,v2=dts$Rep))
sapply(dts[,pen_idx],levels)## $Nursery_Pen
## [1] "ML1_2" "ML1_3" "ML1_4" "ML1_5" "ML1_6" "ML1_7" "ML1_8" "ML1_9"
## [9] "ML2_2" "ML2_3" "ML2_4" "ML2_5" "ML2_6" "ML2_7" "ML2_8" "ML2_9"
## [17] "ML3_2" "ML3_3" "ML3_4" "ML3_5" "ML3_6" "ML3_7" "ML3_8" "ML3_9"
## [25] "ML4_2" "ML4_3" "ML4_4" "ML4_5" "ML4_6" "ML4_7" "ML4_8" "ML4_9"
## [33] "ML5_2" "ML5_3" "ML5_4" "ML5_5" "ML5_6" "ML5_7" "ML5_8" "ML5_9"
## [41] "ML6_2" "ML6_3" "ML6_4" "ML6_5" "ML6_6" "ML6_7" "ML6_8" "ML6_9"
## [49] "ML7_3" "ML7_4" "ML7_5" "ML7_6" "ML7_7" "ML7_8" "ML7_9" "MR1_2"
## [57] "MR1_3" "MR1_4" "MR1_5" "MR1_6" "MR1_7" "MR1_8" "MR1_9" "MR2_2"
## [65] "MR2_3" "MR2_4" "MR2_5" "MR2_6" "MR2_7" "MR2_8" "MR2_9" "MR3_2"
## [73] "MR3_3" "MR3_4" "MR3_5" "MR3_6" "MR3_7" "MR3_8" "MR3_9" "MR4_2"
## [81] "MR4_3" "MR4_4" "MR4_5" "MR4_6" "MR4_7" "MR4_8" "MR4_9" "MR5_2"
## [89] "MR5_3" "MR5_4" "MR5_5" "MR5_6" "MR5_7" "MR5_8" "MR5_9" "MR6_2"
## [97] "MR6_3" "MR6_4" "MR6_5" "MR6_6" "MR6_7" "MR6_8" "MR6_9" "MR7_3"
## [105] "MR7_4" "MR7_5" "MR7_6" "MR7_7" "MR7_8" "MR7_9"
##
## $Finisher_Pen
## [1] "10_2" "10_3" "10_4" "10_5" "10_6" "10_7" "10_8" "10_9" "11_2" "11_3"
## [11] "11_4" "11_5" "11_6" "11_7" "11_8" "11_9" "12_2" "12_3" "12_4" "12_5"
## [21] "12_6" "12_7" "12_8" "12_9" "13_2" "13_3" "13_4" "13_5" "13_6" "13_7"
## [31] "13_8" "13_9" "2_2" "2_3" "2_4" "2_5" "2_6" "2_7" "2_8" "2_9"
## [41] "3_2" "3_3" "3_4" "3_5" "3_6" "3_7" "3_8" "3_9" "4_2" "4_3"
## [51] "4_4" "4_5" "4_6" "4_7" "4_8" "4_9" "5_2" "5_3" "5_4" "5_5"
## [61] "5_6" "5_7" "5_8" "5_9" "6_3" "6_4" "6_5" "6_6" "6_7" "6_8"
## [71] "6_9" "9_3" "9_4" "9_5" "9_6" "9_7" "9_8" "9_9"
##
## $Sow_Pen
## [1] "10_2" "10_6" "10_7" "10_8" "11_3" "11_4" "11_5" "11_9" "12_2" "12_6"
## [11] "12_7" "12_8" "13_3" "13_4" "13_5" "13_9" "2_2" "2_6" "2_7" "2_8"
## [21] "3_3" "3_4" "3_5" "3_7" "3_9" "4_2" "4_3" "4_6" "4_7" "4_8"
## [31] "5_3" "5_4" "5_5" "5_9" "6_6" "6_7" "6_8" "9_3" "9_4" "9_5"
## [41] "9_9"dts[,ob_idx]<-obsf
sapply(dts[,ob_idx],levels)## $Finisher_Post_Obs
## [1] "CO" "JS" "KW"
##
## $Finisher_Pre_Obs
## [1] "CO" "JS" "KW"
##
## $Finisher_Stable_Obs
## [1] "CO" "JS" "KW"
##
## $Nursery_Post_Obs
## [1] "CO" "JS" "KW"
##
## $Nursery_Pre_Obs
## [1] "#N/A" "CO" "JS" "KW"
##
## $Nursery_Stable_Obs
## [1] "CO" "JS" "KW"
##
## $Sow_Post_Obs
## [1] "CO" "JS" "KW"
##
## $Sow_Pre_Obs
## [1] "CO" "JS" "KW"
##
## $Sow_Stable_Obs
## [1] "CO" "JS" "KW"dts[is.na(dts)]<-NA
dts$Rep<-as.factor(dts$Rep)
na<-rowSums(is.na(dts))#
#dts<-dts[na==0,]
lsdata<-dts
save(lsdata,file="lession_plus.Rdata")library(regress)
library(gwaR)
load("Z_and_G_matrices.Rdata")
#----------------------------
sum(!colnames(Z)%in%rownames(map))## [1] 0map<-map[colnames(Z),]
to_include<-as.character(map$chr)%in%(1:18)
map<-map[to_include,]
Z<-Z[,to_include]
map$chr<-droplevels(map$chr)
#Notice an issue with the SNP: many unpositioned SNP and even a Y-chr SNP!
#unmapped SNP are OK for computing G but they should be dropped from GWA
#Y-snp and X-SNP COULD be included in GWA. not in G.cont_struc<-rbind(
c("Finisher_Post_Front_ls",y~Sex+Rep+Finisher_Pre_Front_ls+Finisher_Pre_Obs:Finisher_Post_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Post_Middle_ls",y~Sex+Rep+Finisher_Pre_Middle_ls+Finisher_Pre_Obs:Finisher_Post_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Post_Rear_ls",y~Sex+Rep+Finisher_Pre_Rear_ls+Finisher_Pre_Obs:Finisher_Post_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Post_Total_ls",y~Sex+Rep+Finisher_Pre_Total_ls+Finisher_Pre_Obs:Finisher_Post_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Stable_Front_ls",y~Sex+Rep+Finisher_Stable_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Stable_Middle_ls",y~Sex+Rep+Finisher_Stable_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Stable_Rear_ls",y~Sex+Rep+Finisher_Stable_Obs+Fin_Wt,~Finisher_Pen),
c("Finisher_Stable_Total_ls",y~Sex+Rep+Finisher_Stable_Obs+Fin_Wt,~Finisher_Pen),
c("Nursery_Post_Front_ls",y~Sex+Rep+Nursery_Pre_Front_ls+Nursery_Pre_Obs:Nursery_Post_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Post_Middle_ls",y~Sex+Rep+Nursery_Pre_Middle_ls+Nursery_Pre_Obs:Nursery_Post_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Post_Rear_ls",y~Sex+Rep+Nursery_Pre_Rear_ls+Nursery_Pre_Obs:Nursery_Post_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Post_Total_ls",y~Sex+Rep+Nursery_Pre_Total_ls+Nursery_Pre_Obs:Nursery_Post_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Stable_Front_ls",y~Sex+Rep+Nursery_Stable_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Stable_Middle_ls",y~Sex+Rep+Nursery_Stable_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Stable_Rear_ls",y~Sex+Rep+Nursery_Stable_Obs+Nurs_Wt,~Nursery_Pen),
c("Nursery_Stable_Total_ls",y~Sex+Rep+Nursery_Stable_Obs+Nurs_Wt,~Nursery_Pen),
c("Sow_Post_Front_ls",y~Rep+Sow_Pre_Front_ls+Sow_Pre_Obs:Sow_Post_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Post_Middle_ls",y~Rep+Sow_Pre_Middle_ls+Sow_Pre_Obs:Sow_Post_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Post_Rear_ls",y~Rep+Sow_Pre_Rear_ls+Sow_Pre_Obs:Sow_Post_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Post_Total_ls",y~Rep+Sow_Pre_Total_ls+Sow_Pre_Obs:Sow_Post_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Stable_Front_ls",y~Rep+Sow_Stable_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Stable_Middle_ls",y~Rep+Sow_Stable_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Stable_Rear_ls",y~Rep+Sow_Stable_Obs+Sow_Wt,~Sow_Pen),
c("Sow_Stable_Total_ls",y~Rep+Sow_Stable_Obs+Sow_Wt,~Sow_Pen)
)
i=1
cont_struc[1,]## [[1]]
## [1] "Finisher_Post_Front_ls"
##
## [[2]]
## y ~ Sex + Rep + Finisher_Pre_Front_ls + Finisher_Pre_Obs:Finisher_Post_Obs +
## Fin_Wt
##
## [[3]]
## ~Finisher_Pengb<-gblup(cont_struc[i,][[1]],dts,
c(cont_struc[i,][[2]],cont_struc[i,][[3]]),
G,pos=c(T,T,T))
i=2
gb2<-gblup(cont_struc[i,][[1]],dts,
c(cont_struc[i,][[2]],cont_struc[i,][[3]]),
G,pos=c(T,T,T))
i=9
cont_struc[9,]## [[1]]
## [1] "Nursery_Post_Front_ls"
##
## [[2]]
## y ~ Sex + Rep + Nursery_Pre_Front_ls + Nursery_Pre_Obs:Nursery_Post_Obs +
## Nurs_Wt
##
## [[3]]
## ~Nursery_Pen gb3<-gblup(cont_struc[i,][[1]],dts,
c(cont_struc[i,][[2]],cont_struc[i,][[3]]),
G,pos=c(T,T,T))
i=5
cont_struc[i,]## [[1]]
## [1] "Finisher_Stable_Front_ls"
##
## [[2]]
## y ~ Sex + Rep + Finisher_Stable_Obs + Fin_Wt
##
## [[3]]
## ~Finisher_Pen gb4<-gblup(cont_struc[i,][[1]],dts,
c(cont_struc[i,][[2]],cont_struc[i,][[3]]),
G,pos=c(T,T,T))
i=17
cont_struc[i,]## [[1]]
## [1] "Sow_Post_Front_ls"
##
## [[2]]
## y ~ Rep + Sow_Pre_Front_ls + Sow_Pre_Obs:Sow_Post_Obs + Sow_Wt
##
## [[3]]
## ~Sow_Pen gb5<-gblup(cont_struc[i,][[1]],dts,
c(cont_struc[i,][[2]],cont_struc[i,][[3]]),
G,pos=c(T,T,T))
cont_struc[c(1,2,5,9,17),1] ## [[1]]
## [1] "Finisher_Post_Front_ls"
##
## [[2]]
## [1] "Finisher_Post_Middle_ls"
##
## [[3]]
## [1] "Finisher_Stable_Front_ls"
##
## [[4]]
## [1] "Nursery_Post_Front_ls"
##
## [[5]]
## [1] "Sow_Post_Front_ls"library(magic)## Loading required package: abindsource("build_biv.R")
source("bivar_func(updated5.4.16).R")
rs2<-runbivar(gb,gb2,G)## Warning: solution lies close to zero for some positive variance components, their standard errors may not be validrs3<-runbivar(gb,gb3,G,commonr = F)
rs4<-runbivar(gb,gb4,G)
rs5<-runbivar(gb,gb5,G,commonr = F)
rs6<-runbivar(gb3,gb5,G,commonr = F)
varcomp(gb)## Estimate StdError prop.var se
## G 0.13495722 0.02603022 0.31949855 0.05514736
## Finisher_Pen 0.02495741 0.00892020 0.05908433 0.02068122
## In 0.26248858 0.01862037 0.62141711 0.06361812varcomp(gb2)## Estimate StdError prop.var se
## G 0.05529746 0.01601668 0.1531418 0.04252140
## Finisher_Pen 0.05283027 0.01298733 0.1463091 0.03250257
## In 0.25295895 0.01554128 0.7005491 0.05804723check_cor(gb,gb2)## gen resi phen correcphen
## 0.8036329 0.6776376 0.7614616 0.6984050print.varcov(varcov(rs2))## gen res zef pheno
## Var1 0.10643762 0.2770293 0.02544789 0.4089148
## var2 0.07512811 0.2434068 0.04848427 0.3670192
## Cov 0.07505341 0.1722213 0.02542886 0.2727036
## Cor 0.83930850 0.6632202 0.72393662 0.7039310varcomp(gb3)## Estimate StdError prop.var se
## G 0.14012006 0.03081159 0.25772861 0.05200779
## Nursery_Pen 0.03575283 0.01233925 0.06576166 0.02223680
## In 0.36779999 0.02448852 0.67650973 0.06379739check_cor(gb,gb3)## gen resi phen correcphen
## 0.4881363 0.1688381 0.2506977 0.2528819print.varcov(varcov(rs3,Z = list(V1=c("ZZ1"),V2=c("ZZ2"),cov=NULL)))## gen res zef pheno
## Var1 0.1288448 0.26691812 0.02364684 0.4194097
## var2 0.1353751 0.37219947 0.03406272 0.5416373
## Cov 0.1006790 0.02454693 0.00000000 0.1252259
## Cor 0.7623178 0.07787899 0.00000000 0.2627367varcomp(gb4)## Estimate StdError prop.var se
## G 0.09455502 0.02539940 0.1776062 0.04527629
## Finisher_Pen 0.06278573 0.01667261 0.1179327 0.02912678
## In 0.37504510 0.02349282 0.7044611 0.05978193check_cor(gb,gb4)## gen resi phen correcphen
## 0.114428820 -0.027075955 0.157885634 -0.002481347print.varcov(varcov(rs4))## gen res zef pheno
## Var1 0.13420155 0.26282009 0.025182641 0.422204281
## var2 0.09366933 0.37575287 0.062347142 0.531769346
## Cov 0.00302476 -0.01057457 -0.002194693 -0.009744505
## Cor 0.02697823 -0.03364979 -0.055387884 -0.020565384varcomp(gb5)## Estimate StdError prop.var se
## G 0.18216408 0.04606445 0.39455808 0.09186994
## Sow_Pen 0.01342438 0.01165887 0.02907651 0.02505188
## In 0.26610296 0.03303881 0.57636541 0.10179882check_cor(gb,gb5)## gen resi phen correcphen
## 0.4241279 0.2323130 0.2595509 0.3159587print.varcov(varcov(rs5,Z = list(V1=c("ZZ1"),V2=c("ZZ2"),cov=NULL)))## gen res zef pheno
## Var1 0.1338935 0.26238993 0.02674913 0.4230325
## var2 0.1901715 0.26369888 0.01584251 0.4697129
## Cov 0.1189742 0.04417216 0.00000000 0.1631463
## Cor 0.7455905 0.16792716 0.00000000 0.3659942check_cor(gb3,gb5)## gen resi phen correcphen
## 0.2641671 0.1093516 0.1994008 0.1768464print.varcov(varcov(rs6,Z = list(V1=c("ZZ1"),V2=c("ZZ2"),cov=NULL)))## gen res zef pheno
## Var1 0.1419197 0.365922752 0.03654532 0.5443878
## var2 0.1931754 0.257883638 0.01579557 0.4668547
## Cov 0.1045256 0.002335234 0.00000000 0.1068608
## Cor 0.6312849 0.007601929 0.00000000 0.2119696