Mercurial > repos > iuc > sleuth

comparison sleuth.R @ 3:2c3d294dbe42 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/sleuth commit 4adb19064a6973121680119529444286c24c1ac1
author iuc
date Mon, 02 Jun 2025 21:31:22 +0000
parents d6b5fc94062c
children
comparison
equal deleted inserted replaced
2:d6b5fc94062c 3:2c3d294dbe42
1 library(sleuth, 1 library(sleuth,
2 quietly = TRUE, 2 quietly = TRUE,
3 warn.conflicts = FALSE) 3 warn.conflicts = FALSE
4 )
4 library(annotables, quietly = TRUE, warn.conflicts = FALSE) 5 library(annotables, quietly = TRUE, warn.conflicts = FALSE)
5 library(argparse, quietly = TRUE, warn.conflicts = FALSE) 6 library(argparse, quietly = TRUE, warn.conflicts = FALSE)
6 library(tidyverse) 7 library(tidyverse)
7 8
8 9
9 # setup R error handling to go to stderr 10 # setup R error handling to go to stderr
10 options( 11 options(
11 show.error.messages = FALSE, 12 show.error.messages = FALSE,
12 error = function() { 13 error = function() {
13 cat(geterrmessage(), file = stderr()) 14 cat(geterrmessage(), file = stderr())
14 q("no", 1, FALSE) 15 q("no", 1, FALSE)
15 } 16 }
16 ) 17 )
17 18
18 # we need that to not crash galaxy with an UTF8 error on German LC settings. 19 # we need that to not crash galaxy with an UTF8 error on German LC settings.
19 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 20 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
20 21
26 # Collect arguments from command line 27 # Collect arguments from command line
27 parser <- ArgumentParser(description = "Sleuth R script") 28 parser <- ArgumentParser(description = "Sleuth R script")
28 29
29 parser$add_argument("--factorLevel", action = "append", required = FALSE) 30 parser$add_argument("--factorLevel", action = "append", required = FALSE)
30 parser$add_argument("--factorLevel_counts", 31 parser$add_argument("--factorLevel_counts",
31 action = "append", 32 action = "append",
32 required = FALSE) 33 required = FALSE
33 parser$add_argument("--factorLevel_n", action = "append", required = FALSE) 34 )
34 parser$add_argument("--cores", type = "integer", required = FALSE) 35 parser$add_argument("--factorLevel_n", action = "append", required = FALSE)
36 parser$add_argument("--cores", type = "integer", required = FALSE)
35 parser$add_argument("--normalize", action = "store_true", required = FALSE) 37 parser$add_argument("--normalize", action = "store_true", required = FALSE)
36 parser$add_argument("--nbins", type = "integer", required = FALSE) 38 parser$add_argument("--nbins", type = "integer", required = FALSE)
37 parser$add_argument("--lwr", type = "numeric", required = FALSE) 39 parser$add_argument("--lwr", type = "numeric", required = FALSE)
38 parser$add_argument("--upr", type = "numeric", required = FALSE) 40 parser$add_argument("--upr", type = "numeric", required = FALSE)
39 parser$add_argument("--metadata_file", 41 parser$add_argument("--metadata_file",
40 action = "append", 42 action = "append",
41 required = FALSE) 43 required = FALSE
44 )
42 parser$add_argument("--experiment_design", required = FALSE) 45 parser$add_argument("--experiment_design", required = FALSE)
43 46
44 args <- parser$parse_args() 47 args <- parser$parse_args()
45 48
46 if (args$experiment_design == "complex") { 49 if (args$experiment_design == "complex") {
47 ## Complex experiment design 50 ## Complex experiment design
48 ############################ 51 ############################
49 52
50 s2c <- 53 s2c <-
51 read.table(file = args$metadata_file, 54 read.table(
52 header = TRUE, 55 file = args$metadata_file,
53 sep = "\t") 56 header = TRUE,
57 sep = "\t"
58 )
54 59
55 s2c$path <- file.path("./kallisto_outputs/", paste(s2c$path, ".h5", sep = "")) 60 s2c$path <- file.path("./kallisto_outputs/", paste(s2c$path, ".h5", sep = ""))
56 for (f in args$factorLevel_counts) { 61 for (f in args$factorLevel_counts) {
57 file.rename(f, paste(f, ".h5", sep = "")) 62 file.rename(f, paste(f, ".h5", sep = ""))
58 } 63 }
59 64
60 so <- sleuth_prep(s2c, full_model = ~ condition, num_cores = 1) 65 so <- sleuth_prep(s2c, full_model = ~condition, num_cores = 1)
61 so <- sleuth_fit(so) 66 so <- sleuth_fit(so)
67 } else {
68 ## Simple experiment design
69 ###########################
62 70
63 } else { 71 conditions <- c()
64 ## Simple experiment design 72 for (x in seq_along(args$factorLevel)) {
65 ########################### 73 temp <- append(conditions, rep(args$factorLevel[[x]]))
74 conditions <- temp
75 }
66 76
67 conditions <- c() 77 sample_names <-
68 for (x in seq_along(args$factorLevel)) { 78 gsub(".fastq.+", "", basename(args$factorLevel_counts))
69 temp <- append(conditions, rep(args$factorLevel[[x]]))
70 conditions <- temp
71 }
72 79
73 sample_names <- 80 design <-
74 gsub(".fastq.+", "", basename(args$factorLevel_counts)) 81 data.frame(list(
75 82 sample = sample_names,
76 design <- 83 condition = conditions,
77 data.frame(list( 84 path = args$factorLevel_counts
78 sample = sample_names, 85 ))
79 condition = conditions, 86 so <- sleuth_prep(design,
80 path = args$factorLevel_counts 87 cores = args$cores,
81 )) 88 normalize = args$normalize
82 so <- sleuth_prep(design, 89 )
83 cores = args$cores,
84 normalize = args$normalize)
85 } 90 }
86 91
87 so <- sleuth_fit( 92 so <- sleuth_fit(
88 so, 93 so,
89 ~ condition, 94 ~condition,
90 "full", 95 "full",
91 n_bins = args$nbins, 96 n_bins = args$nbins,
92 lwr = args$lwr, 97 lwr = args$lwr,
93 upr = args$upr 98 upr = args$upr
94 ) 99 )
95 100
96 so <- sleuth_fit( 101 so <- sleuth_fit(
97 so, 102 so,
98 ~ 1, 103 ~1,
99 "reduced", 104 "reduced",
100 n_bins = args$nbins, 105 n_bins = args$nbins,
101 lwr = args$lwr, 106 lwr = args$lwr,
102 upr = args$upr 107 upr = args$upr
103 ) 108 )
104 109
105 so <- sleuth_lrt(so, "reduced", "full") 110 so <- sleuth_lrt(so, "reduced", "full")
106 111
107 sleuth_table <- 112 sleuth_table <-
108 sleuth_results(so, "reduced:full", "lrt", show_all = FALSE) 113 sleuth_results(so, "reduced:full", "lrt", show_all = FALSE)
109 114
110 write.table( 115 write.table(
111 sleuth_table, 116 sleuth_table,
112 file = "sleuth_table.tab", 117 file = "sleuth_table.tab",
113 quote = FALSE, 118 quote = FALSE,
114 sep = "\t", 119 sep = "\t",
115 col.names = TRUE, 120 col.names = TRUE,
116 row.names = FALSE 121 row.names = FALSE
117 ) 122 )
118 123
119 124
120 outputFile <- file.path(getwd(), "pca_plot.pdf") 125 outputFile <- file.path(getwd(), "pca_plot.pdf")
121 pdf(file = outputFile, 126 pdf(
127 file = outputFile,
122 height = 6, 128 height = 6,
123 width = 9) 129 width = 9
130 )
124 plot_pca(so, color_by = "condition") 131 plot_pca(so, color_by = "condition")
125 dev.off() 132 dev.off()
126 133
127 outputFile <- file.path(getwd(), "group_density.pdf") 134 outputFile <- file.path(getwd(), "group_density.pdf")
128 pdf(file = outputFile, 135 pdf(
136 file = outputFile,
129 height = 6, 137 height = 6,
130 width = 9) 138 width = 9
139 )
131 plot_group_density( 140 plot_group_density(
132 so, 141 so,
133 use_filtered = TRUE, 142 use_filtered = TRUE,
134 units = "est_counts", 143 units = "est_counts",
135 trans = "log", 144 trans = "log",
136 grouping = setdiff(colnames(so$sample_to_covariates), 145 grouping = setdiff(
137 "sample"), 146 colnames(so$sample_to_covariates),
138 offset = 1 147 "sample"
148 ),
149 offset = 1
139 ) 150 )
140 dev.off() 151 dev.off()