Mercurial > repos > sblanck > smagexp

comparison MetaMA.R @ 6:3ce32282f6a4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit 228bd64b01f184d5d8ebc9f65fe0add2d45ed4fe
author sblanck
date Tue, 26 Jun 2018 08:54:45 -0400
parents f18413a94742
children
comparison
equal deleted inserted replaced
5:ca46ad51fe5a 6:3ce32282f6a4
33 suppressPackageStartupMessages(require(metaMA)) 33 suppressPackageStartupMessages(require(metaMA))
34 suppressPackageStartupMessages(require(affy)) 34 suppressPackageStartupMessages(require(affy))
35 suppressPackageStartupMessages(require(annaffy)) 35 suppressPackageStartupMessages(require(annaffy))
36 suppressPackageStartupMessages(require(VennDiagram)) 36 suppressPackageStartupMessages(require(VennDiagram))
37 suppressPackageStartupMessages(require(GEOquery)) 37 suppressPackageStartupMessages(require(GEOquery))
38 suppressPackageStartupMessages(require(UpSetR))
38 39
39 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) ) 40 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
40 species=opt$species 41 species=opt$species
41 42
42 rdataList=list() 43 rdataList=list()
54 55
55 result.html<-opt$htmloutput 56 result.html<-opt$htmloutput
56 result.path<-opt$htmloutputpath 57 result.path<-opt$htmloutputpath
57 result.template<-opt$htmltemplate 58 result.template<-opt$htmltemplate
58 59
59 showVenn<-function(res,file) 60 #showVenn<-function(res,file)
60 { 61 #{
61 venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)), 62 # venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),
62 filename = NULL, col = "black", 63 # filename = NULL, col = "black",
63 fill = c(1:(length(res)-2)), 64 # fill = c(1:(length(res)-2)),
64 margin=0.05, alpha = 0.6) 65 # margin=0.05, alpha = 0.6)
65 jpeg(file) 66 # jpeg(file)
66 grid.draw(venn.plot) 67 # grid.draw(venn.plot)
67 dev.off() 68 # dev.off()
68 } 69 #}
69 70
70 if(!species %in% installed.packages()[,"Package"]) { 71 if(!species %in% installed.packages()[,"Package"]) {
71 source("https://bioconductor.org/biocLite.R") 72 source("https://bioconductor.org/biocLite.R")
72 biocLite(species) 73 biocLite(species)
73 } 74 }
74 75
75 library("org.Hs.eg.db") 76 library(species,character.only = TRUE)
76 x <- org.Hs.egUNIGENE 77 x <- org.Hs.egUNIGENE
77 mapped_genes <- mappedkeys(x) 78 mapped_genes <- mappedkeys(x)
78 link <- as.list(x[mapped_genes]) 79 link <- as.list(x[mapped_genes])
79 80
80 probe2unigene<-function(expset){ 81 probe2unigene<-function(expset){
157 158
158 #write(file="~/galaxy-modal/esets.txt",length(esets)) 159 #write(file="~/galaxy-modal/esets.txt",length(esets))
159 #write(file="~/galaxy-modal/classes.txt",length(classes)) 160 #write(file="~/galaxy-modal/classes.txt",length(classes))
160 res=pvalcombination(esets=esets,classes=classes,moderated="limma") 161 res=pvalcombination(esets=esets,classes=classes,moderated="limma")
161 resIDDIRR=IDDIRR(res$Meta,res$AllIndStudies) 162 resIDDIRR=IDDIRR(res$Meta,res$AllIndStudies)
162 length(res$Meta) 163 #length(res$Meta)
163 Hs.Meta=rownames(esets[[1]])[res$Meta] 164 Hs.Meta=rownames(esets[[1]])[res$Meta]
164 origId.Meta=lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta]))) 165 origId.Meta=lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta])))
165 166
166 gpllist <- lapply(rdataList, FUN=function(ann) annotation(ann)) 167 gpllist <- lapply(rdataList, FUN=function(ann) annotation(ann))
167 platflist <- lapply(gpllist, FUN=function(gpl) getGEO(gpl, AnnotGPL=TRUE)) 168 platflist <- lapply(gpllist, FUN=function(gpl) getGEO(gpl, AnnotGPL=TRUE))
180 #venn=paste0("./",vennsplit[length(vennsplit)]) 181 #venn=paste0("./",vennsplit[length(vennsplit)])
181 182
182 183
183 vennFilename="venn.png" 184 vennFilename="venn.png"
184 vennFile=file.path(result.path,vennFilename) 185 vennFile=file.path(result.path,vennFilename)
185 htmlfile=readChar(result.template, file.info(result.template)$size)
186 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
187 htmlfile=gsub(x=htmlfile,pattern = "###SUMMARYJSON###",replacement = summaryjson, fixed = TRUE)
188 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
189 write(htmlfile,result.html)
190
191 library(VennDiagram) 186 library(VennDiagram)
192 flog.threshold(ERROR) 187 flog.threshold(ERROR)
193 188
194 #venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),filename = v, col = "black", fill = c(1:(length(res)-2)), margin=0.05, alpha = 0.6,imagetype = "png") 189 #venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),filename = v, col = "black", fill = c(1:(length(res)-2)), margin=0.05, alpha = 0.6,imagetype = "png")
195 dir.create(result.path, showWarnings = TRUE, recursive = FALSE) 190 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
210 for(i in 1:length(esets)) 205 for(i in 1:length(esets))
211 { 206 {
212 l[[paste("study",i,sep="")]]<-res[[i]] 207 l[[paste("study",i,sep="")]]<-res[[i]]
213 } 208 }
214 l[["Meta"]]=res[[length(res)-1]] 209 l[["Meta"]]=res[[length(res)-1]]
215 showVenn(l,vennFile) 210
211 if (length(l)<5) {
212 title="Venn diagram"
213 width=500
214 showVenn(l,vennFile)
215 } else {
216 title="Upsetr diagram"
217 width=1000
218 png(vennFile,width=width,height=500)
219 upset(fromList(as.list(l)), order.by = "freq")
220 dev.off()
221 }
222
216 file.copy(vennFilename,result.path) 223 file.copy(vennFilename,result.path)
224
225 htmlfile=readChar(result.template, file.info(result.template)$size)
226 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
227 htmlfile=gsub(x=htmlfile,pattern = "###SUMMARYJSON###",replacement = summaryjson, fixed = TRUE)
228 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
229 tmlfile=gsub(x=htmlfile,pattern = "###WIDTH###",replacement = as.character(width), fixed = TRUE)
230 htmlfile=gsub(x=htmlfile,pattern = "###TITLE###",replacement = title, fixed = TRUE)
231 write(htmlfile,result.html)
217 #l=list() 232 #l=list()
218 #for(i in 1:length(esets)) 233 #for(i in 1:length(esets))
219 #{ 234 #{
220 # l[[paste("study",i,sep="")]]<-res[[i]] 235 # l[[paste("study",i,sep="")]]<-res[[i]]
221 #} 236 #}