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
  73   diagnosticsFile  <- grep("anova_diagnostics",
  74                         outputFiles,
  75                         value = TRUE)
  76
  77   anovaOut   <- runAnova(phenoData, trait)
  78
  79   png(diagnosticsFile, 960, 480)
  80   par(mfrow=c(1,2))
  81   plot(fitted(anovaOut), resid(anovaOut),
  82        xlab="Fitted values",
  83        ylab="Residuals",
  84        main="Fitted values vs Residuals")
  85   abline(0,0)
  86   qqnorm(resid(anovaOut))
  87   dev.off()
  88
  89   anovaTable <- getAnovaTable(anovaOut,
  90                               tableType="html",
  91                               traitName=trait,
  92                               out=anovaHtmlFile)
  93
  94   anovaTable <- getAnovaTable(anovaOut,
  95                               tableType="text",
  96                               traitName=trait,
  97                               out=anovaTxtFile)
  98
  99   adjMeans   <- getAdjMeans(phenoData, trait)
 100
 101   fwrite(adjMeans,
 102          file      = adjMeansFile,
 103          row.names = FALSE,
 104          sep       = "\t",
 105          quote     = FALSE,
 106          )
 107
 108   sink(modelSummFile)
 109   print(anovaOut)
 110   sink()
 111
 112 }
 113
 114
 115 q(save = "no", runLast = FALSE)