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