R/solGS/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(phenoAnalysis)
  13 library(dplyr)
  14
  15 allArgs     <- commandArgs()
  16
  17 outputFiles <- scan(grep("output_files", allArgs, value = TRUE),
  18                     what = "character")
  19
  20 inputFiles  <- scan(grep("input_files", allArgs, value = TRUE),
  21                    what = "character")
  22
  23 allTraitsPhenoFile <- grep("phenotype_data", inputFiles, value = TRUE)
  24 message('pheno file: ', allTraitsPhenoFile)
  25
  26 metaDataFile <- grep("metadata_file", inputFiles, value = TRUE)
  27 message('metadata  file: ', metaDataFile)
  28
  29 traitsFile <- grep("traits", inputFiles, value = TRUE)
  30 message('traits file: ', traitsFile)
  31 traits  <- scan(traitsFile,  what = "character")
  32 trait  <- strsplit(traits, "\t")
  33 message("trait: ", trait)
  34
  35 traitPhenoMeansFile <- grep("phenotype_data", outputFiles, value = TRUE)
  36 message('pheno file: ', traitPhenoMeansFile)
  37
  38 traitRawPhenoFile <- grep("trait_raw_phenodata", outputFiles, value = TRUE)
  39 message('raw pheno file: ', traitRawPhenoFile)
  40
  41
  42 if (is.null(grep("phenotype_data", allTraitsPhenoFile))) {
  43   stop("Phenotype dataset missing.")
  44 }
  45
  46 if (is.null(grep("phenotype_trait", traitPhenoMeansFile))) {
  47   stop("Output file is missing.")
  48 }
  49
  50 if (is.null(trait)) {
  51   stop("trait name is missing.")
  52 }
  53
  54 traitRawPhenoData <- extractPhenotypes(inputFiles, metadata_file)
  55
  56 allTraitsPhenoData <- read.table(allTraitsPhenoFile,
  57                                  header = TRUE,
  58                                  row.names = NULL,
  59                                  sep = "\t",
  60                                  na.strings = c("NA", " ", "--", "-", ".", ".."),
  61                                  )
  62
  63 allTraitsPhenoData <- data.frame(fread(allTraitsPhenoFile,
  64                                        na.strings = c("NA", "", "--", "-", ".")
  65                                        ))
  66
  67 traitPhenoMeansData <- getAdjMeans(allTraitsPhenoData,
  68                                traitName=trait,
  69                                calcAverages=TRUE)
  70
  71 keepMetaCols <- c('observationUnitName', 'germplasmName', 'studyDbId', 'locationName',
  72                                               'studyYear', 'replicate', 'blockNumber')
  73
  74 traitRawPhenoData <- allTraitsPhenoData %>%
  75                                     select(c(keepMetaCols, trait))
  76
  77
  78 write.table(traitPhenoMeansData,
  79             file = traitPhenoMeansFile,
  80             sep = "\t",
  81             quote = FALSE,
  82             )
  83
  84 write.table(traitRawPhenoData,
  85                 file = traitRawPhenoFile,
  86                 sep = "\t",
  87                 na = 'NA',
  88                 quote = FALSE,
  89                 )
  90
  91 q(save = "no", runLast = FALSE)