R/solGS/gs.r

   1 #SNOPSIS
   2 #calculates genomic estimated breeding values (GEBVs) using rrBLUP,
   3 #GBLUP method
   4
   5 #AUTHOR
   6 # Isaak Y Tecle (iyt2@cornell.edu)
   7
   8 options(echo = FALSE)
   9
  10 library(methods)
  11 library(rrBLUP)
  12 library(plyr)
  13 library(stringr)
  14 #library(lme4)
  15 library(randomForest)
  16 library(parallel)
  17 library(genoDataFilter)
  18 library(phenoAnalysis)
  19 library(caret)
  20 library(dplyr)
  21 library(tibble)
  22 library(rlang)
  23 library(jsonlite)
  24 library(data.table)
  25
  26
  27
  28 allArgs <- commandArgs()
  29
  30
  31 inputFiles  <- scan(grep("input_files", allArgs, value = TRUE),
  32                    what = "character")
  33
  34 outputFiles <- scan(grep("output_files", allArgs, value = TRUE),
  35                     what = "character")
  36
  37 traitsFile <- grep("traits", inputFiles,  value = TRUE)
  38 modelInfoFile  <- grep("model_info", inputFiles, value = TRUE)
  39 message('model_info_file ', modelInfoFile)
  40
  41 modelInfo  <- read.table(modelInfoFile,
  42                          header=TRUE, sep ="\t",
  43                          as.is = c('Value'))
  44
  45 modelInfo  <- column_to_rownames(modelInfo, var="Name")
  46 traitId    <- modelInfo["trait_id", 1]
  47 traitAbbr  <- modelInfo["trait_abbr", 1]
  48 modelId    <- modelInfo["model_id", 1]
  49 protocolId <- modelInfo["protocol_id", 1]
  50
  51 message('class ', class(traitAbbr))
  52 message('trait_id ', traitId)
  53 message('trait_abbr ', traitAbbr)
  54 message('protocol_id ', protocolId)
  55 message('model_id ', modelId)
  56
  57 datasetInfoFile <- grep("dataset_info", inputFiles, value = TRUE)
  58 datasetInfo     <- c()
  59
  60 if (length(datasetInfoFile) != 0 ) {
  61     datasetInfo <- scan(datasetInfoFile, what = "character")
  62     datasetInfo <- paste(datasetInfo, collapse = " ")
  63   } else {
  64     datasetInfo <- c('single population')
  65   }
  66
  67 #validationTrait <- paste("validation", trait, sep = "_")
  68 validationFile  <- grep('validation', outputFiles, value = TRUE)
  69
  70 if (is.null(validationFile)) {
  71   stop("Validation output file is missing.")
  72 }
  73
  74 #kinshipTrait <- paste("rrblup_training_gebvs", trait, sep = "_")
  75 blupFile     <- grep('rrblup_training_gebvs', outputFiles, value = TRUE)
  76
  77 if (is.null(blupFile)) {
  78   stop("GEBVs file is missing.")
  79 }
  80
  81 #markerTrait <- paste("marker_effects", trait, sep = "_")
  82 markerFile  <- grep('marker_effects', outputFiles, value = TRUE)
  83
  84 #traitPhenoFile <- paste("trait_phenotype_data", traitId, sep = "_")
  85 modelPhenoFile <- grep('model_phenodata', outputFiles, value = TRUE)
  86 message('model input trait pheno file ', modelPhenoFile)
  87 varianceComponentsFile <- grep("variance_components", outputFiles, value = TRUE)
  88 filteredGenoFile       <- grep("filtered_genotype_data", outputFiles, value = TRUE)
  89 formattedPhenoFile     <- grep("formatted_phenotype_data", inputFiles, value = TRUE)
  90
  91 genoFile <- grep("genotype_data_", inputFiles, value = TRUE)
  92
  93 if (is.null(genoFile)) {
  94   stop("genotype data file is missing.")
  95 }
  96
  97 if (file.info(genoFile)$size == 0) {
  98   stop("genotype data file is empty.")
  99 }
 100
 101 readFilteredGenoData <- c()
 102 filteredGenoData <- c()
 103 formattedPhenoData <- c()
 104 phenoData          <- c()
 105 genoData           <- c()
 106
 107 if (length(filteredGenoFile) != 0 && file.info(filteredGenoFile)$size != 0) {
 108     filteredGenoData     <- fread(filteredGenoFile,
 109                                   na.strings = c("NA", "", "--", "-"),
 110                                   header = TRUE)
 111
 112     genoData <-  data.frame(filteredGenoData)
 113     genoData <- column_to_rownames(genoData, 'V1')
 114     readFilteredGenoData <- 1
 115 }
 116
 117
 118 if (is.null(filteredGenoData)) {
 119     genoData <- fread(genoFile,
 120                       na.strings = c("NA", "", "--", "-"),
 121                       header = TRUE)
 122
 123     genoData <- unique(genoData, by='V1')
 124     genoData <- data.frame(genoData)
 125     genoData <- column_to_rownames(genoData, 'V1')
 126
 127   #genoDataFilter::filterGenoData
 128     genoData <- convertToNumeric(genoData)
 129     genoData <- filterGenoData(genoData, maf=0.01)
 130     genoData <- roundAlleleDosage(genoData)
 131
 132     filteredGenoData   <- genoData
 133
 134 }
 135
 136 genoData <- genoData[order(row.names(genoData)), ]
 137
 138 if (length(formattedPhenoFile) != 0 && file.info(formattedPhenoFile)$size != 0) {
 139     formattedPhenoData <- data.frame(fread(formattedPhenoFile,
 140                                            header = TRUE,
 141                                            na.strings = c("NA", "", "--", "-", ".")
 142                                             ))
 143
 144 } else {
 145
 146     if (datasetInfo == 'combined populations') {
 147
 148          phenoFile <- grep("model_phenodata", inputFiles, value = TRUE)
 149     } else {
 150
 151         phenoFile <- grep("\\/phenotype_data", inputFiles, value = TRUE)
 152     }
 153
 154     if (is.null(phenoFile)) {
 155         stop("phenotype data file is missing.")
 156     }
 157
 158     if (file.info(phenoFile)$size == 0) {
 159         stop("phenotype data file is empty.")
 160     }
 161
 162     phenoData <- data.frame(fread(phenoFile,
 163                                   sep = "\t",
 164                                   na.strings = c("NA", "", "--", "-", "."),
 165                                   header = TRUE))
 166 }
 167
 168 phenoTrait <- c()
 169
 170 if (datasetInfo == 'combined populations') {
 171
 172    if (!is.null(formattedPhenoData)) {
 173       phenoTrait <- subset(formattedPhenoData, select = traitAbbr)
 174       phenoTrait <- na.omit(phenoTrait)
 175
 176   } else {
 177
 178       if (any(grepl('Average', names(phenoData)))) {
 179           phenoTrait <- phenoData %>% select(V1, Average) %>% data.frame
 180       } else {
 181           phenoTrait <- phenoData
 182       }
 183
 184       colnames(phenoTrait)  <- c('genotypes', traitAbbr)
 185   }
 186  } else {
 187
 188      if (!is.null(formattedPhenoData)) {
 189          phenoTrait <- subset(formattedPhenoData, select = c('V1', traitAbbr))
 190          phenoTrait <- as.data.frame(phenoTrait)
 191          phenoTrait <- na.omit(phenoTrait)
 192          print(head(phenoTrait))
 193          colnames(phenoTrait)[1] <- 'genotypes'
 194
 195      } else if (length(grep('list', phenoFile)) != 0) {
 196  message('phenoTrait traitAbbr ', traitAbbr)
 197          phenoTrait <- averageTrait(phenoData, traitAbbr)
 198
 199      } else {
 200          print(head(phenoTrait))
 201           print(head(phenoData))
 202          message('phenoTrait trait_abbr ', traitAbbr)
 203          print(class(traitAbbr))
 204          print(traitAbbr)
 205          phenoTrait <- getAdjMeans(phenoData,
 206                                    traitName = traitAbbr,
 207                                    calcAverages = TRUE)
 208      }
 209 }
 210
 211 print('phenoTrait')
 212 print(head(phenoTrait))
 213 meanType <- names(phenoTrait)[2]
 214 names(phenoTrait)  <- c('genotypes', traitAbbr)
 215
 216 selectionTempFile <- grep("selection_population", inputFiles, value = TRUE)
 217
 218 selectionFile       <- c()
 219 filteredPredGenoFile <- c()
 220 selectionAllFiles   <- c()
 221
 222 if (length(selectionTempFile) !=0 ) {
 223   selectionAllFiles <- scan(selectionTempFile, what = "character")
 224
 225   selectionFile <- grep("\\/genotype_data", selectionAllFiles, value = TRUE)
 226
 227   #filteredPredGenoFile   <- grep("filtered_genotype_data_",  selectionAllFiles, value = TRUE)
 228 }
 229
 230 selectionPopGEBVsFile <- grep("rrblup_selection_gebvs", outputFiles, value = TRUE)
 231
 232 selectionData            <- c()
 233 readFilteredPredGenoData <- c()
 234 filteredPredGenoData     <- c()
 235
 236 ## if (length(filteredPredGenoFile) != 0 && file.info(filteredPredGenoFile)$size != 0) {
 237 ##   selectionData <- fread(filteredPredGenoFile, na.strings = c("NA", " ", "--", "-"),)
 238 ##   readFilteredPredGenoData <- 1
 239
 240 ##   selectionData           <- data.frame(selectionData)
 241 ##   rownames(selectionData) <- selectionData[, 1]
 242 ##   selectionData[, 1]      <- NULL
 243
 244 ## } else
 245 if (length(selectionFile) != 0) {
 246
 247     selectionData <- fread(selectionFile,
 248                            header = TRUE,
 249                            na.strings = c("NA", "", "--", "-"))
 250
 251     selectionData <- unique(selectionData, by='V1')
 252     selectionData <- data.frame(selectionData)
 253     selectionData <- column_to_rownames(selectionData, 'V1')
 254
 255     selectionData <- convertToNumeric(selectionData)
 256     selectionData <- filterGenoData(selectionData, maf=0.01)
 257     selectionData <- roundAlleleDosage(selectionData)
 258
 259     filteredPredGenoData <- selectionData
 260 }
 261
 262
 263 #impute genotype values for obs with missing values,
 264 genoDataMissing <- c()
 265
 266 if (sum(is.na(genoData)) > 0) {
 267   genoDataMissing<- c('yes')
 268
 269   genoData <- na.roughfix(genoData)
 270   genoData <- data.frame(genoData)
 271 }
 272
 273 #create phenotype and genotype datasets with
 274 #common stocks only
 275
 276 #extract observation lines with both
 277 #phenotype and genotype data only.
 278 commonObs           <- intersect(phenoTrait$genotypes, row.names(genoData))
 279 commonObs           <- data.frame(commonObs)
 280 rownames(commonObs) <- commonObs[, 1]
 281
 282 #include in the genotype dataset only phenotyped lines
 283 genoDataFilteredObs <- genoData[(rownames(genoData) %in% rownames(commonObs)), ]
 284
 285 #drop phenotyped lines without genotype data
 286 phenoTrait <- phenoTrait[(phenoTrait$genotypes %in% rownames(commonObs)), ]
 287
 288 phenoTraitMarker           <- data.frame(phenoTrait)
 289 rownames(phenoTraitMarker) <- phenoTraitMarker[, 1]
 290 phenoTraitMarker[, 1]      <- NULL
 291
 292 #impute missing data in prediction data
 293 selectionDataMissing <- c()
 294 if (length(selectionData) != 0) {
 295   #purge markers unique to both populations
 296   commonMarkers       <- intersect(names(data.frame(genoDataFilteredObs)), names(selectionData))
 297   selectionData      <- subset(selectionData, select = commonMarkers)
 298   genoDataFilteredObs <- subset(genoDataFilteredObs, select= commonMarkers)
 299
 300   if (sum(is.na(selectionData)) > 0) {
 301     selectionDataMissing <- c('yes')
 302     selectionData <- na.roughfix(selectionData)
 303     selectionData <- data.frame(selectionData)
 304   }
 305 }
 306
 307 #change genotype coding to [-1, 0, 1], to use the A.mat ) if  [0, 1, 2]
 308 genoTrCode <- grep("2", genoDataFilteredObs[1, ], value = TRUE)
 309 if(length(genoTrCode) != 0) {
 310   genoData            <- genoData - 1
 311   genoDataFilteredObs <- genoDataFilteredObs - 1
 312 }
 313
 314 if (length(selectionData) != 0 ) {
 315   genoSlCode <- grep("2", selectionData[1, ], value = TRUE)
 316   if (length(genoSlCode) != 0 ) {
 317     selectionData <- selectionData - 1
 318   }
 319 }
 320
 321 ordered.markerEffects <- c()
 322 trGEBV                <- c()
 323 validationAll         <- c()
 324 combinedGebvsFile     <- c()
 325 allGebvs              <- c()
 326 modelPhenoData        <- c()
 327 relationshipMatrix    <- c()
 328
 329 #additive relationship model
 330 #calculate the inner products for
 331 #genotypes (realized relationship matrix)
 332 relationshipMatrixFile <- grep("relationship_matrix_table", outputFiles, value = TRUE)
 333 relationshipMatrixJsonFile <- grep("relationship_matrix_json", outputFiles, value = TRUE)
 334
 335 traitRelationshipMatrixFile <- grep("relationship_matrix_adjusted_table", outputFiles, value = TRUE)
 336 traitRelationshipMatrixJsonFile <- grep("relationship_matrix_adjusted_json", outputFiles, value = TRUE)
 337
 338 inbreedingFile <- grep('inbreeding_coefficients', outputFiles, value=TRUE)
 339 aveKinshipFile <- grep('average_kinship', outputFiles, value=TRUE)
 340
 341 inbreeding <- c()
 342 aveKinship <- c()
 343
 344 if (length(relationshipMatrixFile) != 0) {
 345   if (file.info(relationshipMatrixFile)$size > 0 ) {
 346       relationshipMatrix <- data.frame(fread(relationshipMatrixFile,
 347                          header = TRUE))
 348
 349       rownames(relationshipMatrix) <- relationshipMatrix[, 1]
 350       relationshipMatrix[, 1]      <- NULL
 351       colnames(relationshipMatrix) <- rownames(relationshipMatrix)
 352       relationshipMatrix           <- data.matrix(relationshipMatrix)
 353
 354   } else {
 355     relationshipMatrix           <- A.mat(genoData)
 356     diag(relationshipMatrix)     <- diag(relationshipMatrix) + 1e-6
 357
 358     inbreeding <- diag(relationshipMatrix)
 359     inbreeding <- inbreeding - 1
 360
 361     inbreeding <- inbreeding %>% replace(., . < 0, 0)
 362     inbreeding <- data.frame(inbreeding)
 363
 364     inbreeding <- inbreeding %>%
 365         rownames_to_column('genotypes') %>%
 366         rename(Inbreeding = inbreeding) %>%
 367         arrange(Inbreeding) %>%
 368         mutate_at('Inbreeding', round, 3) %>%
 369         column_to_rownames('genotypes')
 370   }
 371 }
 372
 373 relationshipMatrix <- data.frame(relationshipMatrix)
 374 colnames(relationshipMatrix) <- rownames(relationshipMatrix)
 375
 376 relationshipMatrix <- rownames_to_column(relationshipMatrix, var="genotypes")
 377 relationshipMatrix <- relationshipMatrix %>% mutate_if(is.numeric, round, 3)
 378 relationshipMatrix <- column_to_rownames(relationshipMatrix, var="genotypes")
 379
 380 traitRelationshipMatrix <- relationshipMatrix[(rownames(relationshipMatrix) %in% rownames(commonObs)), ]
 381 traitRelationshipMatrix <- traitRelationshipMatrix[, (colnames(traitRelationshipMatrix) %in% rownames(commonObs))]
 382
 383 traitRelationshipMatrix <- data.matrix(traitRelationshipMatrix)
 384
 385 #relationshipMatrixFiltered <- relationshipMatrixFiltered + 1e-3
 386
 387 nCores <- detectCores()
 388
 389 if (nCores > 1) {
 390   nCores <- (nCores %/% 2)
 391 } else {
 392   nCores <- 1
 393 }
 394
 395
 396 if (length(selectionData) == 0) {
 397
 398   trModel  <- kin.blup(data   = phenoTrait,
 399                       geno   = 'genotypes',
 400                       pheno  = traitAbbr,
 401                       K      = traitRelationshipMatrix,
 402                       n.core = nCores,
 403                       PEV    = TRUE
 404                      )
 405
 406   trGEBV    <- trModel$g
 407   trGEBVPEV <- trModel$PEV
 408   trGEBVSE  <- sqrt(trGEBVPEV)
 409   trGEBVSE  <- data.frame(round(trGEBVSE, 2))
 410
 411   trGEBV <- data.frame(round(trGEBV, 2))
 412
 413   colnames(trGEBVSE) <- c('SE')
 414   colnames(trGEBV) <- traitAbbr
 415
 416   trGEBVSE <- rownames_to_column(trGEBVSE, var="genotypes")
 417   trGEBV   <- rownames_to_column(trGEBV, var="genotypes")
 418
 419   trGEBVSE <- full_join(trGEBV, trGEBVSE)
 420
 421   trGEBVSE <-  trGEBVSE %>% arrange_(.dots= paste0('desc(', traitAbbr, ')'))
 422
 423   trGEBVSE <- column_to_rownames(trGEBVSE, var="genotypes")
 424
 425   trGEBV <- trGEBV %>% arrange_(.dots = paste0('desc(', traitAbbr, ')'))
 426   trGEBV <- column_to_rownames(trGEBV, var="genotypes")
 427
 428   phenoTraitMarker    <- data.matrix(phenoTraitMarker)
 429   genoDataFilteredObs <- data.matrix(genoDataFilteredObs)
 430
 431   markerEffects <- mixed.solve(y = phenoTraitMarker,
 432                                Z = genoDataFilteredObs
 433                                )
 434
 435   ordered.markerEffects <- data.matrix(markerEffects$u)
 436   ordered.markerEffects <- data.matrix(ordered.markerEffects [order (-ordered.markerEffects[, 1]), ])
 437   ordered.markerEffects <- round(ordered.markerEffects, 5)
 438
 439   colnames(ordered.markerEffects) <- c("Marker Effects")
 440   ordered.markerEffects <- data.frame(ordered.markerEffects)
 441
 442
 443   modelPhenoData   <- data.frame(round(phenoTraitMarker, 2))
 444
 445   heritability  <- round((trModel$Vg/(trModel$Ve + trModel$Vg)), 2)
 446
 447   cat("\n", file = varianceComponentsFile,  append = FALSE)
 448   cat('Additive genetic variance',  trModel$Vg, file = varianceComponentsFile, sep = '\t', append = TRUE)
 449   cat("\n", file = varianceComponentsFile,  append = TRUE)
 450   cat('Error variance', trModel$Ve, file = varianceComponentsFile, sep = "\t", append = TRUE)
 451   cat("\n", file = varianceComponentsFile,  append = TRUE)
 452   cat('SNP heritability (h)', heritability, file = varianceComponentsFile, sep = '\t', append = TRUE)
 453
 454   combinedGebvsFile <- grep('selected_traits_gebv', outputFiles, ignore.case = TRUE,value = TRUE)
 455
 456   if (length(combinedGebvsFile) != 0) {
 457       fileSize <- file.info(combinedGebvsFile)$size
 458       if (fileSize != 0 ) {
 459           combinedGebvs <- data.frame(fread(combinedGebvsFile,
 460                                             header = TRUE))
 461
 462         rownames(combinedGebvs) <- combinedGebvs[,1]
 463           combinedGebvs[,1]       <- NULL
 464
 465           allGebvs <- merge(combinedGebvs, trGEBV,
 466                             by = 0,
 467                             all = TRUE
 468                             )
 469
 470           rownames(allGebvs) <- allGebvs[,1]
 471           allGebvs[,1] <- NULL
 472       }
 473   }
 474
 475 #cross-validation
 476
 477   if (is.null(selectionFile)) {
 478       genoNum <- nrow(phenoTrait)
 479
 480       if (genoNum < 20 ) {
 481           warning(genoNum, " is too small number of genotypes.")
 482       }
 483
 484       set.seed(4567)
 485
 486       k <- 10
 487       times <- 2
 488       cvFolds <- createMultiFolds(phenoTrait[, 2], k=k, times=times)
 489
 490       for ( r in 1:times) {
 491           re <- paste0('Rep', r)
 492
 493           for (i in 1:k) {
 494               fo <- ifelse(i < 10, 'Fold0', 'Fold')
 495
 496               trFoRe <- paste0(fo, i, '.', re)
 497               trG <- cvFolds[[trFoRe]]
 498               slG <- as.numeric(rownames(phenoTrait[-trG,]))
 499
 500               kblup <- paste("rKblup", i, sep = ".")
 501
 502               result <- kin.blup(data  = phenoTrait[trG,],
 503                                  geno  = 'genotypes',
 504                                  pheno = traitAbbr,
 505                                  K     = traitRelationshipMatrix,
 506                                  n.core = nCores,
 507                                  PEV    = TRUE
 508                                  )
 509
 510               assign(kblup, result)
 511
 512                                         #calculate cross-validation accuracy
 513               valBlups   <- result$g
 514
 515               valBlups   <- data.frame(valBlups)
 516
 517               slG <- slG[which(slG <= nrow(phenoTrait))]
 518
 519               slGDf <- phenoTrait[(rownames(phenoTrait) %in% slG),]
 520               rownames(slGDf) <- slGDf[, 1]
 521               slGDf[, 1] <- NULL
 522
 523               valBlups <-  rownames_to_column(valBlups, var="genotypes")
 524               slGDf    <-  rownames_to_column(slGDf, var="genotypes")
 525
 526               valCorData <- inner_join(slGDf, valBlups, by="genotypes")
 527               valCorData$genotypes <- NULL
 528
 529               accuracy   <- try(cor(valCorData))
 530               validation <- paste("validation", trFoRe, sep = ".")
 531               cvTest <- paste("CV", trFoRe, sep = " ")
 532
 533               if ( class(accuracy) != "try-error")
 534               {
 535                   accuracy <- round(accuracy[1,2], digits = 3)
 536                   accuracy <- data.matrix(accuracy)
 537
 538                   colnames(accuracy) <- c("correlation")
 539                   rownames(accuracy) <- cvTest
 540
 541                   assign(validation, accuracy)
 542
 543                   if (!is.na(accuracy[1,1])) {
 544                       validationAll <- rbind(validationAll, accuracy)
 545                   }
 546               }
 547           }
 548       }
 549
 550       validationAll <- data.frame(validationAll[order(-validationAll[, 1]), ])
 551       colnames(validationAll) <- c('Correlation')
 552   }
 553 }
 554
 555 selectionPopResult <- c()
 556 selectionPopGEBVs  <- c()
 557 selectionPopGEBVSE <- c()
 558
 559 if (length(selectionData) != 0) {
 560
 561     genoDataTrSl <- rbind(genoDataFilteredObs, selectionData)
 562     rTrSl <- A.mat(genoDataTrSl)
 563
 564     selectionPopResult <- kin.blup(data   = phenoTrait,
 565                                     geno   = 'genotypes',
 566                                     pheno  = traitAbbr,
 567                                     K      = rTrSl,
 568                                     n.core = nCores,
 569                                     PEV    = TRUE
 570                                     )
 571
 572     selectionPopGEBVs <- round(data.frame(selectionPopResult$g), 2)
 573     colnames(selectionPopGEBVs) <- traitAbbr
 574     selectionPopGEBVs <- rownames_to_column(selectionPopGEBVs, var="genotypes")
 575
 576     selectionPopPEV <- selectionPopResult$PEV
 577     selectionPopSE  <- sqrt(selectionPopPEV)
 578     selectionPopSE  <- data.frame(round(selectionPopSE, 2))
 579     colnames(selectionPopSE) <- 'SE'
 580     genotypesSl     <- rownames(selectionData)
 581
 582     selectionPopSE <- rownames_to_column(selectionPopSE, var="genotypes")
 583     selectionPopSE <-  selectionPopSE %>% filter(genotypes %in% genotypesSl)
 584
 585     selectionPopGEBVs <-  selectionPopGEBVs %>% filter(genotypes %in% genotypesSl)
 586
 587     selectionPopGEBVSE <- inner_join(selectionPopGEBVs, selectionPopSE, by="genotypes")
 588
 589     sortVar <- parse_quosure(traitAbbr)
 590     selectionPopGEBVs <- selectionPopGEBVs %>% arrange(desc((!!sortVar)))
 591     selectionPopGEBVs <- column_to_rownames(selectionPopGEBVs, var="genotypes")
 592
 593     selectionPopGEBVSE <-  selectionPopGEBVSE %>% arrange(desc((!!sortVar)))
 594     selectionPopGEBVSE <- column_to_rownames(selectionPopGEBVSE, var="genotypes")
 595 }
 596
 597 if (!is.null(selectionPopGEBVs) & length(selectionPopGEBVsFile) != 0)  {
 598     fwrite(selectionPopGEBVs,
 599            file  = selectionPopGEBVsFile,
 600            row.names = TRUE,
 601            sep   = "\t",
 602            quote = FALSE,
 603            )
 604 }
 605
 606 if(!is.null(validationAll)) {
 607     fwrite(validationAll,
 608            file  = validationFile,
 609            row.names = TRUE,
 610            sep   = "\t",
 611            quote = FALSE,
 612            )
 613 }
 614
 615
 616 if (!is.null(ordered.markerEffects)) {
 617     fwrite(ordered.markerEffects,
 618            file  = markerFile,
 619            row.names = TRUE,
 620            sep   = "\t",
 621            quote = FALSE,
 622            )
 623 }
 624
 625
 626 if (!is.null(trGEBV)) {
 627     fwrite(trGEBV,
 628            file  = blupFile,
 629            row.names = TRUE,
 630            sep   = "\t",
 631            quote = FALSE,
 632            )
 633 }
 634
 635 if (length(combinedGebvsFile) != 0 ) {
 636     if(file.info(combinedGebvsFile)$size == 0) {
 637         fwrite(trGEBV,
 638                file  = combinedGebvsFile,
 639                row.names = TRUE,
 640                sep   = "\t",
 641                quote = FALSE,
 642                )
 643       } else {
 644       fwrite(allGebvs,
 645              file  = combinedGebvsFile,
 646              row.names = TRUE,
 647              sep   = "\t",
 648              quote = FALSE,
 649              )
 650     }
 651 }
 652
 653
 654 if (!is.null(modelPhenoData) & length(modelPhenoFile) != 0) {
 655
 656     if (!is.null(meanType)) {
 657         colnames(modelPhenoData) <- meanType
 658     }
 659
 660     fwrite(modelPhenoData,
 661            file  = modelPhenoFile,
 662            row.names = TRUE,
 663            sep   = "\t",
 664            quote = FALSE,
 665            )
 666 }
 667
 668
 669 if (!is.null(filteredGenoData) && is.null(readFilteredGenoData)) {
 670   fwrite(filteredGenoData,
 671          file  = filteredGenoFile,
 672          row.names = TRUE,
 673          sep   = "\t",
 674          quote = FALSE,
 675          )
 676
 677 }
 678
 679 ## if (length(filteredPredGenoFile) != 0 && is.null(readFilteredPredGenoData)) {
 680 ##   fwrite(filteredPredGenoData,
 681 ##          file  = filteredPredGenoFile,
 682 ##          row.names = TRUE,
 683 ##          sep   = "\t",
 684 ##          quote = FALSE,
 685 ##          )
 686 ## }
 687
 688 ## if (!is.null(genoDataMissing)) {
 689 ##   write.table(genoData,
 690 ##               file = genoFile,
 691 ##               sep = "\t",
 692 ##               col.names = NA,
 693 ##               quote = FALSE,
 694 ##             )
 695
 696 ## }
 697
 698 ## if (!is.null(predictionDataMissing)) {
 699 ##   write.table(predictionData,
 700 ##               file = predictionFile,
 701 ##               sep = "\t",
 702 ##               col.names = NA,
 703 ##               quote = FALSE,
 704 ##               )
 705 ## }
 706
 707
 708 if (file.info(relationshipMatrixFile)$size == 0) {
 709
 710   fwrite(relationshipMatrix,
 711          file  = relationshipMatrixFile,
 712          row.names = TRUE,
 713          sep   = "\t",
 714          quote = FALSE,
 715          )
 716 }
 717
 718 if (file.info(relationshipMatrixJsonFile)$size == 0) {
 719
 720     relationshipMatrixJson <- relationshipMatrix
 721     relationshipMatrixJson[upper.tri(relationshipMatrixJson)] <- NA
 722
 723
 724     relationshipMatrixJson <- data.frame(relationshipMatrixJson)
 725
 726     relationshipMatrixList <- list(labels = names(relationshipMatrixJson),
 727                                        values = relationshipMatrixJson)
 728
 729     relationshipMatrixJson <- jsonlite::toJSON(relationshipMatrixList)
 730
 731
 732     write(relationshipMatrixJson,
 733                     file  = relationshipMatrixJsonFile,
 734                     )
 735 }
 736
 737
 738 if (file.info(traitRelationshipMatrixFile)$size == 0) {
 739
 740     inbre <- diag(traitRelationshipMatrix)
 741     inbre <- inbre - 1
 742
 743     diag(traitRelationshipMatrix) <- inbre
 744
 745     traitRelationshipMatrix <- data.frame(traitRelationshipMatrix) %>% replace(., . < 0, 0)
 746
 747     fwrite(traitRelationshipMatrix,
 748            file  = traitRelationshipMatrixFile,
 749            row.names = TRUE,
 750            sep   = "\t",
 751            quote = FALSE,
 752            )
 753
 754     if (file.info(traitRelationshipMatrixJsonFile)$size == 0) {
 755
 756         traitRelationshipMatrixJson <- traitRelationshipMatrix
 757         traitRelationshipMatrixJson[upper.tri(traitRelationshipMatrixJson)] <- NA
 758
 759         traitRelationshipMatrixJson <- data.frame(traitRelationshipMatrixJson)
 760
 761         traitRelationshipMatrixList <- list(labels = names(traitRelationshipMatrixJson),
 762                                             values = traitRelationshipMatrixJson)
 763
 764         traitRelationshipMatrixJson <- jsonlite::toJSON(traitRelationshipMatrixList)
 765
 766         write(traitRelationshipMatrixJson,
 767               file  = traitRelationshipMatrixJsonFile,
 768               )
 769     }
 770 }
 771
 772
 773 if (file.info(inbreedingFile)$size == 0) {
 774
 775   fwrite(inbreeding,
 776          file  = inbreedingFile,
 777          row.names = TRUE,
 778          sep   = "\t",
 779          quote = FALSE,
 780          )
 781 }
 782
 783
 784 if (file.info(aveKinshipFile)$size == 0) {
 785
 786     aveKinship <- data.frame(apply(traitRelationshipMatrix, 1, mean))
 787
 788     aveKinship<- aveKinship %>%
 789         rownames_to_column('genotypes') %>%
 790         rename(Mean_kinship = contains('traitRe')) %>%
 791         arrange(Mean_kinship) %>%
 792         mutate_at('Mean_kinship', round, 3) %>%
 793         column_to_rownames('genotypes')
 794
 795     fwrite(aveKinship,
 796            file  = aveKinshipFile,
 797            row.names = TRUE,
 798            sep   = "\t",
 799            quote = FALSE,
 800            )
 801 }
 802
 803
 804 if (file.info(formattedPhenoFile)$size == 0 && !is.null(formattedPhenoData) ) {
 805   fwrite(formattedPhenoData,
 806          file = formattedPhenoFile,
 807          row.names = TRUE,
 808          sep = "\t",
 809          quote = FALSE,
 810          )
 811 }
 812
 813 message("Done.")
 814
 815 q(save = "no", runLast = FALSE)