R/solGS/selection_index.r

   1 # calculates selection index
   2 # and ranks genotypes accordingly
   3 # Isaak Y Tecle iyt2cornell.edu
   4
   5
   6 options(echo = FALSE)
   7
   8 library(data.table)
   9 library(stats)
  10 library(stringi)
  11 library(dplyr)
  12
  13 allArgs <- commandArgs()
  14
  15 inputFiles <- scan(grep("input_files", allArgs, value = TRUE),
  16                    what = "character")
  17
  18 relWeightsFile <- grep("rel_weights", inputFiles, value = TRUE)
  19
  20 outputFiles <- scan(grep("output_files", allArgs, value = TRUE),
  21                     what = "character")
  22
  23 traitsFiles <- grep("gebv_files_of_traits", inputFiles, value = TRUE)
  24
  25 gebvsSelectionIndexFile <- grep("gebvs_selection_index",
  26                                 outputFiles,
  27                                 value = TRUE)
  28
  29 selectionIndexFile <- grep("selection_index_only",
  30                            outputFiles,
  31                            value=TRUE)
  32
  33 inTraitFiles   <- scan(traitsFiles, what = "character")
  34
  35 traitFilesList <- strsplit(inTraitFiles, "\t");
  36 traitsTotal    <- length(traitFilesList)
  37
  38 if (traitsTotal == 0)
  39   stop("There are no traits with GEBV data.")
  40 if (length(relWeightsFile) == 0)
  41   stop("There is no file with relative weights of traits.")
  42
  43
  44 relWeights           <- data.frame(fread(relWeightsFile, header = TRUE))
  45 rownames(relWeights) <- relWeights[, 1]
  46 relWeights[, 1]      <- NULL
  47
  48 if (is.null(relWeights)) {
  49     stop('There were no relative weights for all the traits.')
  50 }
  51
  52 combinedRelGebvs <- c()
  53
  54 for (i in 1:traitsTotal) {
  55   traitFile           <- traitFilesList[[i]]
  56   traitGEBV           <- data.frame(fread(traitFile, header = TRUE))
  57   rownames(traitGEBV) <- traitGEBV[, 1]
  58   traitGEBV[, 1]      <- NULL
  59   traitGEBV           <- traitGEBV[order(rownames(traitGEBV)),,drop=FALSE]
  60   trait               <- colnames(traitGEBV)
  61
  62   relWeight <- relWeights[trait, ]
  63
  64   if (is.na(relWeight) == FALSE && relWeight != 0 ) {
  65
  66       weightedTraitGEBV <- apply(traitGEBV, 1,
  67                                  function(x) x*relWeight)
  68
  69       weightedTraitGEBV <- data.frame(weightedTraitGEBV)
  70       colnames(weightedTraitGEBV) <- paste0(trait, '_weighted')
  71
  72       combinedRelGebvs  <- merge(combinedRelGebvs, weightedTraitGEBV,
  73                                  by = 0,
  74                                  all = TRUE)
  75
  76
  77       rownames(combinedRelGebvs) <- combinedRelGebvs[, 1]
  78       combinedRelGebvs[, 1]      <- NULL
  79     }
  80 }
  81
  82 sumRelWeights <- apply(relWeights, 2, sum)
  83 sumRelWeights <- sumRelWeights[[1]]
  84
  85 combinedRelGebvs$Index <- apply(combinedRelGebvs, 1, function (x) sum(x))
  86
  87 combinedRelGebvs <- combinedRelGebvs[ with(combinedRelGebvs,
  88                                            order(-combinedRelGebvs$Index)
  89                                            ),
  90                                      ]
  91
  92 combinedRelGebvs <- round(combinedRelGebvs, 2)
  93
  94 selectionIndex <-c()
  95
  96 if (!is.null(combinedRelGebvs)) {
  97   selectionIndex <- subset(combinedRelGebvs,
  98                            select = 'Index'
  99                            )
 100 }
 101
 102 if (gebvsSelectionIndexFile != 0) {
 103   if (!is.null(combinedRelGebvs)) {
 104     fwrite(combinedRelGebvs,
 105            file      = gebvsSelectionIndexFile,
 106            sep       = "\t",
 107            row.names = TRUE,
 108            quote     = FALSE,
 109            )
 110       }
 111 }
 112
 113 if (!is.null(selectionIndexFile)) {
 114   if (!is.null(selectionIndex)) {
 115     fwrite(selectionIndex,
 116            file      = selectionIndexFile,
 117            row.names = TRUE,
 118            quote     = FALSE,
 119            sep       = "\t",
 120            )
 121   }
 122 }
 123
 124 q(save = "no", runLast = FALSE)