Mercurial > repos > marie-tremblay-metatoul > asca_w4m
diff galaxy/asca_wrapper.R @ 2:826a2a145147 draft default tip
Uploaded
author | marie-tremblay-metatoul |
---|---|
date | Thu, 09 Aug 2018 04:25:34 -0400 |
parents | 20395c0079ae |
children |
line wrap: on
line diff
--- a/galaxy/asca_wrapper.R Mon Jul 30 07:47:12 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,225 +0,0 @@ -#!/usr/bin/env Rscript - -################################################################################################### -# -# MetStaT ASCA.calculate function -# -# -# R-Package: MetStaT -# -# Version: 1.0 -# -# Author (asca.calculate): Tim Dorscheidt -# Author (wrapper & .r adaptation for workflow4metabolomics.org): M. Tremblay-Franco & Y. Guitton # -# -# Expected parameters from the commandline -# input files: -# dataMatrix -# sampleMetadata -# variableMetadata -# params: -# Factors (Factor1 & Factor2) -# scaling -# Number of permutations -# Significance threshold -# output files: -# sampleMetadata -# variableMetadata -# Graphical outputs -# Information text -################################################################################################### -pkgs=c("MetStaT","batch","pcaMethods") -for(pkg in pkgs) { - suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE))) - cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="") -} - - -listArguments = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects - -#Redirect all stdout to the log file -sink(listArguments$information) - -# ----- PACKAGE ----- -cat("\tPACKAGE INFO\n") -sessionInfo() - -source_local <- function(fname) { - argv <- commandArgs(trailingOnly = FALSE) - base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) - source(paste(base_dir, fname, sep="/")) -} - -#load asca_w4m function -source_local("asca_w4m.R") -source_local("ASCA.Calculate_w4m.R") -source_local("ASCA.PlotScoresPerLevel_w4m.R") -print("first loadings OK") - -## libraries -##---------- - -cat('\n\nRunning ASCA.calculate\n'); -options(warn=-1); -#remove rgl warning -options(rgl.useNULL = TRUE); - - -## constants -##---------- - -modNamC <- "asca" ## module name - -topEnvC <- environment() -flgC <- "\n" - -## functions -##----------For manual input of function -##--end function - -flgF <- function(tesC, - envC = topEnvC, - txtC = NA) { ## management of warning and error messages - - tesL <- eval(parse(text = tesC), envir = envC) - - if(!tesL) { - - sink(NULL) - stpTxtC <- ifelse(is.na(txtC), - paste0(tesC, " is FALSE"), - txtC) - - stop(stpTxtC, - call. = FALSE) - - } - -} ## flgF - - -## log file -##--------- -cat("\nStart of the '", modNamC, "' Galaxy module call: ", - format(Sys.time(), "%a %d %b %Y %X"), "\n", sep="") - - -## arguments -##---------- -## loading files and checks -xMN <- t(as.matrix(read.table(listArguments[["dataMatrix_in"]], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t"))) - -samDF <- read.table(listArguments[["sampleMetadata_in"]], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t") - -varDF <- read.table(listArguments[["variableMetadata_in"]], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t") - -result <- asca_w4m(xMN, samDF, c(listArguments[["factor1"]],listArguments[["factor2"]]), varDF, as.numeric(listArguments[["threshold"]]), - scaling=listArguments[["scaling"]], listArguments[["nPerm"]]) - - -##saving - -if (exists("result")) { - ## writing output files - cat("\n\nWriting output files\n\n"); - write.table(result[[4]], - file = listArguments$sampleMetadata_out, - quote = FALSE, - row.names = TRUE, - sep = "\t") - - write.table(result[[5]], - file = listArguments$variableMetadata_out, - quote = FALSE, - row.names = TRUE, - sep = "\t") - - # Graphical display for each significant parameter - print(result[[3]]) - cat("\n p-value of Residuals must not be taken into account\n") - - if (any(result[[2]] < as.numeric(listArguments[["threshold"]]))) - { - data.asca.permutation <- result[[2]] - design <- data.matrix(samDF[, colnames(samDF) %in% c(listArguments[["factor1"]],listArguments[["factor2"]])]) - - pdf(listArguments$figure, onefile=TRUE) - par(mfrow=c(1,3)) - if (data.asca.permutation[1] < as.numeric(listArguments[["threshold"]])) - { - eigenvalues <- data.frame(1:length(unique(design[,1])), result[[1]]$'1'$svd$var.explained[1:length(unique(design[,1]))]) - colnames(eigenvalues) <- c("PC", "explainedVariance") - barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component") - noms <- levels(as.factor(samDF[, listArguments$factor1])) - ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="1", interaction=0, factorName=listArguments$factor1, factorModalite=noms) - f1.loadings <- data.matrix(result[[5]][,2:3]) - f1.loadings.leverage <- diag(f1.loadings%*%t(f1.loadings)) - names(f1.loadings.leverage) <- colnames(xMN) - f1.loadings.leverage <- sort(f1.loadings.leverage, decreasing=TRUE) - barplot(f1.loadings.leverage[f1.loadings.leverage > 0.001], main="PC1 loadings") - } - if (data.asca.permutation[2] < as.numeric(listArguments[["threshold"]])) - { - eigenvalues <- data.frame(1:length(unique(design[,2])), result[[1]]$'2'$svd$var.explained[1:length(unique(design[,2]))]) - colnames(eigenvalues) <- c("PC", "explainedVariance") - barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component") - noms <- levels(as.factor(samDF[, listArguments$factor2])) - ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="2", interaction=0, factorName=listArguments$factor2, factorModalite=noms) - f2.loadings <- data.matrix(result[[5]][,4:5]) - f2.loadings.leverage <- diag(f2.loadings%*%t(f2.loadings)) - names(f2.loadings.leverage) <- colnames(xMN) - f2.loadings.leverage <- sort(f2.loadings.leverage, decreasing=TRUE) - barplot(f2.loadings.leverage[f2.loadings.leverage > 0.001], main="PC1 loadings") - } - if (data.asca.permutation[3] < as.numeric(listArguments[["threshold"]])) - { - eigenvalues <- data.frame(1:(length(unique(design[,1]))*length(unique(design[,2]))), result[[1]]$'12'$svd$var.explained[1:(length(unique(design[,1]))*length(unique(design[,2])))]) - colnames(eigenvalues) <- c("PC", "explainedVariance") - barplot(eigenvalues[,2], names.arg=eigenvalues[,1], ylab="% of explained variance", xlab="Principal component") - noms1 <- data.matrix(levels(as.factor(samDF[, listArguments$factor1]))) - noms2 <- data.matrix(levels(as.factor(samDF[, listArguments$factor2]))) - noms <- apply(noms1, 1, FUN=function(x){paste(x, "-", noms2, sep="")}) - noms <- apply(noms, 1, FUN=function(x){c(noms)}) - ASCA.PlotScoresPerLevel_w4m(result[[1]], ee="12", interaction=1, factorModalite=noms[,1]) - f1f2.loadings <- data.matrix(result[[5]][,6:7]) - f1f2.loadings.leverage <- diag(f1f2.loadings%*%t(f1f2.loadings)) - names(f1f2.loadings.leverage) <- colnames(xMN) - f1f2.loadings.leverage <- sort(f1f2.loadings.leverage, decreasing=TRUE) - barplot(f1f2.loadings.leverage[f1f2.loadings.leverage > 0.001], main="PC1 loadings") - } - dev.off() - } - - tryCatch({ - save(result, file="asca.RData"); - }, warning = function(w) { - print(paste("Warning: ", w)); - }, error = function(err) { - stop(paste("ERROR saving result RData object:", err)); - }); -} - -## ending -##------- - -cat("\nEnd of the '", modNamC, "' Galaxy module call: ", - format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "") - -sink() - -# options(stringsAsFactors = strAsFacL) - - -rm(list = ls())