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