# NOT RUN {
####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples using
#### command + shift + C |OR| control + shift + C
####=========================================####
data(RICE)
W <- RICE$RiceGenoN; W[1:5,1:5]; dim(W)
#### lets use only 5000 random markers across the genome
#### for demonstration purposes
ww <- sample(1:dim(W)[2],5000)
W <- W[,ww]; dim(W)
Y <- RICE$RicePheno; Y[1:5,1:5]; dim(Y)
(vp.str <- RICE$vp.str[1:25]) # validation population with structure
bases <- seq(50,200,25) # TP sizes required
####=========================================####
####=========================================####
#### a) Random sample based
####=========================================####
####=========================================####
#### rest of the plants not in the validation population (VP)
rest <- (1:dim(W)[1])[-which(rownames(W) %in% vp.str)]
#### select the best TP for the VP, provide:
#### markers, names of VP, size of the TP desired, method
#tp.plant.strR <- TP.prep(markers=W, vp.names = vp.str,
# tp.size = bases, method="random")
####=========================================####
#### b) Genetic simmilarity-dissimilarity based
####=========================================####
#tp.plant.strG1 <- TP.prep(markers=W, vp.names = vp.str,
# tp.size = bases, method="sim-dissim")
#tp.plant.strG2 <- TP.prep(markers=W, vp.names = vp.str,
# tp.size = bases, method="sim")
####=========================================####
####=========================================####
#### c) QTL based
#### we assume real situation that VP
#### wouldn't be phenotyped for GWAS
####=========================================####
####=========================================####
#K1 <- A.mat(W)
#Z1 <- diag(dim(K1)[1])
#image(K1)
#ETA1 <- list(list(Z=Z1,K=K1))
#YY <- Y; YY[vp.str,] <- NA
####=========================================####
#### be patient, the GWAS is for >36,000 SNPs
####=========================================####
#ans.GWAS <- mmer(Y=YY$Protein.content, Z=ETA1, W=W)
#(h2 <- sqrt(ans.GWAS$var.comp[1,]/sum(ans.GWAS$var.comp)))
####=========================================####
#### get marker matrix for the most significant QTLs
#### and do the TP selection using only those markers
####=========================================####
#X <- hits(ans.GWAS, nmar = 25); head(X); dim(X)
#tp.plant.strG3 <- TP.prep(markers=X, vp.names = vp.str,
# tp.size = bases, method="sim-dissim")
####=========================================####
####=========================================####
#### COMPARE METHODS with 50 plants
####=========================================####
####=========================================####
#metho <- list(a=tp.plant.strR[[1]], b=tp.plant.strG1[[1]],
# c=tp.plant.strG2[[1]], d=tp.plant.strG3[[1]])
#resos <- numeric()
#for(i in 1:4){ # for each method
# tpgi <- metho[[i]] # training pop genetic-based
# Wpgi <- W[c(vp.str,tpgi),] # provisional markers based on genes
# Ypgi <- Y[c(vp.str,tpgi),] # provisional phenos based on genes
# Xpgi <- X[c(vp.str,tpgi),] # provisional bag-mat based on genes
# Kpgi <- A.mat(Wpgi)
# Zpgi <- diag(dim(Kpgi)[1])
# Ypgi[vp.str,] <- NA
# yt <- Ypgi$Protein.content
# ETAg <- list(list(Z=Zpgi, K=Kpgi))
# mixg <- mmer(Y=yt, X=Xpgi, Z=ETAg, method="EMMA") #X=Xpgi,
# resos[i] <- cor(Y[vp.str,"Protein.content"],
# fitted(mixg)[(1:length(vp.str)),],
# use="complete")
#}
#resos
### the random method needs more iterations to have a real idea what
### would be the predictive ability of using plants at random,
### this is just to give an idea to users.
# }
Run the code above in your browser using DataLab