Mercurial > repos > iuc > sleuth
comparison sleuth.R @ 3:2c3d294dbe42 draft default tip
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() |