diff univariate_script.R @ 0:ab2ee3414e4e draft

planemo upload for repository https://github.com/workflow4metabolomics/univariate.git commit 98e8f4464b2f7321acb010e26e2a1c82fe37096e

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/univariate_script.R	Tue Oct 24 08:57:25 2017 -0400
@@ -0,0 +1,332 @@
+univariateF <- function(datMN,
+                        samDF,
+                        varDF,
+                        facC,
+                        tesC = c("ttest", "wilcoxon", "anova", "kruskal", "pearson", "spearman")[1],
+                        adjC = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")[7],
+                        thrN = 0.05,
+                        pdfC) {
+
+
+    ## Option
+
+    strAsFacL <- options()$stringsAsFactors
+    options(stingsAsFactors = FALSE)
+    options(warn = -1)
+
+    ## Getting the response (either a factor or a numeric)
+
+    if(mode(samDF[, facC]) == "character") {
+        facFcVn <- factor(samDF[, facC])
+        facLevVc <- levels(facFcVn)
+    } else
+        facFcVn <- samDF[, facC]
+
+    cat("\nPerforming '", tesC, "'\n", sep="")
+
+    varPfxC <- paste0(make.names(facC), "_", tesC, "_")
+
+    
+    if(tesC %in% c("ttest", "wilcoxon", "pearson", "spearman")) {
+
+        
+        switch(tesC,
+               ttest = {
+                   staF <- function(y) diff(tapply(y, facFcVn, function(x) mean(x, na.rm = TRUE)))
+                   tesF <- function(y) t.test(y ~ facFcVn)[["p.value"]]
+               },
+               wilcoxon = {
+                   staF <- function(y) diff(tapply(y, facFcVn, function(x) median(x, na.rm = TRUE)))
+                   tesF <- function(y) wilcox.test(y ~ facFcVn)[["p.value"]]
+               },
+               pearson = {
+                   staF <- function(y) cor(facFcVn, y, method = "pearson", use = "pairwise.complete.obs")
+                   tesF <- function(y) cor.test(facFcVn, y, method = "pearson", use = "pairwise.complete.obs")[["p.value"]]
+               },
+               spearman = {
+                   staF <- function(y) cor(facFcVn, y, method = "spearman", use = "pairwise.complete.obs")
+                   tesF <- function(y) cor.test(facFcVn, y, method = "spearman", use = "pairwise.complete.obs")[["p.value"]]
+               })
+
+        staVn <- apply(datMN, 2, staF)
+
+        adjVn <- p.adjust(apply(datMN,
+                                2,
+                                tesF),
+                          method = adjC)
+
+        sigVn <- as.numeric(adjVn < thrN)
+
+        if(tesC %in% c("ttest", "wilcoxon"))
+            varPfxC <- paste0(varPfxC, paste(rev(facLevVc), collapse = "."), "_")
+
+        varDF[, paste0(varPfxC, ifelse(tesC %in% c("ttest", "wilcoxon"), "dif", "cor"))] <- staVn
+
+        varDF[, paste0(varPfxC, adjC)] <- adjVn
+
+        varDF[, paste0(varPfxC, "sig")] <- sigVn
+
+        ## graphic
+
+        pdf(pdfC, onefile = TRUE)
+
+        varVi <- which(sigVn > 0)
+
+        if(tesC %in% c("ttest", "wilcoxon")) {
+
+            facVc <- as.character(facFcVn)
+            names(facVc) <- rownames(samDF)
+
+            for(varI in varVi) {
+                
+                varC <- rownames(varDF)[varI]
+                
+                boxF(facFcVn,
+                     datMN[, varI],
+                     paste0(varC, " (", adjC, " = ", signif(adjVn[varI], 2), ")"),
+                     facVc)
+                
+            }
+
+        } else { ## pearson or spearman
+
+            for(varI in varVi) {
+
+                varC <- rownames(varDF)[varI]
+
+                mod <- lm(datMN[, varI] ~  facFcVn)
+
+                plot(facFcVn, datMN[, varI],
+                     xlab = facC,
+                     ylab = "",
+                     pch = 18,
+                     main = paste0(varC, " (", adjC, " = ", signif(adjVn[varI], 2), ", R2 = ", signif(summary(mod)$r.squared, 2), ")"))
+            
+                abline(mod, col = "red")
+
+                }
+
+        }
+
+        dev.off()
+
+        
+    } else if(tesC == "anova") {
+
+        
+        ## getting the names of the pairwise comparisons 'class1Vclass2'
+        prwVc <- rownames(TukeyHSD(aov(datMN[, 1] ~ facFcVn))[["facFcVn"]])
+
+        prwVc <- gsub("-", ".", prwVc, fixed = TRUE) ## 2016-08-05: '-' character in dataframe column names seems not to be converted to "." by write.table on ubuntu R-3.3.1
+
+        ## omnibus and post-hoc tests
+        
+        aovMN <- t(apply(datMN, 2, function(varVn) {
+
+            aovMod <- aov(varVn ~ facFcVn)
+            pvaN <- summary(aovMod)[[1]][1, "Pr(>F)"]
+            hsdMN <- TukeyHSD(aovMod)[["facFcVn"]]
+            c(pvaN, c(hsdMN[, c("diff", "p adj")]))
+
+        }))
+
+        difVi <- 1:length(prwVc) + 1
+
+        ## difference of the means for each pairwise comparison
+        
+        difMN <- aovMN[, difVi]
+        colnames(difMN) <- paste0(varPfxC, prwVc, "_dif")
+
+        ## correction for multiple testing
+        
+        aovMN <- aovMN[, -difVi, drop = FALSE]
+        aovMN <- apply(aovMN, 2, function(pvaVn) p.adjust(pvaVn, method = adjC))
+
+        ## significance coding (0 = not significant, 1 = significant)
+        
+        adjVn <- aovMN[, 1]
+        sigVn <-  as.numeric(adjVn < thrN)
+
+        aovMN <- aovMN[, -1, drop = FALSE]
+        colnames(aovMN) <- paste0(varPfxC, prwVc, "_", adjC)
+
+        aovSigMN <- aovMN < thrN
+        mode(aovSigMN) <- "numeric"
+        colnames(aovSigMN) <- paste0(varPfxC, prwVc, "_sig")
+
+        ## final aggregated table
+
+        resMN <- cbind(adjVn, sigVn, difMN, aovMN, aovSigMN)
+        colnames(resMN)[1:2] <- paste0(varPfxC, c(adjC, "sig"))
+
+        varDF <- cbind.data.frame(varDF, as.data.frame(resMN))
+
+        ## graphic
+
+        pdf(pdfC, onefile = TRUE)
+        
+        for(varI in 1:nrow(varDF)) {
+            
+            if(sum(aovSigMN[varI, ]) > 0) {
+                
+                varC <- rownames(varDF)[varI]
+
+                boxplot(datMN[, varI] ~ facFcVn,
+                        main = paste0(varC, " (", adjC, " = ", signif(adjVn[varI], 2), ")"))
+                
+                for(prwI in 1:length(prwVc)) {
+                    
+                    if(aovSigMN[varI, paste0(varPfxC, prwVc[prwI], "_sig")] == 1) {
+                        
+                        claVc <- unlist(strsplit(prwVc[prwI], ".", fixed = TRUE))
+                        aovClaVl <- facFcVn %in% claVc
+                        aovFc <- facFcVn[aovClaVl, drop = TRUE]
+                        aovVc <- as.character(aovFc)
+                        names(aovVc) <- rownames(samDF)[aovClaVl]
+                        boxF(aovFc,
+                             datMN[aovClaVl, varI],
+                             paste0(varC, " (", adjC, " = ", signif(aovMN[varI, paste0(varPfxC, prwVc[prwI], "_", adjC)], 2), ")"),
+                             aovVc)
+                        
+                    }
+                       
+                }
+                
+            }
+            
+        }
+
+        dev.off()
+
+        
+    } else if(tesC == "kruskal") {
+        
+
+        ## getting the names of the pairwise comparisons 'class1.class2'
+        
+        nemMN <- posthoc.kruskal.nemenyi.test(datMN[, 1], facFcVn, "Tukey")[["p.value"]]
+        nemVl <- c(lower.tri(nemMN, diag = TRUE))
+        nemClaMC <- cbind(rownames(nemMN)[c(row(nemMN))][nemVl],
+                          colnames(nemMN)[c(col(nemMN))][nemVl])
+        nemNamVc <- paste0(nemClaMC[, 1], ".", nemClaMC[, 2])
+        pfxNemVc <- paste0(varPfxC, nemNamVc)
+
+        ## omnibus and post-hoc tests
+        
+        nemMN <- t(apply(datMN, 2, function(varVn) {
+
+            pvaN <- kruskal.test(varVn ~ facFcVn)[["p.value"]]
+            varNemMN <- posthoc.kruskal.nemenyi.test(varVn, facFcVn, "Tukey")[["p.value"]]
+            c(pvaN, c(varNemMN))
+
+        }))
+
+        ## correction for multiple testing
+        
+        nemMN <- apply(nemMN, 2,
+                       function(pvaVn) p.adjust(pvaVn, method = adjC))
+        adjVn <- nemMN[, 1]
+        sigVn <- as.numeric(adjVn < thrN)
+        nemMN <- nemMN[, c(FALSE, nemVl)]
+        colnames(nemMN) <- paste0(pfxNemVc, "_", adjC)
+
+        ## significance coding (0 = not significant, 1 = significant)
+        
+        nemSigMN <- nemMN < thrN
+        mode(nemSigMN) <- "numeric"
+        colnames(nemSigMN) <- paste0(pfxNemVc, "_sig")
+
+        ## difference of the medians for each pairwise comparison
+        
+        difMN <- sapply(1:nrow(nemClaMC), function(prwI) {
+            prwVc <- nemClaMC[prwI, ]
+            prwVi <- which(facFcVn %in% prwVc)
+            prwFacFc <- factor(as.character(facFcVn)[prwVi], levels = prwVc)
+            apply(datMN[prwVi, ], 2, function(varVn) -diff(as.numeric(tapply(varVn, prwFacFc, function(x) median(x, na.rm = TRUE)))))
+        })
+        colnames(difMN) <- gsub("_sig", "_dif", colnames(nemSigMN))
+
+        ## final aggregated table
+        
+        resMN <- cbind(adjVn, sigVn, difMN, nemMN, nemSigMN)
+        colnames(resMN)[1:2] <- paste0(varPfxC, c(adjC, "sig"))
+
+        varDF <- cbind.data.frame(varDF, as.data.frame(resMN))
+
+        ## graphic
+
+        pdf(pdfC, onefile = TRUE)
+        
+        for(varI in 1:nrow(varDF)) {
+            
+            if(sum(nemSigMN[varI, ]) > 0) {
+                
+                varC <- rownames(varDF)[varI]
+
+                boxplot(datMN[, varI] ~ facFcVn,
+                        main = paste0(varC, " (", adjC, " = ", signif(adjVn[varI], 2), ")"))
+                
+                for(nemI in 1:length(nemNamVc)) {
+                    
+                    if(nemSigMN[varI, paste0(varPfxC, nemNamVc[nemI], "_sig")] == 1) {
+                        
+                        nemClaVc <- nemClaMC[nemI, ]
+                        nemClaVl <- facFcVn %in% nemClaVc
+                        nemFc <- facFcVn[nemClaVl, drop = TRUE]
+                        nemVc <- as.character(nemFc)
+                        names(nemVc) <- rownames(samDF)[nemClaVl]
+                        boxF(nemFc,
+                             datMN[nemClaVl, varI],
+                             paste0(varC, " (", adjC, " = ", signif(nemMN[varI, paste0(varPfxC, nemNamVc[nemI], "_", adjC)], 2), ")"),
+                             nemVc)
+                        
+                    }
+                       
+                }
+                
+            }
+            
+        }
+
+        dev.off()
+        
+    }
+    
+    names(sigVn) <- rownames(varDF)
+    sigSumN <- sum(sigVn, na.rm = TRUE)
+    if(sigSumN) {
+        cat("\nThe following ", sigSumN, " variable", ifelse(sigSumN > 1, "s", ""), " (", round(sigSumN / length(sigVn) * 100), "%) ", ifelse(sigSumN > 1, "were", "was"), " found significant at the ", thrN, " level:\n", sep = "")
+        cat(paste(rownames(varDF)[sigVn > 0], collapse = "\n"), "\n", sep = "")
+    } else
+        cat("\nNo significant variable found at the selected ", thrN, " level\n", sep = "")
+
+    options(stingsAsFactors = strAsFacL)
+
+    return(varDF)
+
+}
+
+
+boxF <- function(xFc,
+                 yVn,
+                 maiC,
+                 xVc) {
+    
+    boxLs <- boxplot(yVn ~  xFc,
+                     main = maiC)
+    
+    outVn <- boxLs[["out"]]
+    
+    if(length(outVn)) {
+        
+        for(outI in 1:length(outVn)) {
+            levI <- which(levels(xFc) == xVc[names(outVn)[outI]])
+            text(levI,
+                 outVn[outI],
+                 labels = names(outVn)[outI],
+                 pos = ifelse(levI == 2, 2, 4))
+            }
+        
+    }
+    
+}