Mercurial > repos > cristian > rbgoa
comparison GO_MWU.R @ 1:f7287f82602f draft
"planemo upload commit 486235d6560c9e95bd42152ad19bf7c3941cdc1b"
author | cristian |
---|---|
date | Tue, 19 Apr 2022 08:28:43 +0000 |
parents | 91261b42c07e |
children | 5acf9dfdfa27 |
comparison
equal
deleted
inserted
replaced
0:91261b42c07e | 1:f7287f82602f |
---|---|
118 source(paste(base_dir, fname, sep = "/")) | 118 source(paste(base_dir, fname, sep = "/")) |
119 } | 119 } |
120 | 120 |
121 source_local("gomwu.functions.R") | 121 source_local("gomwu.functions.R") |
122 | 122 |
123 | |
124 nn <- strsplit(opt$input, "[/.]") | |
125 if (length(nn[[1]]) == 3) { | |
126 dir <- nn[[1]][1] | |
127 name <- nn[[1]][2] | |
128 ext <- nn[[1]][3] | |
129 } else if (length(nn[[1]]) == 2) { | |
130 dir <- "." | |
131 name <- nn[[1]][1] | |
132 ext <- nn[[1]][2] | |
133 } | |
123 # It might take a few minutes for MF and BP. Do not rerun it if you just want | 134 # It might take a few minutes for MF and BP. Do not rerun it if you just want |
124 # to replot the data with different cutoffs, go straight to gomwuPlot. If you | 135 # to replot the data with different cutoffs, go straight to gomwuPlot. If you |
125 # change any of the numeric values below, delete the files that were generated | 136 # change any of the numeric values below, delete the files that were generated |
126 # in previos runs first. | 137 # in previos runs first. |
127 | 138 |
143 | 154 |
144 | 155 |
145 # ----------- Plotting results | 156 # ----------- Plotting results |
146 | 157 |
147 # change this to a pdf output | 158 # change this to a pdf output |
148 pdf() | 159 pdf(paste0(dir,"/","Rplots.pdf")) |
149 | 160 |
150 results <- gomwuPlot(opt$input, opt$goAnnotations, opt$goDivision, | 161 results <- gomwuPlot(opt$input, opt$goAnnotations, opt$goDivision, |
151 # genes with the measure value exceeding this will be counted as "good genes". | 162 # genes with the measure value exceeding this will be counted as "good genes". |
152 # This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing | 163 # This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing |
153 # Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). | 164 # Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). |
166 ) | 177 ) |
167 # manually rescale the plot so the tree matches the text | 178 # manually rescale the plot so the tree matches the text |
168 # if there are too many categories displayed, try make it more stringent with level1 = 0.05,level2=0.01,level3=0.001. | 179 # if there are too many categories displayed, try make it more stringent with level1 = 0.05,level2=0.01,level3=0.001. |
169 | 180 |
170 # text representation of results, with actual adjusted p-values | 181 # text representation of results, with actual adjusted p-values |
171 write.table(results[[1]], "results.tsv", sep = "\t") | 182 write.table(results[[1]], paste0(dir, "/", "results.tsv"), sep = "\t") |
172 | 183 |
173 | 184 |
174 # ------- extracting representative GOs | 185 # ------- extracting representative GOs |
175 | 186 |
176 # this module chooses GO terms that best represent *independent* groups of significant GO terms | 187 # this module chooses GO terms that best represent *independent* groups of significant GO terms |
210 if (bestrr$pval[best] <= pcut) { | 221 if (bestrr$pval[best] <= pcut) { |
211 annots <- c(annots, sub("\\d+\\/\\d+ ", "", row.names(bestrr)[best])) | 222 annots <- c(annots, sub("\\d+\\/\\d+ ", "", row.names(bestrr)[best])) |
212 } | 223 } |
213 } | 224 } |
214 | 225 |
215 mwus <- read.table(paste("MWU", opt$goDivision, opt$input, sep = "_"), header = T) | 226 mwus <- read.table(paste0(dir,"/",paste("MWU", opt$goDivision, name, sep = "_"), ".", ext), header = T) |
216 best_GOs <- mwus[mwus$name %in% annots, ] | 227 best_GOs <- mwus[mwus$name %in% annots, ] |
217 write.table(best_GOs, "best_go.tsv", sep = "\t") | 228 write.table(best_GOs, paste0(dir, "/","best_go.tsv"), sep = "\t") |
218 dev.off() | 229 dev.off() |