R/gs.r

   1 #SNOPSIS
   2 #calculates genomic estimated breeding values (GEBVs) using rrBLUP
   3
   4 #AUTHOR
   5 # Isaak Y Tecle (iyt2@cornell.edu)
   6
   7 options(echo = FALSE)
   8
   9 library(rrBLUP)
  10 library(plyr)
  11 library(stringr)
  12 library(lme4)
  13 library(randomForest)
  14 library(data.table)
  15
  16
  17 allArgs <- commandArgs()
  18
  19 inFile <- grep("input_files",
  20                allArgs,
  21                ignore.case = TRUE,
  22                perl = TRUE,
  23                value = TRUE
  24                )
  25
  26 outFile <- grep("output_files",
  27                 allArgs,
  28                 ignore.case = TRUE,
  29                 perl = TRUE,
  30                 value = TRUE
  31                 )
  32
  33 outFiles <- scan(outFile,
  34                  what = "character"
  35                  )
  36
  37 inFiles <- scan(inFile,
  38                 what = "character"
  39                 )
  40
  41
  42 traitsFile <- grep("traits",
  43                    inFiles,
  44                    ignore.case = TRUE,
  45                    fixed = FALSE,
  46                    value = TRUE
  47                    )
  48
  49 traitFile <- grep("trait_info",
  50                    inFiles,
  51                    ignore.case = TRUE,
  52                    fixed = FALSE,
  53                    value = TRUE
  54                   )
  55
  56
  57 traitInfo <- scan(traitFile,
  58                   what = "character",
  59                   )
  60
  61 traitInfo <- strsplit(traitInfo, "\t");
  62 traitId   <- traitInfo[[1]]
  63 trait     <- traitInfo[[2]]
  64
  65 datasetInfoFile <- grep("dataset_info",
  66                         inFiles,
  67                         ignore.case = TRUE,
  68                         fixed = FALSE,
  69                         value = TRUE
  70                         )
  71 datasetInfo <- c()
  72
  73 if (length(datasetInfoFile) != 0 ) {
  74     datasetInfo <- scan(datasetInfoFile,
  75                         what= "character"
  76                         )
  77
  78     datasetInfo <- paste(datasetInfo, collapse = " ")
  79
  80   } else {
  81     datasetInfo <- c('single population')
  82
  83   }
  84
  85 validationTrait <- paste("validation", trait, sep = "_")
  86
  87 validationFile  <- grep(validationTrait,
  88                         outFiles,
  89                         ignore.case=TRUE,
  90                         fixed = FALSE,
  91                         value=TRUE
  92                         )
  93
  94 kinshipTrait <- paste("kinship", trait, sep = "_")
  95
  96 blupFile <- grep(kinshipTrait,
  97                  outFiles,
  98                  ignore.case = TRUE,
  99                  fixed = FALSE,
 100                  value = TRUE
 101                  )
 102
 103 markerTrait <- paste("marker", trait, sep = "_")
 104 markerFile  <- grep(markerTrait,
 105                    outFiles,
 106                    ignore.case = TRUE,
 107                    fixed = FALSE,
 108                    value = TRUE
 109                    )
 110
 111 traitPhenoFile <- paste("phenotype_trait", trait, sep = "_")
 112 traitPhenoFile <- grep(traitPhenoFile,
 113                        outFiles,
 114                        ignore.case = TRUE,
 115                        fixed = FALSE,
 116                        value = TRUE
 117                        )
 118
 119 varianceComponentsFile <- grep("variance_components",
 120                                outFiles,
 121                                ignore.case = TRUE,
 122                                fixed = FALSE,
 123                                value = TRUE
 124                                )
 125
 126 formattedPhenoFile <- grep("formatted_phenotype_data",
 127                            inFiles,
 128                            ignore.case = TRUE,
 129                            fixed = FALSE,
 130                            value = TRUE
 131                            )
 132
 133 formattedPhenoData <- c()
 134 phenoData          <- c()
 135
 136 if (length(formattedPhenoFile) != 0 && file.info(formattedPhenoFile)$size != 0) {
 137     formattedPhenoData <- as.data.frame(fread(formattedPhenoFile,
 138                                               na.strings = c("NA", " ", "--", "-", ".")
 139                                               ))
 140
 141
 142     row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 143     formattedPhenoData[, 1]       <- NULL
 144
 145 } else {
 146   phenoFile <- grep("\\/phenotype_data",
 147                     inFiles,
 148                     ignore.case = TRUE,
 149                     fixed = FALSE,
 150                     value = TRUE,
 151                     perl = TRUE,
 152                     )
 153
 154   phenoData <- fread(phenoFile,
 155                      header=TRUE,
 156                      na.strings = c("NA", " ", "--", "-", "."),
 157                      )
 158
 159 }
 160
 161 phenoData  <- as.data.frame(phenoData)
 162 phenoTrait <- c()
 163
 164 if (datasetInfo == 'combined populations') {
 165
 166    if (!is.null(formattedPhenoData)) {
 167       phenoTrait <- subset(formattedPhenoData, select=trait)
 168       phenoTrait <- na.omit(phenoTrait)
 169
 170     } else {
 171       dropColumns <- grep(trait,
 172                           names(phenoData),
 173                           ignore.case = TRUE,
 174                           value = TRUE,
 175                           fixed = FALSE
 176                           )
 177
 178       phenoTrait <- phenoData[,!(names(phenoData) %in% dropColumns)]
 179
 180       phenoTrait            <- as.data.frame(phenoTrait)
 181       row.names(phenoTrait) <- phenoTrait[, 1]
 182       phenoTrait[, 1]       <- NULL
 183       colnames(phenoTrait)  <- trait
 184     }
 185
 186 } else {
 187
 188   if (!is.null(formattedPhenoData)) {
 189     phenoTrait <- subset(formattedPhenoData, select=trait)
 190     phenoTrait <- na.omit(phenoTrait)
 191
 192   } else {
 193     dropColumns <- c("uniquename", "stock_name")
 194     phenoData   <- phenoData[,!(names(phenoData) %in% dropColumns)]
 195
 196     phenoTrait <- subset(phenoData,
 197                          select = c("object_name", "object_id", "design", "block", "replicate", trait)
 198                          )
 199
 200     experimentalDesign <- phenoTrait[2, 'design']
 201
 202     if (class(phenoTrait[, trait]) != 'numeric') {
 203       phenoTrait[, trait] <- as.numeric(as.character(phenoTrait[, trait]))
 204     }
 205
 206     if (is.na(experimentalDesign) == TRUE) {experimentalDesign <- c('No Design')}
 207
 208     if ((experimentalDesign == 'Augmented' || experimentalDesign == 'RCBD')  &&  unique(phenoTrait$block) > 1) {
 209
 210       message("GS experimental design: ", experimentalDesign)
 211
 212       augData <- subset(phenoTrait,
 213                         select = c("object_name", "object_id",  "block",  trait)
 214                         )
 215
 216       colnames(augData)[1] <- "genotypes"
 217       colnames(augData)[4] <- "trait"
 218
 219       model <- try(lmer(trait ~ 0 + genotypes + (1|block),
 220                         augData,
 221                         na.action = na.omit
 222                         ))
 223
 224       if (class(model) != "try-error") {
 225         phenoTrait <- data.frame(fixef(model))
 226
 227         colnames(phenoTrait) <- trait
 228
 229         nn <- gsub('genotypes', '', rownames(phenoTrait))
 230         rownames(phenoTrait) <- nn
 231
 232         phenoTrait <- round(phenoTrait, digits = 2)
 233       }
 234
 235     } else if (experimentalDesign == 'Alpha') {
 236
 237       message("Experimental desgin: ", experimentalDesign)
 238
 239       alphaData <- subset(phenoData,
 240                             select = c("object_name", "object_id","block", "replicate", trait)
 241                             )
 242
 243       colnames(alphaData)[1] <- "genotypes"
 244       colnames(alphaData)[5] <- "trait"
 245
 246       model <- try(lmer(trait ~ 0 + genotypes + (1|replicate/block),
 247                         alphaData,
 248                         na.action = na.omit
 249                         ))
 250
 251       if (class(model) != "try-error") {
 252         phenoTrait <- data.frame(fixef(model))
 253
 254         colnames(phenoTrait) <- trait
 255
 256         nn <- gsub('genotypes', '', rownames(phenoTrait))
 257         rownames(phenoTrait) <- nn
 258
 259         phenoTrait <- round(phenoTrait, digits = 2)
 260
 261       }
 262
 263     } else {
 264
 265       phenoTrait <- subset(phenoData,
 266                            select = c("object_name", "object_id",  trait)
 267                            )
 268
 269       if (sum(is.na(phenoTrait)) > 0) {
 270         message("No. of pheno missing values: ", sum(is.na(phenoTrait)))
 271         phenoTrait <- na.omit(phenoTrait)
 272       }
 273
 274         #calculate mean of reps/plots of the same accession and
 275         #create new df with the accession means
 276
 277       phenoTrait   <- phenoTrait[order(row.names(phenoTrait)), ]
 278       phenoTrait   <- data.frame(phenoTrait)
 279       message('phenotyped lines before averaging: ', length(row.names(phenoTrait)))
 280
 281       phenoTrait<-ddply(phenoTrait, "object_name", colwise(mean))
 282       message('phenotyped lines after averaging: ', length(row.names(phenoTrait)))
 283
 284       phenoTrait <- subset(phenoTrait, select=c("object_name", trait))
 285       row.names(phenoTrait) <- phenoTrait[, 1]
 286       phenoTrait[, 1] <- NULL
 287
 288         #format all-traits population phenotype dataset
 289         ## formattedPhenoData <- phenoData
 290         ## dropColumns <- c("object_id", "stock_id", "design", "block", "replicate" )
 291
 292         ## formattedPhenoData <- formattedPhenoData[, !(names(formattedPhenoData) %in% dropColumns)]
 293         ## formattedPhenoData <- ddply(formattedPhenoData,
 294         ##                             "object_name",
 295         ##                             colwise(mean)
 296         ##                             )
 297
 298         ## row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 299         ## formattedPhenoData[, 1] <- NULL
 300
 301         ## formattedPhenoData <- round(formattedPhenoData,
 302         ##                             digits=3
 303         ##                             )
 304     }
 305   }
 306 }
 307
 308 genoFile <- grep("genotype_data",
 309                  inFiles,
 310                  ignore.case = TRUE,
 311                  fixed = FALSE,
 312                  value = TRUE
 313                  )
 314
 315 genoData <- fread(genoFile,
 316                   header=TRUE,
 317                   na.strings = c("NA", " ", "--", "-"),
 318                   )
 319
 320 genoData           <- as.data.frame(genoData)
 321 rownames(genoData) <- genoData[, 1]
 322 genoData[, 1]      <- NULL
 323 genoData           <- genoData[, colSums(is.na(genoData)) < nrow(genoData) * 0.5]
 324
 325 predictionTempFile <- grep("prediction_population",
 326                        inFiles,
 327                        ignore.case = TRUE,
 328                        fixed = FALSE,
 329                        value = TRUE
 330                        )
 331
 332 predictionFile <- c()
 333
 334 message('prediction temp genotype file: ', predictionTempFile)
 335
 336 if (length(predictionTempFile) !=0 ) {
 337   predictionFile <- scan(predictionTempFile,
 338                          what="character"
 339                          )
 340 }
 341
 342 message('prediction genotype file: ', predictionFile)
 343
 344 predictionPopGEBVsFile <- grep("prediction_pop_gebvs",
 345                                outFiles,
 346                                ignore.case = TRUE,
 347                                fixed = FALSE,
 348                                value = TRUE
 349                                )
 350 message("prediction gebv file: ",  predictionPopGEBVsFile)
 351
 352 predictionData <- c()
 353
 354 if (length(predictionFile) !=0 ) {
 355
 356   predictionData <- fread(predictionFile,
 357                           na.strings = c("NA", " ", "--", "-"),
 358                           )
 359
 360   predictionData           <- as.data.frame(predictionData)
 361   rownames(predictionData) <- predictionData[, 1]
 362   predictionData[, 1]      <- NULL
 363   predictionData           <- predictionData[, colSums(is.na(predictionData)) < nrow(predictionData) * 0.5]
 364 }
 365
 366 #impute genotype values for obs with missing values,
 367 #based on mean of neighbouring 10 (arbitrary) obs
 368 genoDataMissing <- c()
 369
 370 if (sum(is.na(genoData)) > 0) {
 371   genoDataMissing<- c('yes')
 372
 373   message("sum of geno missing values, ", sum(is.na(genoData)) )
 374   genoData <- na.roughfix(genoData)
 375   genoData <- data.matrix(genoData)
 376 }
 377
 378 genoData <- genoData[order(row.names(genoData)), ]
 379
 380 #create phenotype and genotype datasets with
 381 #common stocks only
 382 message('phenotyped lines: ', length(row.names(phenoTrait)))
 383 message('genotyped lines: ', length(row.names(genoData)))
 384
 385 #extract observation lines with both
 386 #phenotype and genotype data only.
 387 commonObs <- intersect(row.names(phenoTrait), row.names(genoData))
 388 commonObs <- data.frame(commonObs)
 389 rownames(commonObs)<-commonObs[, 1]
 390
 391 message('lines with both genotype and phenotype data: ', length(row.names(commonObs)))
 392 #include in the genotype dataset only observation lines
 393 #with phenotype data
 394 message("genotype lines before filtering for phenotyped only: ", length(row.names(genoData)))
 395 genoDataFiltered <- genoData[(rownames(genoData) %in% rownames(commonObs)), ]
 396 message("genotype lines after filtering for phenotyped only: ", length(row.names(genoDataFiltered)))
 397 #drop observation lines without genotype data
 398 message("phenotype lines before filtering for genotyped only: ", length(row.names(phenoTrait)))
 399 phenoTrait <- merge(data.frame(phenoTrait), commonObs, by=0, all=FALSE)
 400 rownames(phenoTrait) <-phenoTrait[, 1]
 401 phenoTrait <- subset(phenoTrait, select=trait)
 402
 403 message("phenotype lines after filtering for genotyped only: ", length(row.names(phenoTrait)))
 404 #a set of only observation lines with genotype data
 405
 406 traitPhenoData   <- data.frame(round(phenoTrait, digits=2))
 407 phenoTrait       <- data.matrix(phenoTrait)
 408 genoDataFiltered <- data.matrix(genoDataFiltered)
 409
 410 #impute missing data in prediction data
 411 predictionDataMissing <- c()
 412 if (length(predictionData) != 0) {
 413   #purge markers unique to both populations
 414   commonMarkers    <- intersect(names(data.frame(genoDataFiltered)), names(predictionData))
 415   predictionData   <- subset(predictionData, select = commonMarkers)
 416   genoDataFiltered <- subset(genoDataFiltered, select= commonMarkers)
 417
 418  # predictionData <- data.matrix(predictionData)
 419
 420   if (sum(is.na(predictionData)) > 0) {
 421     predictionDataMissing <- c('yes')
 422     message("sum of geno missing values, ", sum(is.na(predictionData)) )
 423     predictionData <- data.matrix(na.roughfix(predictionData))
 424
 425   }
 426 }
 427
 428 relationshipMatrixFile <- grep("relationship_matrix",
 429                                outFiles,
 430                                ignore.case = TRUE,
 431                                fixed = FALSE,
 432                                value = TRUE
 433                                )
 434
 435 message("relationship matrix file: ", relationshipMatrixFile)
 436 #message("relationship matrix file size: ", file.info(relationshipMatrixFile)$size)
 437 relationshipMatrix <- c()
 438 if (length(relationshipMatrixFile) != 0) {
 439   if (file.info(relationshipMatrixFile)$size > 0 ) {
 440     relationshipDf <- as.data.frame(fread(relationshipMatrixFile))
 441
 442     rownames(relationshipDf) <- relationshipDf[, 1]
 443     relationshipDf[, 1]      <- NULL
 444     relationshipMatrix       <- data.matrix(relationshipDf)
 445   }
 446 }
 447
 448 #change genotype coding to [-1, 0, 1], to use the A.mat ) if  [0, 1, 2]
 449 genoTrCode <- grep("2", genoDataFiltered[1, ], fixed=TRUE, value=TRUE)
 450 if(length(genoTrCode) != 0) {
 451   genoDataFiltered <- genoDataFiltered - 1
 452 }
 453
 454 if (length(predictionData) != 0 ) {
 455   genoSlCode <- grep("2", predictionData[1, ], fixed=TRUE, value=TRUE)
 456   if (length(genoSlCode) != 0 ) {
 457     predictionData <- predictionData - 1
 458   }
 459 }
 460
 461 ordered.markerEffects <- c()
 462 if ( length(predictionData) == 0 ) {
 463   markerEffects <- mixed.solve(y = phenoTrait,
 464                                Z = genoDataFiltered
 465                                )
 466
 467   ordered.markerEffects <- data.matrix(markerEffects$u)
 468   ordered.markerEffects <- data.matrix(ordered.markerEffects [order (-ordered.markerEffects[, 1]), ])
 469   ordered.markerEffects <- round(ordered.markerEffects,
 470                                  digits=5
 471                                  )
 472
 473   colnames(ordered.markerEffects) <- c("Marker Effects")
 474
 475
 476 #correlation between breeding values based on
 477 #marker effects and relationship matrix
 478 #corGEBVs <- cor(genoDataMatrix %*% markerEffects$u, iGEBV$u)
 479 }
 480
 481 #additive relationship model
 482 #calculate the inner products for
 483 #genotypes (realized relationship matrix)
 484 if (length(relationshipMatrixFile) != 0) {
 485   if (file.info(relationshipMatrixFile)$size == 0) {
 486     relationshipMatrix <- tcrossprod(data.matrix(genoData))
 487   }
 488 }
 489 relationshipMatrixFiltered <- relationshipMatrix[(rownames(relationshipMatrix) %in% rownames(commonObs)),]
 490 relationshipMatrixFiltered <- relationshipMatrixFiltered[, (colnames(relationshipMatrixFiltered) %in% rownames(commonObs))]
 491
 492 #construct an identity matrix for genotypes
 493 identityMatrix <- diag(nrow(phenoTrait))
 494
 495 relationshipMatrixFiltered <- data.matrix(relationshipMatrixFiltered)
 496
 497 iGEBV <- mixed.solve(y = phenoTrait,
 498                      Z = identityMatrix,
 499                      K = relationshipMatrixFiltered
 500                      )
 501
 502 iGEBVu <- iGEBV$u
 503
 504 heritability <- c()
 505
 506 if ( is.null(predictionFile) == TRUE ) {
 507    # heritability <- round((iGEBV$Vu /(iGEBV$Vu + iGEBV$Ve) * 100), digits=3)
 508     cat("\n", file=varianceComponentsFile,  append=TRUE)
 509     cat('Error variance', iGEBV$Ve, file=varianceComponentsFile, sep="\t", append=TRUE)
 510     cat("\n", file=varianceComponentsFile,  append=TRUE)
 511     cat('Additive genetic variance',  iGEBV$Vu, file=varianceComponentsFile, sep='\t', append=TRUE)
 512     cat("\n", file=varianceComponentsFile,  append=TRUE)
 513     cat('&#956;', iGEBV$beta,file=varianceComponentsFile, sep='\t', append=TRUE)
 514     cat("\n", file=varianceComponentsFile,  append=TRUE)
 515   #  cat('Heritability (h, %)', heritability, file=varianceComponentsFile, sep='\t', append=TRUE)
 516 }
 517
 518 iGEBV <- data.matrix(iGEBVu)
 519
 520 ordered.iGEBV <- as.data.frame(iGEBV[order(-iGEBV[, 1]), ] )
 521
 522 ordered.iGEBV <- round(ordered.iGEBV,
 523                        digits = 3
 524                        )
 525
 526 combinedGebvsFile <- grep('selected_traits_gebv',
 527                           outFiles,
 528                           ignore.case = TRUE,
 529                           fixed = FALSE,
 530                           value = TRUE
 531                           )
 532
 533 allGebvs<-c()
 534 if (length(combinedGebvsFile) != 0)
 535   {
 536     fileSize <- file.info(combinedGebvsFile)$size
 537     if (fileSize != 0 )
 538       {
 539         combinedGebvs <- as.data.frame(fread(combinedGebvsFile))
 540
 541         rownames(combinedGebvs) <- combinedGebvs[,1]
 542         combinedGebvs[,1]       <- NULL
 543
 544         colnames(ordered.iGEBV) <- c(trait)
 545
 546         traitGEBV <- as.data.frame(ordered.iGEBV)
 547         allGebvs <- merge(combinedGebvs, traitGEBV,
 548                           by = 0,
 549                           all = TRUE
 550                           )
 551
 552         rownames(allGebvs) <- allGebvs[,1]
 553         allGebvs[,1] <- NULL
 554      }
 555   }
 556
 557 colnames(ordered.iGEBV) <- c(trait)
 558
 559 #cross-validation
 560 validationAll <- c()
 561
 562 if(is.null(predictionFile)) {
 563   genoNum <- nrow(phenoTrait)
 564 if(genoNum < 20 ) {
 565   warning(genoNum, " is too small number of genotypes.")
 566 }
 567 reps <- round_any(genoNum, 10, f = ceiling) %/% 10
 568
 569 genotypeGroups <-c()
 570
 571 if (genoNum %% 10 == 0) {
 572     genotypeGroups <- rep(1:10, reps)
 573   } else {
 574     genotypeGroups <- rep(1:10, reps) [- (genoNum %% 10) ]
 575   }
 576
 577 set.seed(4567)
 578 genotypeGroups <- genotypeGroups[ order (runif(genoNum)) ]
 579
 580 for (i in 1:10) {
 581   tr <- paste("trPop", i, sep = ".")
 582   sl <- paste("slPop", i, sep = ".")
 583
 584   trG <- which(genotypeGroups != i)
 585   slG <- which(genotypeGroups == i)
 586
 587   assign(tr, trG)
 588   assign(sl, slG)
 589
 590   kblup <- paste("rKblup", i, sep = ".")
 591
 592   result <- kinship.BLUP(y = phenoTrait[trG, ],
 593                          G.train = genoDataFiltered[trG, ],
 594                          G.pred = genoDataFiltered[slG, ],
 595                          mixed.method = "REML",
 596                          K.method = "RR",
 597                          )
 598
 599   assign(kblup, result)
 600
 601 #calculate cross-validation accuracy
 602   valCorData <- merge(phenoTrait[slG, ], result$g.pred, by=0, all=FALSE)
 603   rownames(valCorData) <- valCorData[, 1]
 604   valCorData[, 1]      <- NULL
 605
 606   accuracy <- try(cor(valCorData))
 607   validation <- paste("validation", i, sep = ".")
 608
 609   cvTest <- paste("Validation test", i, sep = " ")
 610
 611   if ( class(accuracy) != "try-error")
 612     {
 613       accuracy <- round(accuracy[1,2], digits = 3)
 614       accuracy <- data.matrix(accuracy)
 615
 616       colnames(accuracy) <- c("correlation")
 617       rownames(accuracy) <- cvTest
 618
 619       assign(validation, accuracy)
 620
 621       if (!is.na(accuracy[1,1])) {
 622         validationAll <- rbind(validationAll, accuracy)
 623       }
 624     }
 625 }
 626
 627 validationAll <- data.matrix(validationAll[order(-validationAll[, 1]), ])
 628
 629 if (!is.null(validationAll)) {
 630     validationMean <- data.matrix(round(colMeans(validationAll),
 631                                       digits = 2
 632                                       )
 633                                 )
 634
 635     rownames(validationMean) <- c("Average")
 636
 637     validationAll <- rbind(validationAll, validationMean)
 638     colnames(validationAll) <- c("Correlation")
 639   }
 640 }
 641 #predict GEBVs for selection population
 642 #if (length(predictionData) !=0 ) {
 643 #    predictionData <- data.matrix(round(predictionData, digits = 0 ))
 644 #  }
 645
 646 predictionPopResult <- c()
 647 predictionPopGEBVs  <- c()
 648
 649 if (length(predictionData) != 0) {
 650     message("running prediction for selection candidates...marker data", ncol(predictionData), " vs. ", ncol(genoDataFiltered))
 651
 652     predictionPopResult <- kinship.BLUP(y = phenoTrait,
 653                                         G.train = genoDataFiltered,
 654                                         G.pred = predictionData,
 655                                         mixed.method = "REML",
 656                                         K.method = "RR"
 657                                         )
 658  message("running prediction for selection candidates...DONE!!")
 659
 660     predictionPopGEBVs <- round(data.matrix(predictionPopResult$g.pred), digits = 3)
 661     predictionPopGEBVs <- data.matrix(predictionPopGEBVs[order(-predictionPopGEBVs[, 1]), ])
 662
 663     colnames(predictionPopGEBVs) <- c(trait)
 664
 665 }
 666
 667 if (!is.null(predictionPopGEBVs) & length(predictionPopGEBVsFile) != 0)  {
 668     write.table(predictionPopGEBVs,
 669                 file = predictionPopGEBVsFile,
 670                 sep = "\t",
 671                 col.names = NA,
 672                 quote = FALSE,
 673                 append = FALSE
 674                 )
 675 }
 676
 677 if(!is.null(validationAll)) {
 678     write.table(validationAll,
 679                 file = validationFile,
 680                 sep = "\t",
 681                 col.names = NA,
 682                 quote = FALSE,
 683                 append = FALSE
 684                 )
 685 }
 686
 687 if (!is.null(ordered.markerEffects)) {
 688     write.table(ordered.markerEffects,
 689                 file = markerFile,
 690                 sep = "\t",
 691                 col.names = NA,
 692                 quote = FALSE,
 693                 append = FALSE
 694                 )
 695 }
 696
 697 if (!is.null(ordered.iGEBV)) {
 698     write.table(ordered.iGEBV,
 699                 file = blupFile,
 700                 sep = "\t",
 701                 col.names = NA,
 702                 quote = FALSE,
 703                 append = FALSE
 704                 )
 705 }
 706
 707 if (length(combinedGebvsFile) != 0 ) {
 708     if(file.info(combinedGebvsFile)$size == 0) {
 709         write.table(ordered.iGEBV,
 710                     file = combinedGebvsFile,
 711                     sep = "\t",
 712                     col.names = NA,
 713                     quote = FALSE,
 714                     )
 715       } else {
 716       write.table(allGebvs,
 717                   file = combinedGebvsFile,
 718                   sep = "\t",
 719                   quote = FALSE,
 720                   col.names = NA,
 721                   )
 722     }
 723 }
 724
 725 if (!is.null(traitPhenoData) & length(traitPhenoFile) != 0) {
 726     write.table(traitPhenoData,
 727                 file = traitPhenoFile,
 728                 sep = "\t",
 729                 col.names = NA,
 730                 quote = FALSE,
 731                 )
 732 }
 733
 734
 735
 736 ## if (!is.null(genoDataMissing)) {
 737 ##   write.table(genoData,
 738 ##               file = genoFile,
 739 ##               sep = "\t",
 740 ##               col.names = NA,
 741 ##               quote = FALSE,
 742 ##             )
 743
 744 ## }
 745
 746 ## if (!is.null(predictionDataMissing)) {
 747 ##   write.table(predictionData,
 748 ##               file = predictionFile,
 749 ##               sep = "\t",
 750 ##               col.names = NA,
 751 ##               quote = FALSE,
 752 ##               )
 753 ## }
 754
 755
 756 if (file.info(relationshipMatrixFile)$size == 0) {
 757   write.table(relationshipMatrix,
 758               file = relationshipMatrixFile,
 759               sep = "\t",
 760               col.names = NA,
 761               quote = FALSE,
 762               )
 763 }
 764
 765
 766 if (file.info(formattedPhenoFile)$size == 0 && !is.null(formattedPhenoData) ) {
 767   write.table(formattedPhenoData,
 768               file = formattedPhenoFile,
 769               sep = "\t",
 770               col.names = NA,
 771               quote = FALSE,
 772               )
 773 }
 774
 775
 776 q(save = "no", runLast = FALSE)