R/GCPC.R

   1 ################################################################################
   2 # Genomic prediction of cross performance for YamBase
   3 ################################################################################
   4
   5 # There are ten main steps to this protocol:
   6 # 1. Load the software needed.
   7 # 2. Declare user-supplied variables.
   8 # 3. Read in the genotype data and convert to numeric allele counts.
   9 # 4. Get the genetic predictors needed.
  10 # 5. Process the phenotypic data.
  11 # 6. Fit the mixed models in sommer.
  12 # 7. Backsolve from individual estimates to marker effect estimates / GBLUP -> RR-BLUP
  13 # 8. Weight the marker effects and add them together to form an index of merit.
  14 # 9. Predict the crosses.
  15 # 10. Format the information needed for output.
  16
  17
  18 #Get Arguments
  19 args = commandArgs(trailingOnly=TRUE)
  20 if (length(args) < 3) {
  21   stop('Two or more arguments are required.')
  22 }
  23 phenotypeFile = args[1]
  24 genotypeFile = args[2]
  25 traits = args[3]
  26 weights = args[4]
  27 userSexes = args[5]
  28 userFixed = args[6]
  29 userRandom = args[7]
  30
  31
  32 write(paste('PLANT SEX CVTERM: ', userSexes), stderr())
  33
  34 ################################################################################
  35 # 1. Load software needed
  36 ################################################################################
  37
  38 library(sommer)
  39 library(AGHmatrix)
  40 library(VariantAnnotation) # Bioconductor package
  41 library(tools)
  42 Rcpp::sourceCpp("/home/production/cxgn/QuantGenResources/CalcCrossMeans.cpp") # this is called CalcCrossMean.cpp on Github
  43
  44
  45
  46
  47
  48
  49 ################################################################################
  50 # 2. Declare user-supplied variables
  51 ################################################################################
  52
  53 # a. Define path with internal YamBase instructions such that the object 'userGeno'
  54 #    is defined as a VCF file of genotypes.
  55
  56 #userGeno <- path
  57
  58
  59 # b. Define path2 with internal YamBase instructions such that the object 'userPheno'
  60 #    is defined as the phenotype file.
  61
  62 #userPheno <- path2
  63 #write(paste("READING PHENOTYPEFILE: ",phenotypeFile), stderr())
  64 userPheno <- read.delim(phenotypeFile, header = TRUE, sep="\t", fill=TRUE) #testing only
  65 #write(colnames(userPheno), stderr())
  66 #write(summary(userPheno), stderr())
  67 ##userPheno <- userPheno[userPheno$Trial == "SCG", ] #testing only-- needs to replaced with 2-stage
  68
  69 #write("DONE WITH PHENOTYPEFILE"), stderr())
  70
  71 # c. The user should be able to select their fixed variables from a menu
  72 #    of the column names of the userPheno object. The possible interaction terms
  73 #    also need to be shown somehow. Then, those strings should be passed
  74 #    to this vector, 'userFixed'. Please set userFixed to NA if no fixed effects
  75 #    besides f are requested.
  76 #    f is automatically included as a fixed effect- a note to the user would be good.
  77
  78 #userFixed <- c()
  79 #userFixed <- c("studyYear") # for testing only
  80 userFixed <-unlist(strsplit(userFixed, split=",", fixed=T))
  81
  82
  83 # d. The user should be able to select their random variables from a menu
  84 #    of the column names of the userPheno object. The possible interaction terms
  85 #    also need to be shown somehow. Then, those strings should be passed
  86 #    to this vector, 'userRandom'.
  87
  88 #userRandom <- c()
  89 #userRandom <- "blockNumber" # for testing only
  90 userRandom <-unlist(strsplit(userRandom, split=",", fixed=T))
  91
  92 # e. The user should be able to indicate which of the userPheno column names
  93 #    represents individual genotypes identically as they are represented in the VCF
  94 #    column names. No check to ensure matching at this stage. This single string
  95 #    should be passed to this vector, userID.
  96
  97 #userID <- c()
  98 userID <- "germplasmName" # for testing only
  99
 100
 101 # f. The user must indicate the ploidy level of their organism, and the integer
 102 #    provided should be passed to the vector 'userPloidy'. CalcCrossMeans.cpp
 103 #    currently supports ploidy = {2, 4, 6}. Ideally, the user could select
 104 #    their ploidy from a drop-down to avoid errors here, and there would be a note
 105 #    that other ploidies are not currently supported. If not, a possible error is
 106 #    provided.
 107
 108 userPloidy <- c()
 109 userPloidy <- 2 # for testing only
 110
 111 # if(userPloidy %in% c(2, 4, 6) != TRUE){
 112 #   stop("Only ploidies of 2, 4, and 6 are supported currently. \n
 113 #        Please confirm your ploidy level is supported.")
 114 # }
 115
 116
 117 # g. The user should be able to select their response variables from a drop-down menu
 118 #    of the column names of the userPheno object. Then, those strings should be passed
 119 #    to this vector, 'userResponse'.
 120
 121 write(paste("TRAIT STRING:", traits), stderr())
 122 #userResponse <- c()
 123 #userResponse <- c("YIELD", "DMC", "OXBI") # for testing only
 124 userResponse <- unlist(strsplit(traits, split=",", fixed=T))
 125
 126 write(paste("USER RESPONSE", userResponse),stderr())
 127 write(paste("first element: ", userResponse[1]), stderr())
 128 # h. The user must indicate weights for each response. The order of the vector
 129 #    of response weights must match the order of the responses in userResponse.
 130
 131 #userWeights <- c()
 132 #userWeights <- c(1, 0.8, 0.2) # for YIELD, DMC, and OXBI respectively; for testing only
 133 userWeights  <- as.numeric(unlist(strsplit(weights, split=",", fixed=T)))
 134
 135 write(paste("WEIGHTS", userWeights), stderr())
 136
 137 # i. The user can indicate the number of crosses they wish to output.
 138 #    The maximum possible is a full diallel.
 139
 140 #userNCrosses <- c()
 141 userNCrosses <- 100 # for testing only
 142
 143
 144 # j. The user can (optionally) input the individuals' sexes and indicate the column
 145 #    name of the userPheno object which corresponds to sex. The column name
 146 #   string should be passed to the 'userSexes' object. If the user does not wish
 147 #   to remove crosses with incompatible sexes (e.g. because the information is not available),
 148 #   then userSexes should be set to NA.
 149
 150
 151 #userSexes <- c()
 152 #userSexes <- "Sex" # for testing only
 153 #userPheno$Sex <- sample(c("M", "F"), size = nrow(userPheno), replace = TRUE, prob = c(0.7, 0.3)) # for testing only
 154 # Please note that for the test above, sex is sampled randomly for each entry, so the same accession can have
 155 # different sexes. This does not matter for the code or testing.
 156
 157 ################################################################################
 158 # 3. Read in the genotype data and convert to numeric allele counts.
 159 ################################################################################
 160
 161 # a. The VCF file object 'userGeno' needs to be converted to a numeric matrix
 162 #    of allele counts in whic:
 163 #    Rownames represent the individual genotype IDs
 164 #    Colnames represent the site IDs
 165 #    A cell within a given row and column represents the row individual's
 166 #    genotype at the site in the column.
 167
 168 #   The individual's genotype should be an integer from 0... ploidy to represent
 169 #   counts of the alternate allele at the site. Diploid example:
 170 #    0 = homozygous reference
 171 #    1 = heterozygous
 172 #    2 = homozygous alternate
 173
 174 #    The genotypes must not contain monomorphic or non-biallelic sites.
 175 #    Users need to pre-process their VCF to remove these (e.g. in TASSEL or R)
 176 #    I can put an error message into this script if a user tries to input
 177 #    monomorphic or biallelic sites which could be communicated through the GUI.
 178 #    It's also possible to filter them here.
 179
 180 if (file_ext(genotypeFile) == 'vcf') {
 181    write(paste("READING VARIANT FILE ", genotypeFile), stderr())
 182    #  Import VCF with VariantAnnotation package and extract matrix of dosages
 183    myVCF <- readVcf(genotypeFile)
 184    #G <- t(geno(myVCF)$DS) # Individual in row, genotype in column
 185    mat <- genotypeToSnpMatrix(myVCF)
 186    #G <- t(geno(myVCF)$DS) # Individual in row, genotype in column
 187    G <- as(mat$genotypes, "numeric")
 188    G <- G[,colSums(is.na(G))<nrow(G)]
 189
 190    #   TEST temporarily import the genotypes via HapMap:
 191    #source("R/hapMap2numeric.R") # replace and delete
 192    #G <- hapMap2numeric(genotypeFile) # replace and delete
 193
 194 } else {
 195    write(paste("READING DOSAGE FILE ", genotypeFile), stderr())
 196      GF <- read.delim(genotypeFile)
 197      GT <- as.matrix(GF)
 198      G <- G[-1,]
 199      G <- t(G)
 200
 201
 202 }
 203
 204 write("G Matrix start --------", stderr())
 205 write(G[1:5, 1:5], stderr())
 206 write("G Matrix end =========", stderr())
 207
 208
 209
 210
 211 ################################################################################
 212 # 4. Get the genetic predictors needed.
 213 ################################################################################
 214
 215 write("GENETIC PREDICTIONS...", stderr())
 216 # 4a. Get the inbreeding coefficent, f, as described by Xiang et al., 2016
 217 # The following constructs f as the average heterozygosity of the individual
 218 # The coefficient of f estimated later then needs to be divided by the number of markers
 219 # in the matrix D before adding it to the estimated dominance marker effects
 220 # One unit of change in f represents changing all loci from homozygous to heterozygous
 221
 222 ###GC <- G - (userPloidy/2) #this centers G
 223 GC <- G * (userPloidy - G) * (2 / userPloidy)^2 # center at G
 224 f <- rowSums(GC, na.rm = TRUE) / apply(GC, 1, function(x) sum(!is.na(x)))
 225
 226 # Another alternate way to construct f is the total number of heterozygous loci in the individual
 227 # The coefficient of this construction of f does not need to be divided by the number of markers
 228 # It is simply added to each marker dominance effect
 229 # The coefficient of this construction of f represents the average dominance effect of a marker
 230 # One unit of change in f represents changing one locus from homozygous to heterozygous
 231 # f <- rowSums(D, na.rm = TRUE)
 232
 233
 234 write("DISTANCE MATRIX...", stderr())
 235 # 4b. Get the additive and dominance relationship matrices following Batista et al., 2021
 236 # https://doi.org/10.1007/s00122-021-03994-w
 237
 238 # Additive: this gives a different result than AGHmatrix VanRaden's Gmatrix
 239 # AGHmatrix: Weights are implemented for "VanRaden" method as described in Liu (2020)?
 240 allele_freq = colSums(G) / (userPloidy * nrow(G))
 241 W = t(G) - userPloidy * allele_freq
 242 WWt = crossprod(W)
 243 denom = sum(userPloidy * allele_freq * (1 - allele_freq))
 244 A = WWt / denom
 245
 246 # Check with paper equation:
 247 #w <- G - (userPloidy/2)
 248 #num <- w %*% t(w)
 249 #denom = sum(userPloidy * allele_freq * (1 - allele_freq))
 250 #A2 <- num/denom
 251 #table(A == A2)
 252 #cor(as.vector(A), as.vector(A2)) # 0.9996...
 253
 254
 255 # Dominance or digenic dominance
 256 if(userPloidy == 2){
 257   D <- Gmatrix(G, method = "Su", ploidy = userPloidy, missingValue = NA)
 258 }
 259
 260 if(userPloidy > 2){
 261 # Digenic dominance
 262 C_matrix = matrix(length(combn(userPloidy, 2)) / 2,
 263                   nrow = nrow(t(G)),
 264                   ncol = ncol(t(G)))
 265
 266 Ploidy_matrix = matrix(userPloidy,
 267                        nrow = nrow(t(G)),
 268                        ncol = ncol(t(G)))
 269
 270 Q = (allele_freq^2 * C_matrix) -
 271   (Ploidy_matrix - 1) * allele_freq * t(G) +
 272   0.5 * t(G) * (t(G) - 1)
 273
 274 Dnum = crossprod(Q)
 275 denomDom = sum(C_matrix[ ,1] * allele_freq^2 * (1 - allele_freq)^2)
 276 D = Dnum/denomDom
 277 }
 278
 279
 280
 281
 282
 283
 284 ################################################################################
 285 # 5. Process the phenotypic data.
 286 ################################################################################
 287
 288 #write(summary(userPheno), stderr())
 289
 290 # a. Paste f into the phenotype dataframe
 291 write("processing phenotypic data...", stderr())
 292 userPheno$f <- f[as.character(userPheno[ , userID])]
 293
 294 #write(summary(userPheno), stderr())
 295
 296 # b. Scale the response variables.
 297 write("processing phenotypic data... scaling...", stderr())
 298 write(paste("USER RESPONSE LENGTH = ", length(userResponse)), stderr())
 299 for(i in 1:length(userResponse)){
 300   write(paste("working on user response ", userResponse[i]), stderr())
 301   userPheno[ , userResponse[i]] <- (userPheno[ , userResponse[i]] - mean(userPheno[ , userResponse[i]], na.rm = TRUE)) / sd(userPheno[ , userResponse[i]], na.rm = TRUE)
 302 }
 303
 304 write(paste("accession count: ", length(userPheno[ , userID])), stderr())
 305 write("processing phenotypic data... adding dominance effects", stderr())
 306 # c. Paste in a second ID column for the dominance effects.
 307
 308 #write(summary(userPheno), stderr())
 309
 310 dominanceEffectCol = paste(userID, "2", sep="")
 311 write(paste("NEW COL NAME: ", dominanceEffectCol), stderr())
 312
 313 write(paste("USER_ID COLUMN: ", userPheno[ , userID]), stderr());
 314
 315
 316 userPheno[ , dominanceEffectCol] <- userPheno[ , userID]
 317
 318 write(paste("USER PHENO userID2 COL", userPheno[ , dominanceEffectCol]), stderr())
 319
 320
 321 # Additional steps could be added here to remove outliers etc.
 322
 323
 324
 325
 326
 327 ################################################################################
 328 # 6. Fit the mixed models in sommer.
 329 ################################################################################
 330
 331 write("Fit mixed model in sommer", stderr());
 332
 333 # 6a. Make a list to save the models.
 334
 335 userModels <- list()
 336
 337 for(i in 1:length(userResponse)){
 338
 339   write(paste("User response: ", userResponse[i]), stderr())
 340   # check if fixed effects besides f are requested, then paste together
 341   # response variable and fixed effects
 342   if(!is.na(userFixed[1])){
 343   fixedEff <- paste(userFixed, collapse = " + ")
 344   fixedEff <- paste(fixedEff, "f", sep = " + ")
 345   fixedArg <- paste(userResponse[i], " ~ ", fixedEff, sep = "")
 346   }
 347   if(is.na(userFixed[1])){
 348     fixedArg <- paste(userResponse[i], " ~ ", "f")
 349   }
 350
 351
 352   # check if random effects besides genotypic additive and dominance effects
 353   # are requested, then paste together the formula
 354
 355   write("Generate formula...", stderr())
 356
 357   if(!is.na(userRandom[1])){
 358     randEff <- paste(userRandom, collapse = " + ")
 359     ID2 <- paste(userID, 2, sep = "")
 360     randEff2 <- paste("~vsr(", userID, ", Gu = A) + vsr(", ID2, ", Gu = D)", sep = "")
 361     randArg <- paste(randEff2, randEff, sep = " + ")
 362   }
 363   if(is.na(userRandom[1])){
 364     ID2 <- paste(userID, 2, sep ="")
 365     randArg <- paste("~vsr(", userID, ", Gu = A) + vsr(", ID2, ", Gu = D)", sep = "")
 366   }
 367
 368   write(paste("Fit mixed GBLUP model...", randArg), stderr())
 369
 370   #  write(paste("USER PHENO:", userPheno), stderr())
 371   #  write(paste("COLNAMES: ", colnames(userPheno)), stderr())
 372   # fit the mixed GBLUP model
 373   myMod <- mmer(fixed = as.formula(fixedArg),
 374                 random = as.formula(randArg),
 375                 rcov = ~units,
 376                 getPEV = FALSE,
 377                 data = userPheno)
 378
 379
 380   # save the fit model
 381
 382   write(paste("I = ", i), stderr());
 383
 384   userModels[[i]] <- myMod
 385 }
 386
 387
 388
 389
 390
 391
 392 ######################################################################################
 393 # 7. Backsolve from individual estimates to marker effect estimates / GBLUP -> RR-BLUP
 394 ######################################################################################
 395
 396 # a. Get the matrices and inverses needed
 397 #    This is not correct for polyploids yet.
 398 A.G <- G - (userPloidy / 2) # this is the additive genotype matrix (coded -1 0 1 for diploids)
 399 D.G <- 1 - abs(A.G)     # this is the dominance genotype matrix (coded 0 1 0 for diploids)
 400
 401
 402 A.T <- A.G %*% t(A.G) ## additive genotype matrix
 403 A.Tinv <- solve(A.T) # inverse; may cause an error sometimes, if so, add a small amount to the diag
 404 A.TTinv <- t(A.G) %*% A.Tinv # M'%*% (M'M)-
 405
 406 D.T <- D.G %*% t(D.G) ## dominance genotype matrix
 407 D.Tinv <- solve(D.T) ## inverse
 408 D.TTinv <- t(D.G) %*% D.Tinv # M'%*% (M'M)-
 409
 410
 411 # b. Loop through and backsolve to marker effects.
 412
 413 write("backsolve marker effects...", stderr())
 414
 415 userAddEff <- list() # save them in order
 416 userDomEff <- list() # save them in order
 417
 418 for(i in 1:length(userModels)){
 419
 420   myMod <- userModels[[i]]
 421
 422   # get the additive and dominance effects out of the sommer list
 423   subMod <- myMod$U
 424   subModA <- subMod[[1]]
 425   subModA <- subModA[[1]]
 426   subModD <- subMod[[2]]
 427   subModD <- subModD[[1]]
 428
 429   # backsolve
 430   addEff <- A.TTinv %*% matrix(subModA[colnames(A.TTinv)], ncol = 1) # these must be reordered to match A.TTinv
 431   domEff <- D.TTinv %*% matrix(subModD[colnames(D.TTinv)], ncol=1)   # these must be reordered to match D.TTinv
 432
 433   # add f coefficient back into the dominance effects
 434   subModf <- myMod$Beta
 435   fCoef <- subModf[subModf$Effect == "f", "Estimate"] # raw f coefficient
 436   fCoefScal <- fCoef / ncol(G) # divides f coefficient by number of markers
 437   dirDomEff <- domEff + fCoefScal
 438
 439   # save
 440   userAddEff[[i]] <- addEff
 441   userDomEff[[i]] <- dirDomEff
 442 }
 443
 444
 445
 446
 447
 448
 449 ################################################################################
 450 # 8. Weight the marker effects and add them together to form an index of merit.
 451 ################################################################################
 452
 453 write("weight marker effects...", stderr())
 454
 455 ai <- 0
 456 di <- 0
 457 for(i in 1:length(userWeights)){
 458   write(paste("USER ADD EFF : ", userAddEff[[i]]), stderr())
 459     write(paste("USER DOM EFF : ", userDomEff[[i]]), stderr())
 460   write(paste("USER WEIGHT : ", userWeights[i]), stderr())
 461
 462
 463        ai <- ai + userAddEff[[i]] *  userWeights[i]
 464
 465  write("DONE WITH ADDITIVE EFFECTS!\n", stderr())
 466   di <- di + userDomEff[[i]] * userWeights[i]
 467    write("DONE WITH DOM EFFECTS!\n", stderr())
 468 }
 469
 470
 471
 472
 473
 474
 475 ################################################################################
 476 # 9. Predict the crosses.
 477 ################################################################################
 478
 479 # If the genotype matrix provides information about individuals for which
 480 # cross prediction is not desired, then the genotype matrix must be subset
 481 # for use in calcCrossMean(). calcCrossMean will return predicted cross
 482 # values for all individuals in the genotype file otherwise.
 483
 484 write("Predict crosses...", stderr())
 485
 486 GP <- G[rownames(G) %in% userPheno[ , userID], ]
 487
 488 print("GP:")
 489 print(head(GP))
 490
 491 write("calcCrossMean...", stderr())
 492
 493 crossPlan <- calcCrossMean(GP,
 494                            ai,
 495                            di,
 496                            userPloidy)
 497
 498
 499 write("Done with calcCrossMean!!!!!!", stderr())
 500
 501
 502
 503 ################################################################################
 504 # 10. Format the information needed for output.
 505 ################################################################################
 506
 507 # Add option to remove crosses with incompatible sexes.
 508
 509
 510 #hash <- new.env(hash = TRUE, parent = emptyenv(), size = 100L)
 511
 512 #assign_hash(userPheno$germplasmName, userPheno$userSexes, hash)
 513
 514 if(!is.na(userSexes)){  # "plant sex estimation 0-4"
 515
 516
 517   write(paste("userSexes", head(userSexes)), stderr())
 518
 519   # Reformat the cross plan
 520   crossPlan <- as.data.frame(crossPlan)
 521
 522   write(paste("CROSSPLAN = ", head(crossPlan)), stderr())
 523   crossPlan <- crossPlan[order(crossPlan[,3], decreasing = TRUE), ] # orders the plan by predicted merit
 524   crossPlan[ ,1] <- rownames(GP)[crossPlan[ ,1]] # replaces internal ID with genotye file ID
 525   crossPlan[ ,2] <- rownames(GP)[crossPlan[ ,2]] # replaces internal ID with genotye file ID
 526   colnames(crossPlan) <- c("Parent1", "Parent2", "CrossPredictedMerit")
 527
 528   write(paste("CROSSPLAN REPLACED = ", head(crossPlan)), stderr());
 529
 530   # Look up the parent sexes and subset
 531   crossPlan$P1Sex <- userPheno[match(crossPlan$Parent1, userPheno$germplasmName), userSexes] # get sexes ordered by Parent1
 532
 533   write(paste("PARENTS1 ", head(crossPlan)), stderr())
 534
 535   crossPlan$P2Sex <- userPheno[match(crossPlan$Parent2, userPheno$germplasmName), userSexes] # get sexes ordered by Parent2
 536
 537   write(paste("PARENTS2 ", head(crossPlan)), stderr())
 538
 539   crossPlan <- crossPlan[!(crossPlan$P1Sex==0 | crossPlan$P2Sex==0),] #remove the 0s
 540   crossPlan <- crossPlan[!(crossPlan$P1Sex==1 & crossPlan$P2Sex==1),] #remove same sex crosses with score of 1
 541   crossPlan <- crossPlan[!(crossPlan$P1Sex==2 & crossPlan$P2Sex==2),] #remove same sex crosses with score of 2
 542
 543   write(paste("CROSSPLAN FILTERED = ", head(crossPlan)), stderr())
 544   #crossPlan <- crossPlan[crossPlan$P1Sex != crossPlan$P2Sex, ] # remove crosses with same-sex parents
 545
 546   # subset the number of crosses the user wishes to output
 547   crossPlan[1:userNCrosses, ]
 548   finalcrosses=crossPlan[1:userNCrosses, ]
 549   outputFile= paste(phenotypeFile, ".out", sep="")
 550
 551   write.csv(finalcrosses, outputFile)
 552
 553 }
 554
 555
 556 if(is.na(userSexes)){
 557
 558   # only subset the number of crosses the user wishes to output
 559   crossPlan <- as.data.frame(crossPlan)
 560
 561   crossPlan <- crossPlan[order(crossPlan[,3], decreasing = TRUE), ] # orders the plan by predicted merit
 562   crossPlan[ ,1] <- rownames(GP)[crossPlan[ ,1]] # replaces internal ID with genotye file ID
 563   crossPlan[ ,2] <- rownames(GP)[crossPlan[ ,2]] # replaces internal ID with genotye file ID
 564   colnames(crossPlan) <- c("Parent1", "Parent2", "CrossPredictedMerit")
 565
 566   crossPlan[1:userNCrosses, ]
 567   finalcrosses=crossPlan[1:userNCrosses, ]
 568   outputFile= paste(phenotypeFile, ".out", sep="")
 569
 570   write.csv(finalcrosses, outputFile)
 571
 572
 573 }