R/phenotypic_correlation.r

   1  #SNOPSIS
   2
   3  #runs phenotypic correlation analysis.
   4  #Correlation coeffiecients are stored in tabular and json formats
   5
   6  #AUTHOR
   7  # Isaak Y Tecle (iyt2@cornell.edu)
   8
   9
  10 options(echo = FALSE)
  11
  12 library(gplots)
  13 library(ltm)
  14 library(plyr)
  15 library(rjson)
  16 library(lme4)
  17 library(data.table)
  18 #library(rbenchmark)
  19
  20
  21 allargs<-commandArgs()
  22
  23 refererQtl <- grep("qtl",
  24                    allargs,
  25                    ignore.case=TRUE,
  26                    perl=TRUE,
  27                    value=TRUE
  28                    )
  29
  30 phenoDataFile <- grep("\\/phenotype_data",
  31                       allargs,
  32                       ignore.case=TRUE,
  33                       perl=TRUE,
  34                       value=TRUE
  35                       )
  36
  37 correCoefficientsFile <- grep("corre_coefficients_table",
  38                               allargs,
  39                               ignore.case=TRUE,
  40                               perl=TRUE,
  41                               value=TRUE
  42                               )
  43
  44 correCoefficientsJsonFile <- grep("corre_coefficients_json",
  45                                   allargs,
  46                                   ignore.case=TRUE,
  47                                   perl=TRUE,
  48                                   value=TRUE
  49                                   )
  50
  51
  52 formattedPhenoFile <- grep("formatted_phenotype_data",
  53                            allargs,
  54                            ignore.case = TRUE,
  55                            fixed = FALSE,
  56                            value = TRUE
  57                            )
  58
  59
  60 formattedPhenoData <- c()
  61 phenoData          <- c()
  62
  63 if (file.info(formattedPhenoFile)$size > 0 ) {
  64
  65   formattedPhenoData <- as.data.frame(fread(formattedPhenoFile,
  66                                    na.strings = c("NA", " ", "--", "-", "."),
  67                                    ))
  68
  69   row.names(formattedPhenoData) <- formattedPhenoData[, 1]
  70   formattedPhenoData[, 1]       <- NULL
  71
  72 } else {
  73
  74   if ( length(refererQtl) != 0 ) {
  75
  76     phenoData <- as.data.frame(fread(phenoDataFile,
  77                                      sep=",",
  78                                      na.strings=c("NA", "-", " ", ".", "..")
  79                                      ))
  80   } else {
  81
  82     phenoData <- as.data.frame(fread(phenoDataFile,
  83                                      na.strings = c("NA", " ", "--", "-", ".", "..")
  84                                      ))
  85   }
  86 }
  87
  88 allTraitNames <- c()
  89 nonTraitNames <- c()
  90
  91 if (length(refererQtl) != 0) {
  92
  93   allNames      <- names(phenoData)
  94   nonTraitNames <- c("ID")
  95   allTraitNames <- allNames[! allNames %in% nonTraitNames]
  96
  97 } else if (file.info(formattedPhenoFile)$size == 0 && length(refererQtl) == 0) {
  98
  99   dropColumns <- c("uniquename", "stock_name")
 100   phenoData   <- phenoData[,!(names(phenoData) %in% dropColumns)]
 101
 102   allNames      <- names(phenoData)
 103   nonTraitNames <- c("object_name", "object_id", "stock_id", "design", "block", "replicate")
 104   allTraitNames <- allNames[! allNames %in% nonTraitNames]
 105
 106 }
 107
 108 if (!is.null(phenoData) && length(refererQtl) == 0) {
 109
 110   for (i in allTraitNames) {
 111
 112     if (class(phenoData[, i]) != 'numeric') {
 113       phenoData[, i] <- as.numeric(as.character(phenoData[, i]))
 114     }
 115
 116     if (all(is.nan(phenoData$i))) {
 117       phenoData[, i] <- sapply(phenoData[, i], function(x) ifelse(is.numeric(x), x, NA))
 118     }
 119   }
 120 }
 121
 122 phenoData     <- phenoData[, colSums(is.na(phenoData)) < nrow(phenoData)]
 123 allTraitNames <- names(phenoData)[! names(phenoData) %in% nonTraitNames]
 124
 125 ###############################
 126 if (length(refererQtl) == 0  ) {
 127   if (file.info(formattedPhenoFile)$size == 0) {
 128
 129     cnt   <- 0
 130
 131     for (i in allTraitNames) {
 132
 133       cnt   <- cnt + 1
 134       trait <- i
 135
 136       experimentalDesign <- c()
 137
 138       if ('design' %in% colnames(phenoData)) {
 139
 140         experimentalDesign <- phenoData[2, 'design']
 141
 142         if (is.na(experimentalDesign)) {
 143           experimentalDesign <- c('No Design')
 144         }
 145
 146       } else {
 147         experimentalDesign <- c('No Design')
 148       }
 149
 150       if ((experimentalDesign == 'Augmented' || experimentalDesign == 'RCBD')  &&  unique(phenoData$block) > 1) {
 151
 152       message("GS experimental design: ", experimentalDesign)
 153
 154       augData <- subset(phenoData,
 155                         select = c("object_name", "object_id",  "block",  trait)
 156                         )
 157
 158       colnames(augData)[1] <- "genotypes"
 159       colnames(augData)[4] <- "trait"
 160
 161       model <- try(lmer(trait ~ 0 + genotypes + (1|block),
 162                         augData,
 163                         na.action = na.omit
 164                         ))
 165       genoEffects <- c()
 166
 167       if (class(model) != "try-error") {
 168         genoEffects <- data.frame(fixef(model))
 169
 170         colnames(genoEffects) <- trait
 171
 172         nn <- gsub('genotypes', '', rownames(genoEffects))
 173         rownames(genoEffects) <- nn
 174
 175         genoEffects <- round(genoEffects, digits = 2)
 176       }
 177
 178       if (cnt == 1 ) {
 179         formattedPhenoData <- data.frame(genoEffects)
 180       } else {
 181         formattedPhenoData <-  merge(formattedPhenoData, genoEffects, by=0, all=TRUE)
 182         row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 183         formattedPhenoData[, 1] <- NULL
 184       }
 185
 186     } else if (experimentalDesign == 'Alpha') {
 187       trait <- i
 188
 189       message("Experimental desgin: ", experimentalDesign)
 190
 191       alphaData <- subset(phenoData,
 192                             select = c("object_name", "object_id","block", "replicate", trait)
 193                             )
 194
 195       colnames(alphaData)[1] <- "genotypes"
 196       colnames(alphaData)[5] <- "trait"
 197
 198       model <- try(lmer(trait ~ 0 + genotypes + (1|replicate/block),
 199                         alphaData,
 200                         na.action = na.omit
 201                         ))
 202
 203       if (class(model) != "try-error") {
 204         genoEffects <- data.frame(fixef(model))
 205
 206         colnames(genoEffects) <- trait
 207
 208         nn <- gsub('genotypes', '', rownames(genoEffects))
 209         rownames(genoEffects) <- nn
 210
 211         genoEffects <- round(genoEffects, digits = 2)
 212
 213       }
 214
 215       if(cnt == 1 ) {
 216         formattedPhenoData <- genoEffects
 217       } else {
 218         formattedPhenoData <-  merge(formattedPhenoData, genoEffects, by=0, all=TRUE)
 219         row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 220         formattedPhenoData[, 1] <- NULL
 221       }
 222
 223     } else {
 224       message("GS experimental design: ", experimentalDesign)
 225
 226       dropColumns <- c("object_id", "stock_id", "design",  "block", "replicate")
 227
 228       formattedPhenoData <- phenoData[, !(names(phenoData) %in% dropColumns)]
 229
 230       formattedPhenoData <- ddply(formattedPhenoData,
 231                                   "object_name",
 232                                   colwise(mean, na.rm=TRUE)
 233                                   )
 234
 235       row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 236       formattedPhenoData[, 1] <- NULL
 237
 238     }
 239     }
 240   }
 241 } else {
 242   message("qtl stuff")
 243   formattedPhenoData <- ddply(phenoData,
 244                               "ID",
 245                               colwise(mean, na.rm=TRUE)
 246                               )
 247
 248   row.names(formattedPhenoData) <- formattedPhenoData[, 1]
 249   formattedPhenoData[, 1] <- NULL
 250
 251   formattedPhenoData <- round(formattedPhenoData, digits = 2)
 252 }
 253
 254 coefpvalues <- rcor.test(formattedPhenoData,
 255                          method="pearson",
 256                          use="pairwise"
 257                          )
 258
 259 coefficients <- coefpvalues$cor.mat
 260 allcordata   <- coefpvalues$cor.mat
 261
 262 allcordata[lower.tri(allcordata)] <- coefpvalues$p.values[, 3]
 263 diag(allcordata) <- 1.00
 264
 265 pvalues <- as.matrix(allcordata)
 266
 267 pvalues <- round(pvalues,
 268                  digits=2
 269                  )
 270
 271 coefficients <- round(coefficients,
 272                       digits=3
 273                       )
 274
 275 allcordata   <- round(allcordata,
 276                     digits=3
 277                     )
 278
 279 #remove rows and columns that are all "NA"
 280 if ( apply(coefficients,
 281            1,
 282            function(x)any(is.na(x))
 283            )
 284     ||
 285     apply(coefficients,
 286           2,
 287           function(x)any(is.na(x))
 288           )
 289     )
 290   {
 291
 292     coefficients<-coefficients[-which(apply(coefficients,
 293                                             1,
 294                                             function(x)all(is.na(x)))
 295                                       ),
 296                                -which(apply(coefficients,
 297                                             2,
 298                                             function(x)all(is.na(x)))
 299                                       )
 300                                ]
 301   }
 302
 303
 304 pvalues[upper.tri(pvalues)]           <- NA
 305 coefficients[upper.tri(coefficients)] <- NA
 306
 307 coefficients2json <- function(mat) {
 308   mat <- as.list(as.data.frame(t(mat)))
 309   names(mat) <- NULL
 310   toJSON(mat)
 311 }
 312
 313 traits <- colnames(coefficients)
 314
 315 correlationList <- list(
 316                      "traits" = toJSON(traits),
 317                      "coefficients " =coefficients2json(coefficients)
 318                    )
 319
 320 correlationJson <- paste("{",paste("\"", names(correlationList), "\":", correlationList, collapse=","), "}")
 321
 322 write.table(coefficients,
 323             file=correCoefficientsFile,
 324             col.names=TRUE,
 325             row.names=TRUE,
 326             quote=FALSE,
 327             dec="."
 328             )
 329
 330 write.table(correlationJson,
 331             file=correCoefficientsJsonFile,
 332             col.names=FALSE,
 333             row.names=FALSE,
 334             )
 335
 336
 337 if (file.info(formattedPhenoFile)$size == 0 && !is.null(formattedPhenoData) ) {
 338   write.table(formattedPhenoData,
 339               file = formattedPhenoFile,
 340               sep = "\t",
 341               col.names = NA,
 342               quote = FALSE,
 343               )
 344 }
 345
 346
 347 q(save = "no", runLast = FALSE)