R/anova.r

   1  #SNOPSIS
   2
   3  #runs ANOVA.
   4
   5
   6  #AUTHOR
   7  # Isaak Y Tecle (iyt2@cornell.edu)
   8
   9
  10 options(echo = FALSE)
  11
  12
  13 #library(dplyr)
  14 library(data.table)
  15 library(phenoAnalysis)
  16
  17
  18 allArgs     <- commandArgs()
  19
  20 outputFiles <- scan(grep("output_files", allArgs, value = TRUE),
  21                     what = "character")
  22
  23 inputFiles  <- scan(grep("input_files", allArgs, value = TRUE),
  24                    what = "character")
  25
  26 phenoDataFile <- grep("phenotype_data", inputFiles, value = TRUE)
  27 message('pheno file: ', phenoDataFile)
  28 traitsFile <- grep("traits", inputFiles, value = TRUE)
  29 message('traits file: ', traitsFile)
  30
  31 #phenoDataFile <- '/export/prod/tmp/localhost/GBSApeKIgenotypingv4/solgs/cache/phenotype_data_139.txt';
  32 #phenoDataFile <- '/mnt/hgfs/cxgn/phenotype_data_RCBD_1596.txt';
  33
  34 phenoData <- fread(phenoDataFile,
  35                      na.strings=c("NA", "-", " ", ".", ".."))
  36
  37 phenoData <- data.frame(phenoData)
  38
  39 traits  <- scan(traitsFile,
  40                 what = "character",)
  41
  42 traits  <- strsplit(traits, "\t")
  43
  44 #needs more work for multi traits anova
  45 for (trait in traits) {
  46   anovaFiles     <- grep("anova_table",
  47                          outputFiles,
  48                          value = TRUE)
  49
  50   message('anova file: ', anovaFiles)
  51   anovaHtmlFile  <- grep("html",
  52                          anovaFiles,
  53                          value = TRUE)
  54
  55   message('anova html file: ', anovaHtmlFile)
  56   anovaTxtFile   <- grep("txt",
  57                          anovaFiles,
  58                          value = TRUE)
  59
  60   message('anova txt file: ', anovaTxtFile)
  61   modelSummFile <- grep("anova_model",
  62                         outputFiles,
  63                         value = TRUE)
  64
  65   message('model file: ', modelSummFile)
  66   adjMeansFile  <- grep("adj_means",
  67                         outputFiles,
  68                         value = TRUE)
  69
  70   message('means file: ', adjMeansFile)
  71
  72   anovaOut   <- runAnova(phenoData, trait)
  73
  74   anovaTable <- getAnovaTable(anovaOut,
  75                               tableType="html",
  76                               traitName=trait,
  77                               out=anovaHtmlFile)
  78
  79   anovaTable <- getAnovaTable(anovaOut,
  80                               tableType="text",
  81                               traitName=trait,
  82                               out=anovaTxtFile)
  83
  84   adjMeans   <- getAdjMeans(phenoData, trait)
  85
  86   fwrite(adjMeans,
  87          file      = adjMeansFile,
  88          row.names = FALSE,
  89          sep       = "\t",
  90          quote     = FALSE,
  91          )
  92
  93   sink(modelSummFile)
  94   print(anovaOut)
  95   sink()
  96 }
  97
  98
  99 q(save = "no", runLast = FALSE)