Mercurial > repos > sblanck > smagexp

comparison MetaRNASeq.R @ 6:3ce32282f6a4 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 228bd64b01f184d5d8ebc9f65fe0add2d45ed4fe
author sblanck
date Tue, 26 Jun 2018 08:54:45 -0400
parents 58052f8bc987
children
comparison
equal deleted inserted replaced
5:ca46ad51fe5a 6:3ce32282f6a4
8 library("optparse") 8 library("optparse")
9 9
10 ##### Read options 10 ##### Read options
11 option_list=list( 11 option_list=list(
12 make_option("--input",type="character",default="NULL",help="list of rdata objects containing eset objects"), 12 make_option("--input",type="character",default="NULL",help="list of rdata objects containing eset objects"),
13 make_option("--inputName",type="character",default=NULL,help="filenames of the Rddata objects"),
14 make_option("--nbreplicates",type="character",default=NULL,help="number of replicate per study"),
15 make_option("--fdr",type="character",default=NULL,help="Adjusted p-value threshold to be declared differentially expressed"),
13 make_option("--result",type="character",default=NULL,help="text file containing result of the meta-analysis"), 16 make_option("--result",type="character",default=NULL,help="text file containing result of the meta-analysis"),
14 make_option("--htmloutput",type="character",default=NULL,help="Output html report"), 17 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"), 18 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
16 make_option("--htmltemplate",type="character",default=NULL,help="html template)") 19 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
17 ); 20 );
29 suppressPackageStartupMessages(require(affy)) 32 suppressPackageStartupMessages(require(affy))
30 suppressPackageStartupMessages(require(annaffy)) 33 suppressPackageStartupMessages(require(annaffy))
31 suppressPackageStartupMessages(require(VennDiagram)) 34 suppressPackageStartupMessages(require(VennDiagram))
32 suppressPackageStartupMessages(require(GEOquery)) 35 suppressPackageStartupMessages(require(GEOquery))
33 36
34 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) ) 37 listfiles <- trimws( unlist( strsplit(trimws(opt$input), ";") ) )
38 listfilenames <- trimws( unlist( strsplit(trimws(opt$inputName), ";") ) )
39 nbreplicates <- as.numeric(trimws( unlist( strsplit(trimws(opt$nbreplicates), ";") ) ))
35 40
36 listfiles=vector() 41 #for (i in 1:length(listInput))
37 listfilenames=vector() 42 #{
43 # inputFileInfo <- unlist( strsplit( listInput[i], ';' ) )
44 # listfiles=c(listfiles,inputFileInfo[1])
45 # listfilenames=c(listfilenames,inputFileInfo[2])
46 # nbreplicates[i]=as.numeric(inputFileInfo[3])
47 #}
38 48
39 for (i in 1:length(listInput))
40 {
41 inputFileInfo <- unlist( strsplit( listInput[i], ';' ) )
42 listfiles=c(listfiles,inputFileInfo[1])
43 listfilenames=c(listfilenames,inputFileInfo[2])
44 }
45 49
46 outputfile <- opt$result 50 outputfile <- opt$result
47 result.html = opt$htmloutput 51 result.html = opt$htmloutput
48 html.files.path=opt$htmloutputpath 52 html.files.path=opt$htmloutputpath
49 result.template=opt$htmltemplate 53 result.template=opt$htmltemplate
50 54
51 alpha=0.05 55 alpha=as.numeric(opt$fdr)
52
53 #print(comparison)
54 56
55 listData=lapply(listfiles,read.table) 57 listData=lapply(listfiles,read.table)
56 orderData=lapply(listData, function(x) x[order(x[1]), ]) 58 orderData=lapply(listData, function(x) x[order(x[1]), ])
57 rawpval=lapply(orderData,function(x) x[6]) 59 rawpval=lapply(orderData,function(x) x[6])
58 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x))) 60 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x)))
65 67
66 DE=as.data.frame(DE) 68 DE=as.data.frame(DE)
67 colnames(DE)=listfilenames 69 colnames(DE)=listfilenames
68 FC=as.data.frame(FC) 70 FC=as.data.frame(FC)
69 colnames(FC)=listfilenames 71 colnames(FC)=listfilenames
70 # the comparison must only have two values and the conds must
71 # be a vector from those values, at least one of each.
72
73 #if (length(comparison) != 2) {
74 # stop("Comparison type must be a tuple: ", cargs[length(cargs) - 8])
75 #}
76
77 sink("/dev/null")
78 dir.create(html.files.path, recursive=TRUE)
79 #library(DESeq)
80 #library(HTSFilter)
81
82 #DE=list()
83 #FC=list()
84 #i=1
85
86 # Open the html output file
87 #file.conn = file(diag.html, open="w")
88
89 #writeLines( c("<html><body>"), file.conn)
90
91 # Perform deseq analysis on each study
92 #for(i in 1:length(listfiles))
93 #{
94 # f=listfiles[i]
95 # fname=listfilenames[i]
96 # study_name=unlist(strsplit(fname,"[.]"))[1]
97 # print(paste0("study.name ",study_name))
98 # d <- read.table(f, sep=" ", header=TRUE, row.names=1)
99 # conds<-sapply(strsplit(colnames(d),"[.]"),FUN=function(x) x[1])
100 # if (length(unique(conds)) != 2) {
101 # warning(as.data.frame(strsplit(colnames(d),"[.]")))
102 # stop("You can only have two columns types: ", paste(conds,collapse=" "))
103 # }
104 # if (!identical(sort(comparison), sort(unique(conds)))) {
105 # 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])
106 # }
107 # if (length(d) != length(conds)) {
108 # 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.")
109 # }
110 #
111 # cds <- newCountDataSet(d, conds)
112 # cds <- estimateSizeFactors(cds)
113 #
114 # cdsBlind <- estimateDispersions( cds, method="blind" )
115 #
116 # if (length(conds) != 2) {
117 # cds <- estimateDispersions( cds )
118 # norep = FALSE
119 # }
120 #
121 # if (length(conds) == 2) {
122 # cds <- estimateDispersions( cds, method=method, sharingMode=mod, fitType="parametric" )
123 # norep = TRUE
124 # }
125 #
126 # filter<-HTSFilter(cds, plot=FALSE)
127 # cds.filter<-filter$filteredData
128 # on.index<-which(filter$on==1)
129 #
130 # res<-as.data.frame(matrix(NA,nrow=nrow(cds),ncol=ncol(cds)))
131 # nbT <- nbinomTest(cds.filter, comparison[1], comparison[2])
132 # colnames(res)<-colnames(nbT)
133 # res[on.index,]<-nbT
134 # #write.table(res[order(res$padj), ], file=outputfile, quote=FALSE, row.names=FALSE, sep="\t")
135 #
136 #
137 # temp.pval.plot = file.path( html.files.path, paste("PvalHist",i,".png",sep=""))
138 # png( temp.pval.plot, width=500, height=500 )
139 # hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="")
140 # dev.off()
141 #
142 # writeLines( c("<h2>P-value histogram for ",study_name,"</h2>"), file.conn)
143 # writeLines( c("<img src='PvalHist",i,".png'><br/><br/>"), file.conn)
144 #
145 # #on enregistre la p-value
146 # rawpval[[study_name]]<-res$pval
147 # DE[[study_name]]<-ifelse(res$padj<=alpha,1,0)
148 # FC[[study_name]]<-res$log2FoldChange
149 #
150 # i=i+1
151 #}
152
153 72
154 # combinations 73 # combinations
155 library(metaRNASeq) 74 library(metaRNASeq)
75 library(UpSetR)
156 fishcomb<-fishercomb(rawpval, BHth=alpha) 76 fishcomb<-fishercomb(rawpval, BHth=alpha)
157 warning(length(rawpval)) 77 warning(length(rawpval))
158 invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha) 78 invnormcomb<-invnorm(rawpval, nrep=nbreplicates, BHth=alpha)
159 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0) 79 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0)
160 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0) 80 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0)
161 81
162 signsFC<-mapply(FC,FUN=function(x) sign(x)) 82 signsFC<-mapply(FC,FUN=function(x) sign(x))
163 sumsigns<-apply(signsFC,1,sum) 83 sumsigns<-apply(signsFC,1,sum)
178 } 98 }
179 99
180 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de) 100 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de)
181 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de) 101 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de)
182 102
183 #conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],Fishercomb=DEresults[["DE.fishercomb"]],Invnormcomb=DEresults[["DE.invnorm"]],sign=commonsgnFC)
184 conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC) 103 conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC)
185 #write DE outputfile 104 #write DE outputfile
186 write.table(conflits, outputfile,sep="\t",,row.names=FALSE) 105 write.table(conflits, outputfile,sep="\t",,row.names=FALSE)
187 library(VennDiagram) 106 library(VennDiagram)
188 DE_num=apply(DEresults, 2, FUN=function(x) which(x==1)) 107 DE_num=apply(keepDE[,1:(length(listfiles)+2)], 2, FUN=function(x) which(x==1))
189 venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6) 108 #DE_num=apply(DEresults, 2, FUN=function(x) which(x==1))
109 dir.create(html.files.path, showWarnings = TRUE, recursive = FALSE)
190 temp.venn.plot = file.path( html.files.path, paste("venn.png")) 110 temp.venn.plot = file.path( html.files.path, paste("venn.png"))
191 png(temp.venn.plot,width=500,height=500) 111 if (length(listfiles)<=2) {
192 grid.draw(venn.plot) 112 title="VENN DIAGRAM"
193 dev.off() 113 width=500
114 venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6)
115 png(temp.venn.plot,width=width,height=500)
116 grid.draw(venn.plot)
117 dev.off()
118 } else {
119 title="UPSETR DIAGRAM"
120 width=1000
121 png(temp.venn.plot,width=width,height=500)
122 upset(fromList(as.list(DE_num)), order.by = "freq")
123 dev.off()
124
125 }
126
194 127
195 library(jsonlite) 128 library(jsonlite)
196 matrixConflits=as.matrix(conflits) 129 matrixConflits=as.matrix(conflits)
197 datajson=toJSON(matrixConflits,pretty = TRUE) 130 datajson=toJSON(matrixConflits,pretty = TRUE)
198 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE) 131 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE)
199 summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE) 132 summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE)
200
201
202 #vennsplit=strsplit(result.venn,split="/")[[1]]
203 #venn=paste0("./",vennsplit[length(vennsplit)])
204
205 133
206 vennFilename="venn.png" 134 vennFilename="venn.png"
207 vennFile=file.path(html.files.path,vennFilename) 135 vennFile=file.path(html.files.path,vennFilename)
208 htmlfile=readChar(result.template, file.info(result.template)$size) 136 htmlfile=readChar(result.template, file.info(result.template)$size)
209 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE) 137 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
210 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE) 138 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE)
211 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE) 139 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE)
212 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE) 140 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
141 htmlfile=gsub(x=htmlfile,pattern = "###WIDTH###",replacement = as.character(width), fixed = TRUE)
142 htmlfile=gsub(x=htmlfile,pattern = "###TITLE###",replacement = title, fixed = TRUE)
213 write(htmlfile,result.html) 143 write(htmlfile,result.html)
214 144
215 #library(VennDiagram)
216 #flog.threshold(ERROR)
217 #
218 ##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")
219 #dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
220 #
221 #showVenn<-function(liste,file)
222 #{
223 # venn.plot<-venn.diagram(x = liste,
224 # filename = vennFilename, col = "black",
225 # fill = 1:length(liste)+1,
226 # margin=0.05, alpha = 0.6,imagetype = "png")
227 ## png(file);
228 ## grid.draw(venn.plot);
229 ## dev.off();
230 #
231 #}
232 #
233 #l=list()
234 #for(i in 1:length(esets))
235 #{
236 # l[[paste("study",i,sep="")]]<-res[[i]]
237 #}
238 #l[["Meta"]]=res[[length(res)-1]]
239 #showVenn(l,vennFile)
240 #file.copy(vennFilename,result.path)
241
242
243 #writeLines( c("<h2>Venn Plot</h2>"), file.conn)
244 #writeLines( c("<img src='venn.png'><br/><br/>"), file.conn)
245 #writeLines( c("</body></html>"), file.conn)
246 #close(file.conn)
247 #print("passe6")
248 #sink(NULL)