R/histogram.r

   1  #SNOPSIS
   2
   3  #prepares trait phenotype data for histogram plotting
   4
   5
   6  #AUTHOR
   7  # Isaak Y Tecle (iyt2@cornell.edu)
   8
   9
  10 options(echo = FALSE)
  11
  12 library(plyr)
  13
  14 allArgs<-commandArgs(trailingOnly=TRUE)
  15
  16
  17 allTraitsPhenoFile <- sub('input_file=', "", allArgs[1])
  18 trait              <- sub('trait_name=', "", allArgs[2])
  19 traitPhenoFile     <- sub('output_file=', "", allArgs[3])
  20
  21 message("population phenotype file: ", allTraitsPhenoFile)
  22 message("pheno data file: ", traitPhenoFile)
  23 message("trait: ", trait)
  24
  25
  26 if (is.null(grep("phenotype_data", allTraitsPhenoFile)))
  27 {
  28   stop("Phenotype dataset missing.")
  29 }
  30
  31 if (is.null(grep("phenotype_trait", traitPhenoFile)))
  32 {
  33   stop("Output file is missing.")
  34 }
  35
  36 if (is.null(grep("trait_name", allArgs[3])))
  37 {
  38   stop("trait name is missing.")
  39 }
  40
  41 allTraitsPhenoData <- read.table(allTraitsPhenoFile,
  42                         header = TRUE,
  43                         row.names = NULL,
  44                         sep = "\t",
  45                         na.strings = c("NA", " ", "--", "-", ".", ".."),
  46                         dec = "."
  47                         )
  48
  49 selectColumns <- c("object_name", "object_id", "stock_id",  trait)
  50 traitPhenoData  <- allTraitsPhenoData[selectColumns]
  51
  52 dropColumns <- c("object_id", "stock_id")
  53 traitPhenoData <- traitPhenoData[, !(names(traitPhenoData) %in% dropColumns)]
  54
  55 if (class(traitPhenoData[, trait]) != 'numeric') {
  56   traitPhenoData[, trait] <- as.numeric(as.character(traitPhenoData[, trait]))
  57 }
  58
  59 if (!all(is.numeric(traitPhenoData[, trait]))) {
  60   traitPhenoData[, trait] <- sapply(traitPhenoData[, trait], function(x) ifelse(is.numeric(x), x, NA))
  61 }
  62
  63 if (!all(is.na(traitPhenoData[, trait]))) {
  64   traitPhenoData <- ddply(traitPhenoData,
  65                           "object_name",
  66                           colwise(mean, na.rm = TRUE)
  67                           )
  68
  69
  70   row.names(traitPhenoData) <- traitPhenoData[, 1]
  71   traitPhenoData[, 1] <- NULL
  72
  73   traitPhenoData <- round(traitPhenoData,
  74                           digits=2
  75                           )
  76 } else {
  77   traitPhenoData <- NULL
  78 }
  79
  80 write.table(traitPhenoData,
  81             file = traitPhenoFile,
  82             sep = "\t",
  83             col.names = NA,
  84             quote = FALSE,
  85             append = FALSE
  86             )
  87
  88 q(save = "no", runLast = FALSE)