blupADC-Genomic prediction(selection)

Single-trait BLUP model

library(blupSUP)
library(blupADC)
data_path=system.file("extdata", package = "blupSUP")
#prepare parameters the BLUP analysis
blup_par=BLUPpar$new(
phe_file=paste0(data_path,"/phenotype.txt"),                                            
phe_colnames=c("Id","Mean","Sex","Herd_Year_Season","Litter","Trait1","Trait2","Age"),
dmu_integer_n=5, #number of integer variable in the phenotype file (required by DMU)
analysis_model="PBLUP_A",
kinship_file=paste0(data_path,"/pedigree.txt"))
#define the model effect  of Trait1
myblup=BLUP$new(fixed=Trait1~Sex+Herd_Year_Season, #fixed effect
                random=~Id+Litter, #random effect 
                pars=blup_par) #blup_par has been defined previously
myblup

## <BLUP::DMU> 
## <ADCmodel::1 trait> 
## ----Trait:Trait1~Sex(fixed)+Herd_Year_Season(fixed)+Id(random)+Litter(random)

#BLUP analysis of Trait1
myblup$run(output_path="D:/test",task_type="DMU") # default using DMU 
#myblup$run(output_path="D:/test",task_type="BLUPF90") #using BLUPF90 for analysis 

#output of BLUP analysis
ebv=myblup$ebv         #breeding value
h2=myblup$h2   #heritability and its SE
user_defined_h2=myblup$cal_se(~Id1/(Id1+Litter1+Res1)) #user defined expression of calculating h2 based on delta method

#GBLUP-A analysis 
myblup$pars$kinship_file=paste0(data_path,"/G_Ainv_col_three.txt")
myblup$pars$analysis_model="GBLUP_A"
myblup$run(task_type="DMU")

#ssGBLUP-A analysis 
myblup$pars$kinship_file=paste0(data_path,"/",c("pedigree.txt","G_A_col_three.txt"))
myblup$pars$analysis_model="SSBLUP_A"
myblup$run(task_type="DMU")

Multiple traits BLUP model

library(blupSUP)
library(blupADC)
data_path=system.file("extdata", package = "blupSUP")
#prepare parameters the BLUP analysis
blup_par=BLUPpar$new(
phe_file=paste0(data_path,"/phenotype.txt"),                                               phe_colnames=c("Id","Mean","Sex","Herd_Year_Season","Litter","Trait1","Trait2","Age"),
dmu_integer_n=5, #number of integer variable in the phenotype file (required by DMU)
analysis_model="PBLUP_A",
kinship_file=paste0(data_path,"/pedigree.txt"))

#define the effect of Trait1 for BLUP analysis
myblup=BLUP$new(fixed=Trait1~Sex+Herd_Year_Season, #fixed effect
               random=~Id+Litter, #random effect 
               pars=blup_par) #blup_par has been defined previously
myblup$run(output_path="D:/test",task_type="DMU") # default using DMU 
myblup

myblup

## <BLUP::DMU> 
## <ADCmodel::1 trait> 
## ----Trait:Trait1~Sex(fixed)+Herd_Year_Season(fixed)+Id(random)+Litter(random)

#Multiple traits analysis
#add model effect of Trait2 in BLUP analysis
myblup$add_model(fixed=Trait2~Herd_Year_Season,random=~Id,covariate=~Age)
myblup

## <BLUP::DMU> 
## <ADCmodel::2 traits> 
## ----Trait:Trait1~Sex(fixed)+Herd_Year_Season(fixed)+Id(random)+Litter(random)
## ----Trait:Trait2~Herd_Year_Season(fixed)+Id(random)+Age(covariate)

myblup$run() #Multiple traits analysis                 

#remove Trait2 in BLUP analysis 
#myblup$rm_model("Trait2") 

#output of BLUP analysis
ebv=myblup$ebv         #breeding value
h2=myblup$h2   #heritability and its SE
genR2=myblup$r2 #genetic correlation and its SE 
user_defined_h2=myblup$cal_se(~Id1/(Id1+Litter1+Res1)) #user defined expression of calculating h2 based on delta method
user_defined_genR2=myblup$cal_se(~Id1_2/sqrt(Id1 * Id2)) #user defined expression of calculating genR2 based on delta method
#random effect names from multiple traits model were defined in: rownames(myblup$vars_se$vars_mat)         

Cross Validation

library(blupSUP)
library(blupADC)
data_path=system.file("extdata", package = "blupSUP")

#prepare parameters the BLUP analysis
blup_par=BLUPpar$new(
phe_file=paste0(data_path,"/phenotype.txt"),                                               phe_colnames=c("Id","Mean","Sex","Herd_Year_Season","Litter","Trait1","Trait2","Age"),
dmu_integer_n=5, #number of integer variable in the phenotype file (required by DMU)
analysis_model="PBLUP_A",
kinship_file=paste0(data_path,"/pedigree.txt"))

#define the effect of Trait1 for BLUP analysis
myblup=BLUP$new(fixed=Trait1~Sex+Herd_Year_Season, #fixed effect
                random=~Id+Litter, #random effect 
                pars=blup_par) #blup_par has been defined previously


#define the training and test datasets by holdout 
myCV=CVpredi$new(myblup,  #has been defined previously 
                 partion_type="holdout", #the type of data partioning
                 ratio=0.6) #the ratio of  train_size/(train_size+test_size)

#get the prediction accuracy 
myCV$predict_sets()
accu=myCV$accuracy
bias=myCV$bias
train_set=myCV$datasets$train_sets[[1]] #train set was used in analysis
test_set=myCV$datasets$test_sets[[1]]   #train set was used in analysis 
ebv_set=myCV$datasets$ebv_sets[[1]]   #train set was used in analysis 

#define the training and test datasets by kfold 
myCV=CVpredi$new(myblup,  #has been defined previously 
                 partion_type="kfold", #the type of data partioning
                 ratio=0.6,            #the ratio of  train_size/(train_size+test_size)
                 kfold=5,              #number of fold  
                 seed=9894)            #random seed number for partioning data,default is 9894     

#get the prediction accuracy 
myCV$predict_sets() #perform analysis in all folds 
#myCV$predict_sets(2) # only perform analysis in the first fold 
accu=myCV$accuracy
bias=myCV$bias
train_set2=myCV$datasets$train_sets[[2]] #train set was used in analysis of fold2
test_set2=myCV$datasets$test_sets[[2]]   #train set was used in analysis  of fold2
ebv_set2=myCV$datasets$ebv_sets[[2]]   #train set was used in analysis  of fold2