Mercurial > repos > iuc > masigpro

comparison masigpro.R @ 5:401e7c9e09f2 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/masigpro commit 7dc16d5d4169a8fbfb3fad3b2212fd8b6d481f5f"
author iuc
date Sat, 27 Nov 2021 09:57:16 +0000
parents 0e5ddeb06378
children
comparison
equal deleted inserted replaced
4:21a81335cd13 5:401e7c9e09f2
4 # written by Clemens Blank. 4 # written by Clemens Blank.
5 # Thanks to Bjoern Gruening and Michael Love for their DESeq2 5 # Thanks to Bjoern Gruening and Michael Love for their DESeq2
6 # wrapper as a basis to build upon. 6 # wrapper as a basis to build upon.
7 7
8 # setup R error handling to go to stderr 8 # setup R error handling to go to stderr
9 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) 9 error_foo <- function() {
10 cat(geterrmessage(), file = stderr());
11 q("no", 1, F)
12 }
13 options(show.error.messages = F, error = error_foo)
10 14
11 # we need that to not crash galaxy with an UTF8 error on German LC settings. 15 # we need that to not crash galaxy with an UTF8 error on German LC settings.
12 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 16 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
13 17
14 suppressPackageStartupMessages({ 18 suppressPackageStartupMessages({
16 library("optparse") 20 library("optparse")
17 library("mclust") 21 library("mclust")
18 }) 22 })
19 23
20 # The following code fixes an error in the stepback function 24 # The following code fixes an error in the stepback function
21 # of the maSigPro package. This code is hopefully temporary 25 # of the maSigPro package. This code is hopefully temporary
22 # and can be removed if the fix is included in a future 26 # and can be removed if the fix is included in a future
23 # version. The stepback function in the maSigPro namespace 27 # version. The stepback function in the maSigPro namespace
24 # will be overwritten by the following function. 28 # will be overwritten by the following function.
25 stepback <- function (y = y, d = d, alfa = 0.05, family = gaussian() , epsilon=0.00001) 29 stepback <- function(y = y, d = d, alfa = 0.05, family = gaussian(), epsilon=0.00001) {
26 { 30 lm1 <- glm(y ~ ., data = d, family = family, epsilon = epsilon)
27 lm1 <- glm(y ~ ., data = d, family=family, epsilon=epsilon)
28 result <- summary(lm1) 31 result <- summary(lm1)
29 max <- max(result$coefficients[, 4][-1], na.rm = TRUE) 32 max <- max(result$coefficients[, 4][-1], na.rm = TRUE)
30 if (length(result$coefficients[, 4][-1]) == 1) { 33 if (length(result$coefficients[, 4][-1]) == 1) {
31 if (max > alfa) { 34 if (max > alfa) {
32 max = 0 35 max <- 0
33 lm1 <- glm(y ~ 1, family=family, epsilon=epsilon) 36 lm1 <- glm(y ~ 1, family = family, epsilon = epsilon)
34 } 37 }
35 } 38 }
36 while (max > alfa) { 39 while (max > alfa) {
37 varout <- names(result$coefficients[, 4][-1])[result$coefficients[, 40 varout <- names(result$coefficients[, 4][-1])[result$coefficients[, 4][-1] == max][1]
38 4][-1] == max][1]
39 pos <- position(matrix = d, vari = varout) 41 pos <- position(matrix = d, vari = varout)
40 d <- d[, -pos] 42 d <- d[, -pos]
41 if (length(result$coefficients[, 4][-1]) == 2) { 43 if (length(result$coefficients[, 4][-1]) == 2) {
42 min <- min(result$coefficients[, 4][-1], na.rm = TRUE) 44 min <- min(result$coefficients[, 4][-1], na.rm = TRUE)
43 lastname <- names(result$coefficients[, 4][-1])[result$coefficients[,4][-1] == min] 45 lastname <- names(result$coefficients[, 4][-1])[result$coefficients[, 4][-1] == min]
44 } 46 }
45 if (is.null(dim(d))) { 47 if (is.null(dim(d))) {
46 d <- as.data.frame(d) 48 d <- as.data.frame(d)
47 colnames(d) <- lastname 49 colnames(d) <- lastname
48 } 50 }
49 lm1 <- glm(y ~ ., data = d, family=family, epsilon=epsilon) 51 lm1 <- glm(y ~ ., data = d, family = family, epsilon = epsilon)
50 result <- summary(lm1) 52 result <- summary(lm1)
51 max <- max(result$coefficients[, 4][-1], na.rm = TRUE) 53 max <- max(result$coefficients[, 4][-1], na.rm = TRUE)
52 if (length(result$coefficients[, 4][-1]) == 1) { 54 if (length(result$coefficients[, 4][-1]) == 1) {
53 max <- result$coefficients[, 4][-1] 55 max <- result$coefficients[, 4][-1]
54 if (max > alfa) { 56 if (max > alfa) {
55 max = 0 57 max <- 0
56 lm1 <- glm(y ~ 1, family=family, epsilon=epsilon) 58 lm1 <- glm(y ~ 1, family = family, epsilon = epsilon)
57 } 59 }
58 } 60 }
59 } 61 }
60 return(lm1) 62 return(lm1)
61 } 63 }
62 64
63 unlockBinding("stepback", as.environment("package:maSigPro")) 65 unlockBinding("stepback", as.environment("package:maSigPro"))
64 assignInNamespace("stepback", stepback, ns="maSigPro", envir=as.environment("package:maSigPro")) 66 assignInNamespace("stepback", stepback, ns = "maSigPro", envir = as.environment("package:maSigPro"))
65 assign("stepback", stepback, as.environment("package:maSigPro")) 67 assign("stepback", stepback, as.environment("package:maSigPro"))
66 lockBinding("stepback", as.environment("package:maSigPro")) 68 lockBinding("stepback", as.environment("package:maSigPro"))
67 # End of temporary code to fix stepback.R 69 # End of temporary code to fix stepback.R
68 70
69 options(stringAsFactors = FALSE, useFancyQuotes = FALSE) 71 options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
72 # specify our desired options in a list 74 # specify our desired options in a list
73 # by default OptionParser will add an help option equivalent to 75 # by default OptionParser will add an help option equivalent to
74 # make_option(c("-h", "--help"), action="store_true", default=FALSE, 76 # make_option(c("-h", "--help"), action="store_true", default=FALSE,
75 # help="Show this help message and exit") 77 # help="Show this help message and exit")
76 option_list <- list( 78 option_list <- list(
77 make_option(c("-q", "--quiet"), action="store_false", 79 make_option(c("-q", "--quiet"), action = "store_false",
78 dest="verbose", help="Print little output"), 80 dest = "verbose", help = "Print little output"),
79 make_option(c("-e", "--edesign"), type="character"), 81 make_option(c("-e", "--edesign"), type = "character"),
80 make_option(c("-d", "--data"), type="character"), 82 make_option(c("-d", "--data"), type = "character"),
81 make_option(c("-o", "--outfile"), type="character"), 83 make_option(c("-o", "--outfile"), type = "character"),
82 make_option("--degree", type="integer", default=1), 84 make_option("--degree", type = "integer", default = 1),
83 make_option("--time_col", type="integer", default=1), 85 make_option("--time_col", type = "integer", default = 1),
84 make_option("--repl_col", type="integer", default=2), 86 make_option("--repl_col", type = "integer", default = 2),
85 make_option("--qvalue", type="double", default=0.05), 87 make_option("--qvalue", type = "double", default = 0.05),
86 make_option("--min_obs", type="integer", default=6), 88 make_option("--min_obs", type = "integer", default = 6),
87 make_option("--step_method", type="character", default="backward"), 89 make_option("--step_method", type = "character", default = "backward"),
88 make_option("--nvar_correction", type="logical", default=FALSE), 90 make_option("--nvar_correction", type = "logical", default = FALSE),
89 make_option("--alfa", type="double", default=0.05), 91 make_option("--alfa", type = "double", default = 0.05),
90 make_option("--rsq", type="double", default=0.7), 92 make_option("--rsq", type = "double", default = 0.7),
91 make_option("--vars", type="character", default="groups"), 93 make_option("--vars", type = "character", default = "groups"),
92 make_option("--significant_intercept", type="character", default="dummy"), 94 make_option("--significant_intercept", type = "character", default = "dummy"),
93 make_option("--cluster_data", type="integer", default=1), 95 make_option("--cluster_data", type = "integer", default = 1),
94 make_option(c("-k", "--k"), type="integer", default=9), 96 make_option(c("-k", "--k"), type = "integer", default = 9),
95 make_option("--print_cluster", type="logical", default=FALSE), 97 make_option("--print_cluster", type = "logical", default = FALSE),
96 make_option("--cluster_method", type="character", default="hclust"), 98 make_option("--cluster_method", type = "character", default = "hclust"),
97 make_option("--distance", type="character", default="cor"), 99 make_option("--distance", type = "character", default = "cor"),
98 make_option("--agglo_method", type="character", default="ward.D"), 100 make_option("--agglo_method", type = "character", default = "ward.D"),
99 make_option("--iter_max", type="integer", default=500), 101 make_option("--iter_max", type = "integer", default = 500),
100 make_option("--color_mode", type="character", default="rainbow"), 102 make_option("--color_mode", type = "character", default = "rainbow"),
101 make_option("--show_fit", type="logical", default=TRUE), 103 make_option("--show_fit", type = "logical", default = TRUE),
102 make_option("--show_lines", type="logical", default=TRUE), 104 make_option("--show_lines", type = "logical", default = TRUE),
103 make_option("--cexlab", type="double", default=0.8), 105 make_option("--cexlab", type = "double", default = 0.8),
104 make_option("--legend", type="logical", default=TRUE) 106 make_option("--legend", type = "logical", default = TRUE)
105 ) 107 )
106 108
107 # get command line options, if help option encountered print help and exit, 109 # get command line options, if help option encountered print help and exit,
108 # otherwise if options not found on command line then set defaults 110 # otherwise if options not found on command line then set defaults
109 opt <- parse_args(OptionParser(option_list=option_list)) 111 opt <- parse_args(OptionParser(option_list = option_list))
110 112
111 # enforce the following required arguments 113 # enforce the following required arguments
112 if (is.null(opt$edesign)) { 114 if (is.null(opt$edesign)) {
113 cat("'edesign' is required\n") 115 cat("'edesign' is required\n")
114 q(status=1) 116 q(status = 1)
115 } 117 }
116 if (is.null(opt$data)) { 118 if (is.null(opt$data)) {
117 cat("'data' is required\n") 119 cat("'data' is required\n")
118 q(status=1) 120 q(status = 1)
119 } 121 }
120 if (is.null(opt$outfile)) { 122 if (is.null(opt$outfile)) {
121 cat("'outfile' is required\n") 123 cat("'outfile' is required\n")
122 q(status=1) 124 q(status = 1)
123 } 125 }
124 126
125 verbose <- if (is.null(opt$quiet)) { 127 verbose <- if (is.null(opt$quiet)) {
126 TRUE 128 TRUE
127 } else { 129 } else {
128 FALSE 130 FALSE
129 } 131 }
130 132
131 edesign <- as.matrix(read.table(opt$edesign, header=TRUE, row.names = 1)) 133 edesign <- as.matrix(read.table(opt$edesign, header = TRUE, row.names = 1))
132 134
133 data <- read.table(opt$data, header=TRUE, check.names=FALSE) 135 data <- read.table(opt$data, header = TRUE, check.names = FALSE)
134 136
135 results <- maSigPro(data, edesign, degree = opt$degree, time.col = opt$time_col, 137 results <- maSigPro(data, edesign, degree = opt$degree, time.col = opt$time_col,
136 repl.col = opt$repl_col, Q = opt$qvalue, min.obs = opt$min_obs, 138 repl.col = opt$repl_col, Q = opt$qvalue, min.obs = opt$min_obs,
137 step.method = opt$step_method, nvar.correction = opt$nvar_correction, 139 step.method = opt$step_method, nvar.correction = opt$nvar_correction,
138 alfa = opt$alfa, rsq = opt$rsq, vars = opt$vars, 140 alfa = opt$alfa, rsq = opt$rsq, vars = opt$vars,
143 color.mode = opt$color_mode, show.fit = opt$show_fit, 145 color.mode = opt$color_mode, show.fit = opt$show_fit,
144 show.lines = opt$show_lines, cexlab = opt$cexlab, 146 show.lines = opt$show_lines, cexlab = opt$cexlab,
145 legend = opt$legend) 147 legend = opt$legend)
146 148
147 if (opt$print_cluster) { 149 if (opt$print_cluster) {
148 for (i in 1:length(results$sig.genes)) { 150 for (i in seq_len(length(results$sig.genes))) {
149 151 colname <- paste(names(results$sig.genes)[i], "cluster", sep = "_")
150 colname <- paste(names(results$sig.genes)[i], "cluster", sep = "_") 152 results$summary[colname] <- ""
151 153 results$summary[[colname]][seq_len(length(results$sig.genes[[i]]$sig.profiles$`cluster$cut`))] <-
152 results$summary[colname] <- "" 154 results$sig.genes[[i]]$sig.profiles$`cluster$cut`
153 results$summary[[colname]][1:length(results$sig.genes[[i]]$sig.profiles$`cluster$cut`)] <-
154 results$sig.genes[[i]]$sig.profiles$`cluster$cut`
155 } 155 }
156 } 156 }
157 157
158 filename <- opt$outfile 158 filename <- opt$outfile
159 159
160 write.table((results$summary), file=filename, sep="\t", quote=FALSE, 160 write.table((results$summary), file = filename, sep = "\t", quote = FALSE,
161 row.names=FALSE, col.names=TRUE) 161 row.names = FALSE, col.names = TRUE)
162 162
163 cat("Session information:\n\n") 163 cat("Session information:\n\n")
164 164
165 sessionInfo() 165 sessionInfo()