# HG changeset patch
# User sblanck
# Date 1529583936 14400
# Node ID 2c0df656094a6ad44ff0114652c8afe1c04ceb86
# Parent 8191cc63589cecae943b3f0f08cc83df493d1102
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit 071df0c8b84f98fb643b2b9eaa5e238b94ffa254
diff -r 8191cc63589c -r 2c0df656094a MetaRNASeq.R
--- a/MetaRNASeq.R Thu Jun 21 08:05:35 2018 -0400
+++ b/MetaRNASeq.R Thu Jun 21 08:25:36 2018 -0400
@@ -54,8 +54,6 @@
alpha=as.numeric(opt$fdr)
-#print(comparison)
-
listData=lapply(listfiles,read.table)
orderData=lapply(listData, function(x) x[order(x[1]), ])
rawpval=lapply(orderData,function(x) x[6])
@@ -71,99 +69,12 @@
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(""), 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("P-value histogram for ",study_name,"
"), file.conn)
-# writeLines( c("![](PvalHist",i,".png)
"), 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))
-
-#nrep=c(length(listFiles))
-
-
invnormcomb<-invnorm(rawpval, nrep=nbreplicates, BHth=alpha)
#DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0)
#DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0)
@@ -189,18 +100,17 @@
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(keepDE[,1:(length(listfiles)+2)], 2, FUN=function(x) which(x==1))
#DE_num=apply(DEresults, 2, FUN=function(x) which(x==1))
-if (length(listfiles<=2)) {
+temp.venn.plot = file.path( html.files.path, paste("venn.png"))
+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)
- temp.venn.plot = file.path( html.files.path, paste("venn.png"))
png(temp.venn.plot,width=width,height=500)
grid.draw(venn.plot)
dev.off()
@@ -220,11 +130,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)
@@ -232,41 +137,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)
+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("Venn Plot
"), file.conn)
-#writeLines( c("
"), file.conn)
-#writeLines( c(""), file.conn)
-#close(file.conn)
-#print("passe6")
-#sink(NULL)