Mercurial > repos > sblanck > smagexp

comparison MetaMA.R @ 0:1024245abc70 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 5974f806f344dbcc384b931492d7f023bfbbe03b
author sblanck
date Thu, 22 Feb 2018 08:38:22 -0500
parents
children f18413a94742
comparison
equal deleted inserted replaced
-1:000000000000 0:1024245abc70
1 #!/usr/bin/env Rscript
2 # setup R error handling to go to stderr
3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
4
5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
7
8 library("optparse")
9
10 ##### Read options
11 option_list=list(
12 make_option("--input",type="character",default="NULL",help="List of rdata objects containing eset objects"),
13 make_option("--species",type="character",default="NULL",help="Species for annotations"),
14 make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
15 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
16 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
17 );
18
19 opt_parser = OptionParser(option_list=option_list);
20 opt = parse_args(opt_parser);
21
22 if(is.null(opt$input)){
23 print_help(opt_parser)
24 stop("input required.", call.=FALSE)
25 }
26 if(is.null(opt$species)){
27 print_help(opt_parser)
28 stop("input required.", call.=FALSE)
29 }
30
31 #loading libraries
32
33 suppressPackageStartupMessages(require(metaMA))
34 suppressPackageStartupMessages(require(affy))
35 suppressPackageStartupMessages(require(annaffy))
36 suppressPackageStartupMessages(require(VennDiagram))
37 suppressPackageStartupMessages(require(GEOquery))
38
39 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
40 species=opt$species
41
42 rdataList=list()
43 condition1List=list()
44 condition2List=list()
45 for (input in listInput)
46 {
47 load(input)
48
49 rdataList=c(rdataList,(eset))
50 condition1List=c(condition1List,saveConditions[1])
51 condition2List=c(condition2List,saveConditions[2])
52
53 }
54
55 result.html<-opt$htmloutput
56 result.path<-opt$htmloutputpath
57 result.template<-opt$htmltemplate
58
59 showVenn<-function(res,file)
60 {
61 venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),
62 filename = NULL, col = "black",
63 fill = c(1:(length(res)-2)),
64 margin=0.05, alpha = 0.6)
65 jpeg(file)
66 grid.draw(venn.plot)
67 dev.off()
68 }
69
70 if(!species %in% installed.packages()[,"Package"]) {
71 source("https://bioconductor.org/biocLite.R")
72 biocLite(species)
73 }
74
75 #library("org.Hs.eg.db")
76 x <- org.Hs.egUNIGENE
77 mapped_genes <- mappedkeys(x)
78 link <- as.list(x[mapped_genes])
79
80 probe2unigene<-function(expset){
81 #construction of the map probe->unigene
82 probes=rownames(exprs(expset))
83 gene_id=fData(expset)[probes,"ENTREZ_GENE_ID"]
84 unigene=link[gene_id]
85 names(unigene)<-probes
86 probe_unigene=unigene
87 }
88
89 unigene2probe<-function(map)
90 {
91 suppressWarnings(x <- cbind(unlist(map), names(map)))
92 unigene_probe=split(x[,2], x[,1])
93 }
94
95 convert2metaMA<-function(listStudies,mergemeth=mean)
96 {
97 if (!(class(listStudies) %in% c("list"))) {
98 stop("listStudies must be a list")
99 }
100 conv_unigene=lapply(listStudies,
101 FUN=function(x) unigene2probe(probe2unigene(x)))
102
103 id=lapply(conv_unigene,names)
104 inter=Reduce(intersect,id)
105 if(length(inter)<=0){stop("no common genes")}
106 print(paste(length(inter),"genes in common"))
107 esets=lapply(1:length(listStudies),FUN=function(i){
108 l=lapply(conv_unigene[[i]][inter],
109 FUN=function(x) exprs(listStudies[[i]])[x,,drop=TRUE])
110 esetsgr=t(sapply(l,FUN=function(ll) if(is.null(dim(ll))){ll}
111 else{apply(ll,2,mergemeth)}))
112 esetsgr
113 })
114 return(list(esets=esets,conv.unigene=conv_unigene))
115 }
116
117 normalization<-function(data){
118 ex <- exprs(data)
119 qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
120 LogC <- (qx[5] > 100) ||
121 (qx[6]-qx[1] > 50 && qx[2] > 0) ||
122 (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
123 if (LogC) { ex[which(ex <= 0)] <- NaN
124 return (log2(ex)) } else {
125 return (ex)
126 }
127 }
128
129
130 filterCondition<-function(gset,condition1, condition2){
131 selected=c(which((tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition1)),
132 which(tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition2))
133
134 return(gset[,selected])
135 }
136
137 rdatalist <- lapply(rdataList, FUN=function(datalist) normalization(datalist))
138
139 classes=list()
140 filteredRdataList=list()
141 for (i in 1:length(rdatalist))
142 {
143 currentData=rdataList[[i]]
144 currentCondition1=condition1List[[i]]
145 currentCondition2=condition2List[[i]]
146 #currentData=filterCondition(currentData,currentCondition1,currentCondition2)
147 currentClasses=as.numeric(tolower(as.character(pData(currentData)["source_name_ch1"][,1]))==currentCondition1)
148 filteredRdataList=c(filteredRdataList,currentData)
149 classes=c(classes,list(currentClasses))
150 #write(file="~/galaxy-modal/classes.txt",classes)
151 }
152
153 #rdataList=filteredRdataList
154 conv=convert2metaMA(rdataList)
155 esets=conv$esets
156 conv_unigene=conv$conv.unigene
157
158 #write(file="~/galaxy-modal/esets.txt",length(esets))
159 #write(file="~/galaxy-modal/classes.txt",length(classes))
160 res=pvalcombination(esets=esets,classes=classes,moderated="limma")
161 resIDDIRR=IDDIRR(res$Meta,res$AllIndStudies)
162 length(res$Meta)
163 Hs.Meta=rownames(esets[[1]])[res$Meta]
164 origId.Meta=lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta])))
165
166 gpllist <- lapply(rdataList, FUN=function(ann) annotation(ann))
167 platflist <- lapply(gpllist, FUN=function(gpl) getGEO(gpl, AnnotGPL=TRUE))
168 ncbifdlist <- lapply(platflist, FUN=function(data) data.frame(attr(dataTable(data), "table")))
169 ncbifdresult=lapply(1:length(origId.Meta), FUN=function(i) ncbifdlist[[i]][which(ncbifdlist[[i]]$ID %in% origId.Meta[[i]]),])
170 ncbidfannot=do.call(rbind,ncbifdresult)
171 ncbidfannot <- subset(ncbidfannot, select=c("Platform_SPOTID","ID","Gene.symbol","Gene.title","Gene.ID","Chromosome.annotation","GO.Function.ID"))
172
173 library(jsonlite)
174 matrixncbidfannot=as.matrix(ncbidfannot)
175 datajson=toJSON(matrixncbidfannot,pretty = TRUE)
176 summaryjson=toJSON(as.matrix(t(resIDDIRR)),pretty = TRUE)
177
178
179 #vennsplit=strsplit(result.venn,split="/")[[1]]
180 #venn=paste0("./",vennsplit[length(vennsplit)])
181
182
183 vennFilename="venn.png"
184 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)
192 flog.threshold(ERROR)
193
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")
195 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
196
197 showVenn<-function(liste,file)
198 {
199 venn.plot<-venn.diagram(x = liste,
200 filename = vennFilename, col = "black",
201 fill = 1:length(liste)+1,
202 margin=0.05, alpha = 0.6,imagetype = "png")
203 # png(file);
204 # grid.draw(venn.plot);
205 # dev.off();
206
207 }
208
209 l=list()
210 for(i in 1:length(esets))
211 {
212 l[[paste("study",i,sep="")]]<-res[[i]]
213 }
214 l[["Meta"]]=res[[length(res)-1]]
215 showVenn(l,vennFile)
216 file.copy(vennFilename,result.path)
217 #l=list()
218 #for(i in 1:length(esets))
219 #{
220 # l[[paste("study",i,sep="")]]<-res[[i]]
221 #}
222 #l[["Meta"]]=res[[length(res)-1]]
223 #showVenn(res,result.venn)
224 #writeLines(c("<h2>Venn diagram</h2>"),file.conn)
225 #writeLines(c("<img src='venn.png'><br/><br/>"),file.conn)
226 #writeLines(c("</body></html>"),file.conn)
227 #close(file.conn)