R/combine_populations.r

   1 #formats and combines phenotype (of a single trait)
   2 #and genotype datasets of multiple
   3 #populations
   4
   5 options(echo = FALSE)
   6
   7 library(stats)
   8 library(stringr)
   9 library(randomForest)
  10 library(plyr)
  11 library(lme4)
  12
  13
  14 allArgs <- commandArgs()
  15
  16 inFile <- grep("input_files",
  17                allArgs,
  18                ignore.case = TRUE,
  19                perl = TRUE,
  20                value = TRUE
  21                )
  22
  23 outFile <- grep("output_files",
  24                 allArgs,
  25                 ignore.case = TRUE,
  26                 perl = TRUE,
  27                 value = TRUE
  28                 )
  29
  30 outFiles <- scan(outFile,
  31                  what = "character"
  32                  )
  33
  34 combinedGenoFile <- grep("genotype_data",
  35                          outFiles,
  36                          ignore.case = TRUE,
  37                          fixed = FALSE,
  38                          value = TRUE
  39                          )
  40
  41 combinedPhenoFile <- grep("phenotype_data",
  42                           outFiles,
  43                           ignore.case = TRUE,
  44                           fixed = FALSE,
  45                           value = TRUE
  46                           )
  47
  48 inFiles <- scan(inFile,
  49                 what = "character"
  50                 )
  51 print(inFiles)
  52
  53 traitFile <- grep("trait_",
  54                   inFiles,
  55                   ignore.case = TRUE,
  56                   fixed = FALSE,
  57                   value = TRUE
  58                   )
  59
  60 trait <- scan(traitFile,
  61               what = "character",
  62               )
  63
  64 traitInfo<-strsplit(trait, "\t");
  65 traitId<-traitInfo[[1]]
  66 traitName<-traitInfo[[2]]
  67
  68 #extract trait phenotype data from all populations
  69 #and combine them into one dataset
  70
  71 allPhenoFiles <- grep("phenotype_data",
  72                   inFiles,
  73                   ignore.case = TRUE,
  74                   fixed = FALSE,
  75                   value = TRUE
  76                   )
  77 message("phenotype files: ", allPhenoFiles)
  78
  79 allGenoFiles <- grep("genotype_data",
  80                   inFiles,
  81                   ignore.case = TRUE,
  82                   fixed = FALSE,
  83                   value = TRUE
  84                   )
  85
  86 popsPhenoSize     <- length(allPhenoFiles)
  87 popsGenoSize      <- length(allGenoFiles)
  88 popIds            <- c()
  89 combinedPhenoPops <- c()
  90
  91 for (popPhenoNum in 1:popsPhenoSize)
  92   {
  93     popId <- str_extract(allPhenoFiles[[popPhenoNum]], "\\d+")
  94     popIds <- append(popIds, popId)
  95
  96     phenoData <- read.table(allPhenoFiles[[popPhenoNum]],
  97                             header = TRUE,
  98                             row.names = 1,
  99                             sep = "\t",
 100                             na.strings = c("NA", " ", "--", "-", "."),
 101                             dec = "."
 102                            )
 103
 104
 105     phenoTrait <- subset(phenoData,
 106                          select = c("object_name", "object_id", "design", "block", "replicate", traitName)
 107                          )
 108
 109     experimentalDesign <- phenoTrait[2, 'design']
 110
 111     if (is.na(experimentalDesign) == TRUE) {experimentalDesign <- c('No Design')}
 112
 113     if ((experimentalDesign == 'Augmented' || experimentalDesign == 'RCBD')  &&  unique(phenoTrait$block) > 1) {
 114
 115       message("experimental design: ", experimentalDesign)
 116
 117       augData <- subset(phenoTrait,
 118                         select = c("object_name", "object_id",  "block",  traitName)
 119                         )
 120
 121       colnames(augData)[1] <- "genotypes"
 122       colnames(augData)[4] <- "trait"
 123
 124       model <- try(lmer(trait ~ 0 + genotypes + (1|block),
 125                         augData,
 126                         na.action = na.omit
 127                         ))
 128
 129       if (class(model) != "try-error") {
 130         phenoTrait <- data.frame(fixef(model))
 131
 132         colnames(phenoTrait) <- traitName
 133
 134         nn <- gsub('genotypes', '', rownames(phenoTrait))
 135         rownames(phenoTrait) <- nn
 136
 137         phenoTrait <- round(phenoTrait, digits = 2)
 138
 139         formattedPhenoData[, traitName] <- phenoTrait
 140       }
 141
 142     } else if (experimentalDesign == 'Alpha') {
 143
 144       alphaData <-  phenoTrait
 145
 146       colnames(alphaData)[1] <- "genotypes"
 147       colnames(alphaData)[5] <- "trait"
 148
 149       model <- try(lmer(trait ~ 0 + genotypes + (1|replicate/block),
 150                         alphaData,
 151                         na.action = na.omit
 152                         ))
 153
 154       if (class(model) != "try-error") {
 155         phenoTrait <- data.frame(fixef(model))
 156
 157         colnames(phenoTrait) <- traitName
 158
 159         nn <- gsub('genotypes', '', rownames(phenotrait))
 160         rownames(phenoTrait) <- nn
 161
 162         phenoTrait <- round(phenoTrait, digits = 2)
 163
 164         formattedPhenoData[, i] <- phenoTrait
 165       }
 166
 167   } else {
 168
 169     phenoTrait <- subset(phenoData,
 170                          select = c("object_name", "stock_id", traitName)
 171                          )
 172
 173     if (sum(is.na(phenoTrait)) > 0) {
 174       message("No. of pheno missing values: ", sum(is.na(phenoTrait)))
 175
 176       phenoTrait <- na.omit(phenoTrait)
 177
 178       #calculate mean of reps/plots of the same accession and
 179       #create new df with the accession means
 180       phenoTrait$stock_id <- NULL
 181       phenoTrait   <- phenoTrait[order(row.names(phenoTrait)), ]
 182
 183       print('phenotyped lines before averaging')
 184       print(length(row.names(phenoTrait)))
 185
 186       phenoTrait<-ddply(phenoTrait, "object_name", colwise(mean))
 187
 188       print('phenotyped lines after averaging')
 189       print(length(row.names(phenoTrait)))
 190
 191       row.names(phenoTrait) <- phenoTrait[, 1]
 192       phenoTrait[, 1] <- NULL
 193
 194       phenoTrait <- round(phenoTrait, digits = 2)
 195
 196     } else {
 197       print ('No missing data')
 198       phenoTrait$stock_id <- NULL
 199       phenoTrait   <- phenoTrait[order(row.names(phenoTrait)), ]
 200
 201       print('phenotyped lines before averaging')
 202       print(length(row.names(phenoTrait)))
 203
 204       phenoTrait<-ddply(phenoTrait, "object_name", colwise(mean))
 205
 206       print('phenotyped lines after averaging')
 207       print(length(row.names(phenoTrait)))
 208
 209       row.names(phenoTrait) <- phenoTrait[, 1]
 210       phenoTrait[, 1] <- NULL
 211
 212       phenoTrait <- round(phenoTrait, digits = 2)
 213
 214     }
 215   }
 216     newTraitName = paste(traitName, popId, sep = "_")
 217     colnames(phenoTrait)[1] <- newTraitName
 218
 219     if (popPhenoNum == 1 )
 220       {
 221         print('no need to combine, yet')
 222         combinedPhenoPops <- phenoTrait
 223
 224       } else {
 225       print('combining...')
 226       combinedPhenoPops <- merge(combinedPhenoPops, phenoTrait,
 227                             by = 0,
 228                             all=TRUE,
 229                             )
 230
 231       rownames(combinedPhenoPops) <- combinedPhenoPops[, 1]
 232       combinedPhenoPops$Row.names <- NULL
 233
 234     }
 235 }
 236
 237 #fill in missing data in combined phenotype dataset
 238 #using row means
 239 naIndices <- which(is.na(combinedPhenoPops), arr.ind=TRUE)
 240 combinedPhenoPops <- as.matrix(combinedPhenoPops)
 241 combinedPhenoPops[naIndices] <- rowMeans(combinedPhenoPops, na.rm=TRUE)[naIndices[,1]]
 242 combinedPhenoPops <- as.data.frame(combinedPhenoPops)
 243
 244 message("combined total number of stocks in phenotype dataset (before averaging): ", length(rownames(combinedPhenoPops)))
 245
 246 combinedPhenoPops$Average<-round(apply(combinedPhenoPops,
 247                                        1,
 248                                        function(x)
 249                                        { mean(x) }
 250                                        ),
 251                                  digits = 2
 252                                  )
 253
 254 markersList      <- c()
 255 combinedGenoPops <- c()
 256
 257 for (popGenoNum in 1:popsGenoSize)
 258   {
 259     popId <- str_extract(allGenoFiles[[popGenoNum]], "\\d+")
 260     popIds <- append(popIds, popId)
 261
 262     genoData <- read.table(allGenoFiles[[popGenoNum]],
 263                             header = TRUE,
 264                             row.names = 1,
 265                             sep = "\t",
 266                             na.strings = c("NA", " ", "--", "-"),
 267                             dec = "."
 268                            )
 269
 270     popMarkers <- colnames(genoData)
 271     message("No of markers from population ", popId, ": ", length(popMarkers))
 272     #print(popMarkers)
 273
 274     if (sum(is.na(genoData)) > 0)
 275       {
 276         message("sum of geno missing values: ", sum(is.na(genoData)))
 277         genoData <- na.roughfix(genoData)
 278         message("total number of stocks for pop ", popId,": ", length(rownames(genoData)))
 279       }
 280
 281     if (popGenoNum == 1 )
 282       {
 283         print('no need to combine, yet')
 284         combinedGenoPops <- genoData
 285
 286       } else {
 287         print('combining genotype datasets...')
 288         combinedGenoPops <-rbind(combinedGenoPops, genoData)
 289       }
 290
 291
 292   }
 293 message("combined total number of stocks in genotype dataset: ", length(rownames(combinedGenoPops)))
 294 #discard duplicate clones
 295 combinedGenoPops <- unique(combinedGenoPops)
 296 message("combined unique number of stocks in genotype dataset: ", length(rownames(combinedGenoPops)))
 297
 298 message("writing data into files...")
 299 if(length(combinedPhenoFile) != 0 )
 300   {
 301       write.table(combinedPhenoPops,
 302                   file = combinedPhenoFile,
 303                   sep = "\t",
 304                   quote = FALSE,
 305                   col.names = NA,
 306                   )
 307   }
 308
 309 if(length(combinedGenoFile) != 0 )
 310   {
 311       write.table(combinedGenoPops,
 312                   file = combinedGenoFile,
 313                   sep = "\t",
 314                   quote = FALSE,
 315                   col.names = NA,
 316                   )
 317   }
 318
 319 q(save = "no", runLast = FALSE)