R/anova.r

   1  #SNOPSIS
   2
   3  #runs ANOVA.
   4
   5  #AUTHOR
   6  # Isaak Y Tecle (iyt2@cornell.edu)
   7
   8
   9 options(echo = FALSE)
  10
  11
  12 #library(dplyr)
  13 library(data.table)
  14 library(phenoAnalysis)
  15 library(methods)
  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
  29 traitsFile <- grep("traits", inputFiles, value = TRUE)
  30 message('traits file: ', traitsFile)
  31
  32
  33 phenoData <- fread(phenoDataFile,
  34                    na.strings=c("NA", "-", " ", ".", ".."))
  35
  36 phenoData <- data.frame(phenoData)
  37
  38 traits  <- scan(traitsFile,
  39                 what = "character")
  40
  41 traits  <- strsplit(traits, "\t")
  42
  43
  44 #needs more work for multi traits anova
  45 for (trait in traits) {
  46
  47     message('trait: ', trait)
  48     anovaFiles     <- grep("anova_table",
  49                            outputFiles,
  50                            value = TRUE)
  51
  52     message('anova file: ', anovaFiles)
  53     anovaHtmlFile  <- grep("html",
  54                            anovaFiles,
  55                            value = TRUE)
  56
  57     message('anova html file: ', anovaHtmlFile)
  58     anovaTxtFile   <- grep("txt",
  59                            anovaFiles,
  60                            value = TRUE)
  61
  62     message('anova txt file: ', anovaTxtFile)
  63     modelSummFile <- grep("anova_model",
  64                           outputFiles,
  65                           value = TRUE)
  66
  67     message('model file: ', modelSummFile)
  68     adjMeansFile  <- grep("adj_means",
  69                           outputFiles,
  70                           value = TRUE)
  71
  72     message('means file: ', adjMeansFile)
  73
  74
  75     diagnosticsFile  <- grep("anova_diagnostics",
  76                              outputFiles,
  77                              value = TRUE)
  78
  79     errorFile  <- grep("anova_error",
  80                        outputFiles,
  81                        value = TRUE)
  82
  83     anovaOut <- runAnova(phenoData, trait)
  84
  85     if (class(anovaOut)[1] == 'merModLmerTest') {
  86
  87         png(diagnosticsFile, 960, 480)
  88         par(mfrow=c(1,2))
  89         plot(fitted(anovaOut), resid(anovaOut),
  90              xlab="Fitted values",
  91              ylab="Residuals",
  92              main="Fitted values vs Residuals")
  93         abline(0,0)
  94         qqnorm(resid(anovaOut))
  95         dev.off()
  96
  97         anovaTable <- getAnovaTable(anovaOut,
  98                                     tableType="html",
  99                                     traitName=trait,
 100                                     out=anovaHtmlFile)
 101
 102         anovaTable <- getAnovaTable(anovaOut,
 103                                     tableType="text",
 104                                     traitName=trait,
 105                                     out=anovaTxtFile)
 106
 107         adjMeans   <- getAdjMeans(phenoData, trait)
 108
 109         fwrite(adjMeans,
 110                file      = adjMeansFile,
 111                row.names = FALSE,
 112                sep       = "\t",
 113                quote     = FALSE,
 114                )
 115
 116         sink(modelSummFile)
 117         print(anovaOut)
 118         sink()
 119     } else {
 120         cat(anovaOut, file=errorFile)
 121     }
 122
 123 }
 124
 125
 126 q(save = "no", runLast = FALSE)