diff MetaMA.R @ 0:1024245abc70 draft

planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit 5974f806f344dbcc384b931492d7f023bfbbe03b

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MetaMA.R	Thu Feb 22 08:38:22 2018 -0500
@@ -0,0 +1,227 @@
+#!/usr/bin/env Rscript
+# setup R error handling to go to stderr
+options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
+
+# we need that to not crash galaxy with an UTF8 error on German LC settings.
+loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+
+library("optparse")
+
+##### Read options
+option_list=list(
+		make_option("--input",type="character",default="NULL",help="List of rdata objects containing eset objects"),
+		make_option("--species",type="character",default="NULL",help="Species for annotations"),
+		make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
+		make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
+		make_option("--htmltemplate",type="character",default=NULL,help="html template)")
+);
+
+opt_parser = OptionParser(option_list=option_list);
+opt = parse_args(opt_parser);
+
+if(is.null(opt$input)){
+	print_help(opt_parser)
+	stop("input required.", call.=FALSE)
+}
+if(is.null(opt$species)){
+	print_help(opt_parser)
+	stop("input required.", call.=FALSE)
+}
+
+#loading libraries
+
+suppressPackageStartupMessages(require(metaMA))
+suppressPackageStartupMessages(require(affy))
+suppressPackageStartupMessages(require(annaffy))
+suppressPackageStartupMessages(require(VennDiagram))
+suppressPackageStartupMessages(require(GEOquery))
+
+listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
+species=opt$species
+
+rdataList=list()
+condition1List=list()
+condition2List=list()
+for (input in listInput)
+{
+	load(input)
+	
+	rdataList=c(rdataList,(eset))
+	condition1List=c(condition1List,saveConditions[1])
+	condition2List=c(condition2List,saveConditions[2])
+	
+}
+
+result.html<-opt$htmloutput
+result.path<-opt$htmloutputpath
+result.template<-opt$htmltemplate
+
+showVenn<-function(res,file)
+{
+	venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),
+			filename = NULL, col = "black", 
+			fill = c(1:(length(res)-2)),
+			margin=0.05, alpha = 0.6)
+	jpeg(file)
+	grid.draw(venn.plot)
+	dev.off()
+}
+
+if(!species %in% installed.packages()[,"Package"]) {
+	source("https://bioconductor.org/biocLite.R")
+	biocLite(species)
+}
+
+#library("org.Hs.eg.db")
+x <- org.Hs.egUNIGENE
+mapped_genes <- mappedkeys(x)
+link <- as.list(x[mapped_genes])
+
+probe2unigene<-function(expset){
+	#construction of the map probe->unigene
+	probes=rownames(exprs(expset))
+	gene_id=fData(expset)[probes,"ENTREZ_GENE_ID"]
+	unigene=link[gene_id]
+	names(unigene)<-probes
+	probe_unigene=unigene
+}
+
+unigene2probe<-function(map)
+{
+	suppressWarnings(x <- cbind(unlist(map), names(map)))
+	unigene_probe=split(x[,2], x[,1])
+}
+
+convert2metaMA<-function(listStudies,mergemeth=mean)
+{
+	if (!(class(listStudies) %in% c("list"))) {
+		stop("listStudies must be a list")
+	}
+	conv_unigene=lapply(listStudies,
+			FUN=function(x) unigene2probe(probe2unigene(x)))
+	
+	id=lapply(conv_unigene,names)
+	inter=Reduce(intersect,id)
+	if(length(inter)<=0){stop("no common genes")}
+	print(paste(length(inter),"genes in common"))
+	esets=lapply(1:length(listStudies),FUN=function(i){
+				l=lapply(conv_unigene[[i]][inter],
+						FUN=function(x) exprs(listStudies[[i]])[x,,drop=TRUE])
+				esetsgr=t(sapply(l,FUN=function(ll) if(is.null(dim(ll))){ll}
+									else{apply(ll,2,mergemeth)}))
+				esetsgr
+			})
+	return(list(esets=esets,conv.unigene=conv_unigene))
+}
+
+normalization<-function(data){
+	ex <- exprs(data)
+	qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
+	LogC <- (qx[5] > 100) ||
+			(qx[6]-qx[1] > 50 && qx[2] > 0) ||
+			(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
+	if (LogC) { ex[which(ex <= 0)] <- NaN
+		return (log2(ex)) } else {
+		return (ex)
+	}
+}
+
+
+filterCondition<-function(gset,condition1, condition2){
+	selected=c(which((tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition1)), 
+			which(tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition2))
+	
+	return(gset[,selected])
+	}
+
+rdatalist <- lapply(rdataList, FUN=function(datalist) normalization(datalist))
+
+classes=list()
+filteredRdataList=list()
+for (i in 1:length(rdatalist))
+{
+	currentData=rdataList[[i]]
+	currentCondition1=condition1List[[i]]
+	currentCondition2=condition2List[[i]]
+	#currentData=filterCondition(currentData,currentCondition1,currentCondition2)
+	currentClasses=as.numeric(tolower(as.character(pData(currentData)["source_name_ch1"][,1]))==currentCondition1)
+	filteredRdataList=c(filteredRdataList,currentData)
+	classes=c(classes,list(currentClasses))
+	#write(file="~/galaxy-modal/classes.txt",classes)
+}
+
+#rdataList=filteredRdataList
+conv=convert2metaMA(rdataList)
+esets=conv$esets
+conv_unigene=conv$conv.unigene
+
+#write(file="~/galaxy-modal/esets.txt",length(esets))
+#write(file="~/galaxy-modal/classes.txt",length(classes))
+res=pvalcombination(esets=esets,classes=classes,moderated="limma")
+resIDDIRR=IDDIRR(res$Meta,res$AllIndStudies)
+length(res$Meta)
+Hs.Meta=rownames(esets[[1]])[res$Meta]
+origId.Meta=lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta])))
+
+gpllist <- lapply(rdataList, FUN=function(ann) annotation(ann))
+platflist <- lapply(gpllist, FUN=function(gpl) getGEO(gpl, AnnotGPL=TRUE))
+ncbifdlist <- lapply(platflist, FUN=function(data) data.frame(attr(dataTable(data), "table")))
+ncbifdresult=lapply(1:length(origId.Meta), FUN=function(i) ncbifdlist[[i]][which(ncbifdlist[[i]]$ID %in% origId.Meta[[i]]),])
+ncbidfannot=do.call(rbind,ncbifdresult)
+ncbidfannot <- subset(ncbidfannot, select=c("Platform_SPOTID","ID","Gene.symbol","Gene.title","Gene.ID","Chromosome.annotation","GO.Function.ID"))
+
+library(jsonlite)
+matrixncbidfannot=as.matrix(ncbidfannot)
+datajson=toJSON(matrixncbidfannot,pretty = TRUE)
+summaryjson=toJSON(as.matrix(t(resIDDIRR)),pretty = TRUE)
+
+
+#vennsplit=strsplit(result.venn,split="/")[[1]]
+#venn=paste0("./",vennsplit[length(vennsplit)])
+
+
+vennFilename="venn.png"
+vennFile=file.path(result.path,vennFilename)
+htmlfile=readChar(result.template, file.info(result.template)$size)
+htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
+htmlfile=gsub(x=htmlfile,pattern = "###SUMMARYJSON###",replacement = summaryjson, fixed = TRUE)
+htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
+write(htmlfile,result.html)
+
+library(VennDiagram)
+flog.threshold(ERROR)
+
+#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")
+dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
+
+showVenn<-function(liste,file)
+{
+	venn.plot<-venn.diagram(x = liste,
+			filename = vennFilename, col = "black",
+			fill = 1:length(liste)+1,
+			margin=0.05, alpha = 0.6,imagetype = "png")
+#	png(file);
+#	grid.draw(venn.plot);
+#	dev.off();
+	
+}
+
+l=list()
+for(i in 1:length(esets))
+{
+	l[[paste("study",i,sep="")]]<-res[[i]]
+}
+l[["Meta"]]=res[[length(res)-1]]
+showVenn(l,vennFile)
+file.copy(vennFilename,result.path)
+#l=list()
+#for(i in 1:length(esets))
+#{
+#	l[[paste("study",i,sep="")]]<-res[[i]]
+#}
+#l[["Meta"]]=res[[length(res)-1]]
+#showVenn(res,result.venn)
+#writeLines(c("<h2>Venn diagram</h2>"),file.conn)
+#writeLines(c("<img src='venn.png'><br/><br/>"),file.conn)
+#writeLines(c("</body></html>"),file.conn)
+#close(file.conn)
\ No newline at end of file