Load data

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: xlsxjars
all_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] 0
mp<-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] 0
rownames(dts)<-dts$id
sum(rownames(dts)!=dts$Animal)
## [1] 0
pen_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")

Start to prepare data for gwa and gblup

library(regress)
library(gwaR)

load("Z_and_G_matrices.Rdata")

#----------------------------
sum(!colnames(Z)%in%rownames(map))
## [1] 0
map<-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.

now create a dataframe with analysis details and run gwa

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_Pen
gb<-gblup(cont_struc[i,][[1]],dts,
          c(cont_struc[i,][[2]]),
          G,pos=c(T,T,T))
i=2
gb2<-gblup(cont_struc[i,][[1]],dts,
          c(cont_struc[i,][[2]]),
          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]]),
             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]]),
             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]]),
             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: abind
source("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 valid
rs3<-runbivar(gb,gb3,G)
rs4<-runbivar(gb,gb4,G)
rs5<-runbivar(gb,gb5,G)
rs6<-runbivar(gb3,gb5,G)

varcomp(gb)
##     Estimate   StdError  prop.var         se
## G  0.1276567 0.02565998 0.3073761 0.05533503
## In 0.2876544 0.01914469 0.6926239 0.06597559
varcomp(gb2)
##      Estimate   StdError  prop.var         se
## G  0.05351716 0.01661001 0.1515546 0.04492852
## In 0.29960411 0.01726594 0.8484454 0.06102665
check_cor(gb,gb2)
##        gen       resi       phen correcphen 
##  0.7926603  0.6779212  0.7614616  0.6982741
print.varcov(varcov(rs2))
##             gen       res zef     pheno
## Var1 0.10554635 0.2992527   0 0.4047991
## var2 0.07596806 0.2882853   0 0.3642534
## Cov  0.07596805 0.1959486   0 0.2719167
## Cor  0.84838675 0.6671321   0 0.7081320
varcomp(gb3)
##     Estimate   StdError  prop.var         se
## G  0.1458432 0.03149922 0.2697796 0.05294572
## In 0.3947582 0.02503755 0.7302204 0.06504981
check_cor(gb,gb3)
##        gen       resi       phen correcphen 
##  0.5121797  0.1709223  0.2506977  0.2523813
print.varcov(varcov(rs3))
##            gen        res zef     pheno
## Var1 0.1176432 0.29405572   0 0.4116990
## var2 0.1391285 0.39948698   0 0.5386155
## Cov  0.1045987 0.02085281   0 0.1254516
## Cor  0.8175887 0.06084131   0 0.2664077
varcomp(gb4)
##      Estimate   StdError prop.var         se
## G  0.09406119 0.02606037 0.180582 0.04723531
## In 0.42681695 0.02517325 0.819418 0.06252186
check_cor(gb,gb4)
##          gen         resi         phen   correcphen 
##  0.136650097 -0.027626040  0.157885634 -0.003810716
print.varcov(varcov(rs4))
##             gen         res zef        pheno
## Var1 0.12564521  0.28882257   0  0.414467782
## var2 0.09198137  0.42818901   0  0.520170380
## Cov  0.01150882 -0.01717796   0 -0.005669147
## Cor  0.10705519 -0.04884705   0 -0.012209547
varcomp(gb5)
##     Estimate   StdError  prop.var         se
## G  0.1862763 0.04663343 0.4051798 0.09274816
## In 0.2734612 0.03309471 0.5948202 0.10259819
check_cor(gb,gb5)
##        gen       resi       phen correcphen 
##  0.4023669  0.2368411  0.2595509  0.3147959
print.varcov(varcov(rs5))
##            gen        res zef     pheno
## Var1 0.1267288 0.28800480   0 0.4147336
## var2 0.1937629 0.27195416   0 0.4657170
## Cov  0.1126972 0.04381003   0 0.1565073
## Cor  0.7191847 0.15654020   0 0.3561139
check_cor(gb3,gb5)
##        gen       resi       phen correcphen 
##  0.2727429  0.1104116  0.1994008  0.1770295
print.varcov(varcov(rs6))
##            gen          res zef     pheno
## Var1 0.1455564  0.394520001   0 0.5400764
## var2 0.1918530  0.269743894   0 0.4615969
## Cov  0.1039650 -0.001544696   0 0.1024203
## Cor  0.6221384 -0.004735138   0 0.2051289