Mercurial > repos > sblanck > smagexp
diff MetaRNASeq.R @ 6:3ce32282f6a4 draft
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 | 58052f8bc987 |
| children |
line wrap: on
line diff
--- a/MetaRNASeq.R Thu Mar 22 11:16:24 2018 -0400 +++ b/MetaRNASeq.R Tue Jun 26 08:54:45 2018 -0400 @@ -10,6 +10,9 @@ ##### Read options option_list=list( make_option("--input",type="character",default="NULL",help="list of rdata objects containing eset objects"), + make_option("--inputName",type="character",default=NULL,help="filenames of the Rddata objects"), + make_option("--nbreplicates",type="character",default=NULL,help="number of replicate per study"), + make_option("--fdr",type="character",default=NULL,help="Adjusted p-value threshold to be declared differentially expressed"), make_option("--result",type="character",default=NULL,help="text file containing result of the meta-analysis"), make_option("--htmloutput",type="character",default=NULL,help="Output html report"), make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"), @@ -31,26 +34,25 @@ suppressPackageStartupMessages(require(VennDiagram)) suppressPackageStartupMessages(require(GEOquery)) -listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) ) - -listfiles=vector() -listfilenames=vector() +listfiles <- trimws( unlist( strsplit(trimws(opt$input), ";") ) ) +listfilenames <- trimws( unlist( strsplit(trimws(opt$inputName), ";") ) ) +nbreplicates <- as.numeric(trimws( unlist( strsplit(trimws(opt$nbreplicates), ";") ) )) -for (i in 1:length(listInput)) -{ - inputFileInfo <- unlist( strsplit( listInput[i], ';' ) ) - listfiles=c(listfiles,inputFileInfo[1]) - listfilenames=c(listfilenames,inputFileInfo[2]) -} +#for (i in 1:length(listInput)) +#{ +# inputFileInfo <- unlist( strsplit( listInput[i], ';' ) ) +# listfiles=c(listfiles,inputFileInfo[1]) +# listfilenames=c(listfilenames,inputFileInfo[2]) +# nbreplicates[i]=as.numeric(inputFileInfo[3]) +#} + outputfile <- opt$result result.html = opt$htmloutput html.files.path=opt$htmloutputpath result.template=opt$htmltemplate -alpha=0.05 - -#print(comparison) +alpha=as.numeric(opt$fdr) listData=lapply(listfiles,read.table) orderData=lapply(listData, function(x) x[order(x[1]), ]) @@ -67,95 +69,13 @@ colnames(DE)=listfilenames FC=as.data.frame(FC) colnames(FC)=listfilenames -# the comparison must only have two values and the conds must -# be a vector from those values, at least one of each. - -#if (length(comparison) != 2) { -# stop("Comparison type must be a tuple: ", cargs[length(cargs) - 8]) -#} - -sink("/dev/null") -dir.create(html.files.path, recursive=TRUE) -#library(DESeq) -#library(HTSFilter) - -#DE=list() -#FC=list() -#i=1 - -# Open the html output file -#file.conn = file(diag.html, open="w") - -#writeLines( c("<html><body>"), file.conn) - -# Perform deseq analysis on each study -#for(i in 1:length(listfiles)) -#{ -# f=listfiles[i] -# fname=listfilenames[i] -# study_name=unlist(strsplit(fname,"[.]"))[1] -# print(paste0("study.name ",study_name)) -# d <- read.table(f, sep=" ", header=TRUE, row.names=1) -# conds<-sapply(strsplit(colnames(d),"[.]"),FUN=function(x) x[1]) -# if (length(unique(conds)) != 2) { -# warning(as.data.frame(strsplit(colnames(d),"[.]"))) -# stop("You can only have two columns types: ", paste(conds,collapse=" ")) -# } -# if (!identical(sort(comparison), sort(unique(conds)))) { -# stop("Column types must use the two names from Comparison type, and vice versa. Must have at least one of each in the Column types.\nColumn types: ", cargs[2], "\n", "Comparison type: ", cargs[3]) -# } -# if (length(d) != length(conds)) { -# stop("Number of total sample columns in counts file must correspond to the columns types field. E.g. if column types is 'kidney,kidney,liver,liver' then number of sample columns in counts file must be 4 as well.") -# } -# -# cds <- newCountDataSet(d, conds) -# cds <- estimateSizeFactors(cds) -# -# cdsBlind <- estimateDispersions( cds, method="blind" ) -# -# if (length(conds) != 2) { -# cds <- estimateDispersions( cds ) -# norep = FALSE -# } -# -# if (length(conds) == 2) { -# cds <- estimateDispersions( cds, method=method, sharingMode=mod, fitType="parametric" ) -# norep = TRUE -# } -# -# filter<-HTSFilter(cds, plot=FALSE) -# cds.filter<-filter$filteredData -# on.index<-which(filter$on==1) -# -# res<-as.data.frame(matrix(NA,nrow=nrow(cds),ncol=ncol(cds))) -# nbT <- nbinomTest(cds.filter, comparison[1], comparison[2]) -# colnames(res)<-colnames(nbT) -# res[on.index,]<-nbT -# #write.table(res[order(res$padj), ], file=outputfile, quote=FALSE, row.names=FALSE, sep="\t") -# -# -# temp.pval.plot = file.path( html.files.path, paste("PvalHist",i,".png",sep="")) -# png( temp.pval.plot, width=500, height=500 ) -# hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="") -# dev.off() -# -# writeLines( c("<h2>P-value histogram for ",study_name,"</h2>"), file.conn) -# writeLines( c("<img src='PvalHist",i,".png'><br/><br/>"), file.conn) -# -# #on enregistre la p-value -# rawpval[[study_name]]<-res$pval -# DE[[study_name]]<-ifelse(res$padj<=alpha,1,0) -# FC[[study_name]]<-res$log2FoldChange -# -# i=i+1 -#} - # combinations library(metaRNASeq) +library(UpSetR) fishcomb<-fishercomb(rawpval, BHth=alpha) warning(length(rawpval)) -invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha) +invnormcomb<-invnorm(rawpval, nrep=nbreplicates, BHth=alpha) #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0) #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0) @@ -180,17 +100,30 @@ IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de) IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de) -#conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],Fishercomb=DEresults[["DE.fishercomb"]],Invnormcomb=DEresults[["DE.invnorm"]],sign=commonsgnFC) conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC) #write DE outputfile write.table(conflits, outputfile,sep="\t",,row.names=FALSE) library(VennDiagram) -DE_num=apply(DEresults, 2, FUN=function(x) which(x==1)) -venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6) +DE_num=apply(keepDE[,1:(length(listfiles)+2)], 2, FUN=function(x) which(x==1)) +#DE_num=apply(DEresults, 2, FUN=function(x) which(x==1)) +dir.create(html.files.path, showWarnings = TRUE, recursive = FALSE) temp.venn.plot = file.path( html.files.path, paste("venn.png")) -png(temp.venn.plot,width=500,height=500) -grid.draw(venn.plot) -dev.off() +if (length(listfiles)<=2) { + title="VENN DIAGRAM" + width=500 + venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6) + png(temp.venn.plot,width=width,height=500) + grid.draw(venn.plot) + dev.off() +} else { + title="UPSETR DIAGRAM" + width=1000 + png(temp.venn.plot,width=width,height=500) + upset(fromList(as.list(DE_num)), order.by = "freq") + dev.off() + +} + library(jsonlite) matrixConflits=as.matrix(conflits) @@ -198,11 +131,6 @@ summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE) summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE) - -#vennsplit=strsplit(result.venn,split="/")[[1]] -#venn=paste0("./",vennsplit[length(vennsplit)]) - - vennFilename="venn.png" vennFile=file.path(html.files.path,vennFilename) htmlfile=readChar(result.template, file.info(result.template)$size) @@ -210,39 +138,7 @@ htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE) htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE) htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE) +htmlfile=gsub(x=htmlfile,pattern = "###WIDTH###",replacement = as.character(width), fixed = TRUE) +htmlfile=gsub(x=htmlfile,pattern = "###TITLE###",replacement = title, 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) - - -#writeLines( c("<h2>Venn Plot</h2>"), file.conn) -#writeLines( c("<img src='venn.png'><br/><br/>"), file.conn) -#writeLines( c("</body></html>"), file.conn) -#close(file.conn) -#print("passe6") -#sink(NULL)