### Load packages

```{r}
remove(list = ls())

library(rrBLUP)
library(dplyr)
library(qqman)
library(CMplot)
library(ggplot2)
library(reshape)
library(grDevices)
library(RColorBrewer)
```


Genome Wide Association

```{r}

M <- matrix(rep(0,500*1100),1100,500) # matrix of genotypes dimensions

for (i in 1:500) { # fill the genotype matrix
M[,i] <- ifelse(runif(1100)<0.5,-1,1)
}


geno <- data.frame(SNP=1:1100, # create the data frame with chromosomes
                   Chromosome=rep(1:11,each=100),
                   Position=rep(1:100, times=11),
                   M,check.names=FALSE)

geno[1090:1100,1:5]

QTL <- 55 * (1:20) #pick QTL
length(QTL)
u <- rep(0,1100) 
u[QTL] <- c(20,	40,	60,	80,	100, #marker effects
            120,	140,	160,	180, 
            200,	220,	240,	260, 
            280,	300,	320, 340,	
            360, 380,400)
g <- as.vector(crossprod(M,u)) # get the cross product between Markers and marker effects
h2 <- 0.6 # heritability
y <- g + rnorm(500,mean=4000,sd=sqrt((1-h2)/h2*var(g))) # phenotype in a random distribution on 500 samples  with  the mean and sd desired

hist(y, col = "cornflowerblue", breaks = 30, 
     main="Yield", xlab = "kg/ha")

pheno <- data.frame(line=1:500,y=y) # get the data frame input ofr the GWA 

```


### Genomic Selection Random & Fixed variables

```{r warning=F, warning=F}

phenoGS <- scale(pheno[,2]) # sacale data (subtract the mean and divide by sd). This could help to increase the prediction accuracy

dim(M) # M is the markers Matrix
dim(phenoGS) # number of rows must be equal with M

markers <- t(M) # Transpose of M
dim(markers) # 500 Inv and 1100 markers
dim(phenoGS)

QTL # QTL positions
u[QTL] # QTL effects

# Objects with QTL positions in the genome to use them as fixed variables
indexQTL1 = 1100 
indexQTL2 = c(1045, 1100)
indexQTL3 = c(990, 1045, 1100)
indexQTL4 = c(935 , 990, 1045, 1100)
indexQTL5 = c(880, 935 , 990, 1045, 1100)
indexQTL6 = c(825 ,880, 935, 990, 1045, 1100)
indexQTL7 = c(770, 825 ,880, 935, 990, 1045, 1100)
indexQTL8 = c(715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL9 = c(660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL10 = c(605, 660,715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL11 = c(550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL12 = c(495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL13 = c(440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL14 = c(385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL15 = c(330, 385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL16 = c(275, 330, 385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL17 = c(220, 275, 330, 385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990,  1045, 1100)
indexQTL18 = c(165, 220, 275, 330,385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL19 = c(110, 165, 220, 275, 330, 385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)
indexQTL20 = c(55, 110, 165, 220, 275, 330, 385, 440, 495, 550, 605, 660, 715, 770, 825 ,880, 935, 990, 1045, 1100)

# Exclude fixed variables for each data set
markersNoQTL1 = markers[,-indexQTL1] 
markersNoQTL2 = markers[,-indexQTL2]
markersNoQTL3 = markers[,-indexQTL3]
markersNoQTL4 = markers[,-indexQTL4]
markersNoQTL5 = markers[,-indexQTL5]
markersNoQTL6 = markers[,-indexQTL6]
markersNoQTL7 = markers[,-indexQTL7]
markersNoQTL8 = markers[,-indexQTL8]
markersNoQTL9 = markers[,-indexQTL9]
markersNoQTL10 = markers[,-indexQTL10]
markersNoQTL11 = markers[,-indexQTL11]
markersNoQTL12 = markers[,-indexQTL12]
markersNoQTL13 = markers[,-indexQTL13]
markersNoQTL14 = markers[,-indexQTL14]
markersNoQTL15 = markers[,-indexQTL15]
markersNoQTL16 = markers[,-indexQTL16]
markersNoQTL17 = markers[,-indexQTL17]
markersNoQTL18 = markers[,-indexQTL18]
markersNoQTL19 = markers[,-indexQTL19]
markersNoQTL20 = markers[,-indexQTL20]

# Now we need to create the matrix with the fixed variables for each data set
bQTL1 = as.matrix(markers[,indexQTL1]) 
bQTL2 = as.matrix(markers[,indexQTL2])
bQTL3 = as.matrix(markers[,indexQTL3])
bQTL4 = as.matrix(markers[,indexQTL4])
bQTL5 = as.matrix(markers[,indexQTL5])
bQTL6 = as.matrix(markers[,indexQTL6])
bQTL7 = as.matrix(markers[,indexQTL7])
bQTL8 = as.matrix(markers[,indexQTL8])
bQTL9 = as.matrix(markers[,indexQTL9])
bQTL10 = as.matrix(markers[,indexQTL10])
bQTL11 = as.matrix(markers[,indexQTL11])
bQTL12 = as.matrix(markers[,indexQTL12])
bQTL13 = as.matrix(markers[,indexQTL13])
bQTL14 = as.matrix(markers[,indexQTL14])
bQTL15 = as.matrix(markers[,indexQTL15])
bQTL16 = as.matrix(markers[,indexQTL16])
bQTL17 = as.matrix(markers[,indexQTL17])
bQTL18 = as.matrix(markers[,indexQTL18])
bQTL19 = as.matrix(markers[,indexQTL19])
bQTL20 = as.matrix(markers[,indexQTL20])
head(bQTL20)

cycles = 10 # number of cycles 
accuracy = matrix(nrow = cycles, ncol=21) # empty matrix to store results

for(r in 1:cycles) 
  { # open loop to run GP model
train = as.matrix(sample(1:500, 250)) # select 250 lines to train the model
test = setdiff(1:500, train) #select remainder samples to validate the model

yTRN  = phenoGS[train,] # select 250 lines to train the model

# extract samples in the training population for each data set of markers
xTRN = markers[train,]
x1TRN = markersNoQTL1[train,]
x2TRN = markersNoQTL2[train,]
x3TRN = markersNoQTL3[train,]
x4TRN = markersNoQTL4[train,]
x5TRN = markersNoQTL5[train,]
x6TRN = markersNoQTL6[train,]
x7TRN = markersNoQTL7[train,]
x8TRN = markersNoQTL8[train,]
x9TRN = markersNoQTL9[train,]
x10TRN = markersNoQTL10[train,]
x11TRN = markersNoQTL11[train,]
x12TRN = markersNoQTL12[train,]
x13TRN = markersNoQTL13[train,]
x14TRN = markersNoQTL14[train,]
x15TRN = markersNoQTL15[train,]
x16TRN = markersNoQTL16[train,]
x17TRN = markersNoQTL17[train,]
x18TRN = markersNoQTL18[train,]
x19TRN = markersNoQTL19[train,]
x20TRN = markersNoQTL20[train,]

# extract samples in the training population for each data set of fixed variables
bTRN1 = bQTL1[train, ,drop=F]
bTRN2 = bQTL2[train, ,drop=F]
bTRN3 = bQTL3[train, ,drop=F]
bTRN4 = bQTL4[train, ,drop=F]
bTRN5 = bQTL5[train, ,drop=F]
bTRN6 = bQTL6[train, ,drop=F]
bTRN7 = bQTL7[train, ,drop=F]
bTRN8 = bQTL8[train, ,drop=F]
bTRN9 = bQTL9[train, ,drop=F]
bTRN10 = bQTL10[train, ,drop=F]
bTRN11 = bQTL11[train, ,drop=F]
bTRN12 = bQTL12[train, ,drop=F]
bTRN13 = bQTL13[train, ,drop=F]
bTRN14 = bQTL14[train, ,drop=F]
bTRN15 = bQTL15[train, ,drop=F]
bTRN16 = bQTL16[train, ,drop=F]
bTRN17 = bQTL17[train, ,drop=F]
bTRN18 = bQTL18[train, ,drop=F]
bTRN19 = bQTL19[train, ,drop=F]
bTRN20 = bQTL20[train, ,drop=F]

# extract samples in the testing population for each data set of fixed variables

bTST1 = bQTL1[test, ,drop=F]
bTST2 = bQTL2[test, ,drop=F]
bTST3 = bQTL3[test, ,drop=F]
bTST4 = bQTL4[test, ,drop=F]
bTST5 = bQTL5[test, ,drop=F]
bTST6 = bQTL6[test, ,drop=F]
bTST7 = bQTL7[test, ,drop=F]
bTST8 = bQTL8[test, ,drop=F]
bTST9 = bQTL9[test, ,drop=F]
bTST10 = bQTL10[test, ,drop=F]
bTST11 = bQTL11[test, ,drop=F]
bTST12 = bQTL12[test, ,drop=F]
bTST13 = bQTL13[test, ,drop=F]
bTST14 = bQTL14[test, ,drop=F]
bTST15 = bQTL15[test, ,drop=F]
bTST16 = bQTL16[test, ,drop=F]
bTST17 = bQTL17[test, ,drop=F]
bTST18 = bQTL18[test, ,drop=F]
bTST19 = bQTL19[test, ,drop=F]
bTST20= bQTL20[test, ,drop=F]


yTST  = phenoGS[test,] # select remainder samples to validate the model in the testing population

# extract samples in the testing population for each data set of markers
xTST = markers[test,]
x1TST = markersNoQTL1[test,]
x2TST = markersNoQTL2[test,]
x3TST = markersNoQTL3[test,]
x4TST = markersNoQTL4[test,]
x5TST = markersNoQTL5[test,]
x6TST = markersNoQTL6[test,]
x7TST = markersNoQTL7[test,]
x8TST = markersNoQTL8[test,]
x9TST = markersNoQTL9[test,]
x10TST = markersNoQTL10[test,]
x11TST = markersNoQTL11[test,]
x12TST = markersNoQTL12[test,]
x13TST = markersNoQTL13[test,]
x14TST = markersNoQTL14[test,]
x15TST = markersNoQTL15[test,]
x16TST = markersNoQTL16[test,]
x17TST = markersNoQTL17[test,]
x18TST = markersNoQTL18[test,]
x19TST = markersNoQTL19[test,]
x20TST = markersNoQTL20[test,]

yNA = phenoGS # phenotype 
yNA[test] <- NA # data in testing pop will be "NA"

# rrBLUP package
fmGS0 = mixed.solve(yTRN, # only random effects
                    Z=xTRN)

fmGS1 = mixed.solve(yTRN, # QTL fixed  + random (1 QTL)
                    X=bTRN1, 
                    Z=x1TRN)

fmGS2 = mixed.solve(yTRN, # QTL fixed  + random (2 QTL)
                    X=bTRN2, 
                    Z=x2TRN)

fmGS3 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN3, 
                    Z=x3TRN)

fmGS4 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN4, 
                    Z=x4TRN)

fmGS5 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN5, 
                    Z=x5TRN)

fmGS6 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN6, 
                    Z=x6TRN)
fmGS7 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN7, 
                    Z=x7TRN)
fmGS8 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN8, 
                    Z=x8TRN)
fmGS9 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN9, 
                    Z=x9TRN)
fmGS10 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN10, 
                    Z=x10TRN)
fmGS11 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN11, 
                    Z=x11TRN)
fmGS12 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN12, 
                    Z=x12TRN)
fmGS13 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN13, 
                    Z=x13TRN)
fmGS14 = mixed.solve(yTRN, # QTL fixed  + random
                    X=bTRN14, 
                    Z=x14TRN)
fmGS15 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN15, 
                    Z=x15TRN)
fmGS16 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN16, 
                    Z=x16TRN)
fmGS17 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN17, 
                    Z=x17TRN)
fmGS18 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN18, 
                    Z=x18TRN)
fmGS19 = mixed.solve(yTRN, # QTL fixed  + random 
                    X=bTRN19, 
                    Z=x19TRN)
fmGS20 = mixed.solve(yTRN, # QTL fixed  + random (20 QTL)
                    X=bTRN20, 
                    Z=x20TRN)

yHat0 = xTST %*% fmGS0$u  # Predicted genetic value

# Predicted value with fixed variables 
yHat1 = (bTST1 %*% fmGS1$beta) + (x1TST %*% fmGS1$u) 
yHat2 = (bTST2 %*% fmGS2$beta) + (x2TST %*% fmGS2$u)
yHat3 = (bTST3 %*% fmGS3$beta) + (x3TST %*% fmGS3$u)
yHat4 = (bTST4 %*% fmGS4$beta) + (x4TST %*% fmGS4$u)
yHat5 = (bTST5 %*% fmGS5$beta) + (x5TST %*% fmGS5$u)
yHat6 = (bTST6 %*% fmGS6$beta) + (x6TST %*% fmGS6$u)
yHat7 = (bTST7 %*% fmGS7$beta) + (x7TST %*% fmGS7$u)
yHat8 = (bTST8 %*% fmGS8$beta) + (x8TST %*% fmGS8$u)
yHat9 = (bTST9 %*% fmGS9$beta) + (x9TST %*% fmGS9$u)
yHat10 = (bTST10 %*% fmGS10$beta) + (x10TST %*% fmGS10$u)
yHat11 = (bTST11 %*% fmGS11$beta) + (x11TST %*% fmGS11$u)
yHat12 = (bTST12 %*% fmGS12$beta) + (x12TST %*% fmGS12$u)
yHat13 = (bTST13 %*% fmGS13$beta) + (x13TST %*% fmGS13$u)
yHat14 = (bTST14 %*% fmGS14$beta) + (x14TST %*% fmGS14$u)
yHat15 = (bTST15 %*% fmGS15$beta) + (x15TST %*% fmGS15$u)
yHat16 = (bTST16 %*% fmGS16$beta) + (x16TST %*% fmGS16$u)
yHat17 = (bTST17 %*% fmGS17$beta) + (x17TST %*% fmGS17$u)
yHat18 = (bTST18 %*% fmGS18$beta) + (x18TST %*% fmGS18$u)
yHat19 = (bTST19 %*% fmGS19$beta) + (x19TST %*% fmGS19$u)
yHat20 = (bTST20 %*% fmGS20$beta) + (x20TST %*% fmGS20$u)

# store the correlation of each cycle
accuracy[r,1] = cor(yHat0, yTST, use = "complete" ) 
accuracy[r,2] = cor(yHat1, yTST, use = "complete" )
accuracy[r,3] = cor(yHat2, yTST, use = "complete" )
accuracy[r,4] = cor(yHat3, yTST, use = "complete" )
accuracy[r,5] = cor(yHat4, yTST, use = "complete" )
accuracy[r,6] = cor(yHat5, yTST, use = "complete" )
accuracy[r,7] = cor(yHat6, yTST, use = "complete" )
accuracy[r,8] = cor(yHat7, yTST, use = "complete" )
accuracy[r,9] = cor(yHat8, yTST, use = "complete" )
accuracy[r,10] = cor(yHat9, yTST, use = "complete" )
accuracy[r,11] = cor(yHat10, yTST, use = "complete" )
accuracy[r,12] = cor(yHat11, yTST, use = "complete" )
accuracy[r,13] = cor(yHat12, yTST, use = "complete" )
accuracy[r,14] = cor(yHat13, yTST, use = "complete" )
accuracy[r,15] = cor(yHat14, yTST, use = "complete" )
accuracy[r,16] = cor(yHat15, yTST, use = "complete" )
accuracy[r,17] = cor(yHat16, yTST, use = "complete" )
accuracy[r,18] = cor(yHat17, yTST, use = "complete" )
accuracy[r,19] = cor(yHat18, yTST, use = "complete" )
accuracy[r,20] = cor(yHat19, yTST, use = "complete" )
accuracy[r,21] = cor(yHat20, yTST, use = "complete" )
}

apply(accuracy,2,mean) # mean accuracy of cycles

predGS0 <- data.frame(yTST, yHat0) # create the data frame input for GP plot

p_GS0 <- ggplot(predGS0, 
               aes(x=yHat0, y=yTST)) +
  geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) 

p_GS0


colnames(accuracy) = c("0", "1","2", "3", "4", "5",
                       "6","7","8","9","10",
                    "11","12","13","14","15",
                    "16","17","18","19","20")

boxplot(accuracy, col=c(2:23), 
        ylab= "Accuracy",
        xlab= "Fixed QTL", cex.axis=0.8)

```



GS only with QTLs

```{r }

QTL # QTL regions
markers_qtl = markers[,QTL] # select QTL markers
dim(markers_qtl)

traits=1
cycles=10
accuracy2= matrix(nrow = cycles, ncol=traits)

for(r in 1:cycles) {
train = as.matrix(sample(1:500, 350))
test = setdiff(1:500, train)

pheno_train_gq  = phenoGS[train,]
marker_train_gq = markers_qtl[train,]

pheno_valid_gq  = phenoGS[test,]
marker_valid_gq = markers_qtl[test,]

yield_answer_gq = mixed.solve(pheno_train_gq, 
                           Z=marker_train_gq,
                           K=NULL,
                           SE= F,
                           return.Hinv = F)

YLD_gq = yield_answer_gq$u # marker effects
e_q = as.matrix(YLD_gq)

pred_yield_valid_gq = marker_valid_gq %*% e_q
pred_yield_gq = as.vector(pred_yield_valid_gq) + as.vector(yield_answer_gq$beta)
accuracy2[r,1] = cor(pred_yield_gq, 
                    pheno_valid_gq, 
                   use = "complete" )
}

mean(accuracy2)

pred2 <- data.frame(pheno_valid_gq, pred_yield_gq)

p_GS3 <- ggplot(pred2, aes(x=pheno_valid_gq,
                          y=pred_yield_gq)) +
  geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) 

cor(pheno_valid_gq, pred_yield_gq)
p_GS3

```


