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