Mercurial > repos > sblanck > smagexp

comparison MetaRNASeq.R @ 0:1024245abc70 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 5974f806f344dbcc384b931492d7f023bfbbe03b
author sblanck
date Thu, 22 Feb 2018 08:38:22 -0500
parents
children 58052f8bc987
comparison
equal deleted inserted replaced
-1:000000000000 0:1024245abc70
1 #!/usr/bin/env Rscript
2 # setup R error handling to go to stderr
3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
4
5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
7
8 library("optparse")
9
10 ##### Read options
11 option_list=list(
12 make_option("--input",type="character",default="NULL",help="list of rdata objects containing eset objects"),
13 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"),
15 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
16 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
17 );
18
19 opt_parser = OptionParser(option_list=option_list);
20 opt = parse_args(opt_parser);
21
22 if(is.null(opt$input)){
23 print_help(opt_parser)
24 stop("input required.", call.=FALSE)
25 }
26
27 #loading libraries
28
29 suppressPackageStartupMessages(require(affy))
30 suppressPackageStartupMessages(require(annaffy))
31 suppressPackageStartupMessages(require(VennDiagram))
32 suppressPackageStartupMessages(require(GEOquery))
33
34 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
35
36 listfiles=vector()
37 listfilenames=vector()
38
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
46 cargs <- commandArgs()
47 cargs <- cargs[(which(cargs == "--args")+1):length(cargs)]
48 nbargs=length(cargs)
49 listfiles=vector()
50 listfilenames=vector()
51 for (i in seq(1,nbargs-6,2)) {
52 listfiles=c(listfiles,cargs[[i]])
53 listfilenames=c(listfilenames,cargs[[i+1]])
54 }
55 #mod<-cargs[[length(cargs) - 6]]
56 outputfile <- opt$result
57 result.html = opt$htmloutput
58 html.files.path=opt$htmloutputpath
59 result.template=opt$htmltemplate
60
61 alpha=0.05
62
63 #print(comparison)
64
65 listData=lapply(listfiles,read.table)
66 orderData=lapply(listData, function(x) x[order(x[1]), ])
67 rawpval=lapply(orderData,function(x) x[6])
68 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x)))
69
70 DE=list()
71 DE=lapply(orderData, function(x) ifelse(x[7]<=0.05,1,0))
72
73 FC=list()
74 FC=lapply(orderData, function(x) x[3])
75
76 DE=as.data.frame(DE)
77 colnames(DE)=listfilenames
78 FC=as.data.frame(FC)
79 colnames(FC)=listfilenames
80 # the comparison must only have two values and the conds must
81 # be a vector from those values, at least one of each.
82
83 #if (length(comparison) != 2) {
84 # stop("Comparison type must be a tuple: ", cargs[length(cargs) - 8])
85 #}
86
87 sink("/dev/null")
88 dir.create(html.files.path, recursive=TRUE)
89 #library(DESeq)
90 #library(HTSFilter)
91
92 #DE=list()
93 #FC=list()
94 #i=1
95
96 # Open the html output file
97 #file.conn = file(diag.html, open="w")
98
99 #writeLines( c("<html><body>"), file.conn)
100
101 # Perform deseq analysis on each study
102 #for(i in 1:length(listfiles))
103 #{
104 # f=listfiles[i]
105 # fname=listfilenames[i]
106 # study_name=unlist(strsplit(fname,"[.]"))[1]
107 # print(paste0("study.name ",study_name))
108 # d <- read.table(f, sep=" ", header=TRUE, row.names=1)
109 # conds<-sapply(strsplit(colnames(d),"[.]"),FUN=function(x) x[1])
110 # if (length(unique(conds)) != 2) {
111 # warning(as.data.frame(strsplit(colnames(d),"[.]")))
112 # stop("You can only have two columns types: ", paste(conds,collapse=" "))
113 # }
114 # if (!identical(sort(comparison), sort(unique(conds)))) {
115 # 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])
116 # }
117 # if (length(d) != length(conds)) {
118 # 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.")
119 # }
120 #
121 # cds <- newCountDataSet(d, conds)
122 # cds <- estimateSizeFactors(cds)
123 #
124 # cdsBlind <- estimateDispersions( cds, method="blind" )
125 #
126 # if (length(conds) != 2) {
127 # cds <- estimateDispersions( cds )
128 # norep = FALSE
129 # }
130 #
131 # if (length(conds) == 2) {
132 # cds <- estimateDispersions( cds, method=method, sharingMode=mod, fitType="parametric" )
133 # norep = TRUE
134 # }
135 #
136 # filter<-HTSFilter(cds, plot=FALSE)
137 # cds.filter<-filter$filteredData
138 # on.index<-which(filter$on==1)
139 #
140 # res<-as.data.frame(matrix(NA,nrow=nrow(cds),ncol=ncol(cds)))
141 # nbT <- nbinomTest(cds.filter, comparison[1], comparison[2])
142 # colnames(res)<-colnames(nbT)
143 # res[on.index,]<-nbT
144 # #write.table(res[order(res$padj), ], file=outputfile, quote=FALSE, row.names=FALSE, sep="\t")
145 #
146 #
147 # temp.pval.plot = file.path( html.files.path, paste("PvalHist",i,".png",sep=""))
148 # png( temp.pval.plot, width=500, height=500 )
149 # hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="")
150 # dev.off()
151 #
152 # writeLines( c("<h2>P-value histogram for ",study_name,"</h2>"), file.conn)
153 # writeLines( c("<img src='PvalHist",i,".png'><br/><br/>"), file.conn)
154 #
155 # #on enregistre la p-value
156 # rawpval[[study_name]]<-res$pval
157 # DE[[study_name]]<-ifelse(res$padj<=alpha,1,0)
158 # FC[[study_name]]<-res$log2FoldChange
159 #
160 # i=i+1
161 #}
162
163
164 # combinations
165 library(metaRNASeq)
166 fishcomb<-fishercomb(rawpval, BHth=alpha)
167 warning(length(rawpval))
168 invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha)
169 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0)
170 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0)
171
172 signsFC<-mapply(FC,FUN=function(x) sign(x))
173 sumsigns<-apply(signsFC,1,sum)
174 commonsgnFC<-ifelse(abs(sumsigns)==dim(signsFC)[2],sign(sumsigns),0)
175
176 DEresults <- data.frame(DE=DE,"DE.fishercomb"=ifelse(fishcomb$adjpval<=alpha,1,0),"DE.invnorm"=ifelse(invnormcomb$adjpval<=alpha,1,0))
177
178 unionDE <- unique(c(fishcomb$DEindices,invnormcomb$DEindices))
179 FC.selecDE <- data.frame(DEresults[unionDE,],FC[unionDE,],signFC=commonsgnFC[unionDE])
180 keepDE <- FC.selecDE[which(abs(FC.selecDE$signFC)==1),]
181
182 fishcomb_de <- rownames(keepDE)[which(keepDE[,"DE.fishercomb"]==1)]
183 invnorm_de <- rownames(keepDE)[which(keepDE[,"DE.invnorm"]==1)]
184 indstudy_de = list()
185 for (i in 1:length(listfiles)) {
186 currentIndstudy_de = rownames(keepDE)[which(keepDE[,i]==1)]
187 indstudy_de[[listfilenames[i]]]=currentIndstudy_de
188 }
189
190 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de)
191 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de)
192
193 #conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],Fishercomb=DEresults[["DE.fishercomb"]],Invnormcomb=DEresults[["DE.invnorm"]],sign=commonsgnFC)
194 conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC)
195 #write DE outputfile
196 write.table(conflits, outputfile,sep="\t",,row.names=FALSE)
197 library(VennDiagram)
198 DE_num=apply(DEresults, 2, FUN=function(x) which(x==1))
199 venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6)
200 temp.venn.plot = file.path( html.files.path, paste("venn.png"))
201 png(temp.venn.plot,width=500,height=500)
202 grid.draw(venn.plot)
203 dev.off()
204
205 library(jsonlite)
206 matrixConflits=as.matrix(conflits)
207 datajson=toJSON(matrixConflits,pretty = TRUE)
208 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE)
209 summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE)
210
211
212 #vennsplit=strsplit(result.venn,split="/")[[1]]
213 #venn=paste0("./",vennsplit[length(vennsplit)])
214
215
216 vennFilename="venn.png"
217 vennFile=file.path(html.files.path,vennFilename)
218 htmlfile=readChar(result.template, file.info(result.template)$size)
219 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
220 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE)
221 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE)
222 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
223 write(htmlfile,result.html)
224
225 #library(VennDiagram)
226 #flog.threshold(ERROR)
227 #
228 ##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")
229 #dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
230 #
231 #showVenn<-function(liste,file)
232 #{
233 # venn.plot<-venn.diagram(x = liste,
234 # filename = vennFilename, col = "black",
235 # fill = 1:length(liste)+1,
236 # margin=0.05, alpha = 0.6,imagetype = "png")
237 ## png(file);
238 ## grid.draw(venn.plot);
239 ## dev.off();
240 #
241 #}
242 #
243 #l=list()
244 #for(i in 1:length(esets))
245 #{
246 # l[[paste("study",i,sep="")]]<-res[[i]]
247 #}
248 #l[["Meta"]]=res[[length(res)-1]]
249 #showVenn(l,vennFile)
250 #file.copy(vennFilename,result.path)
251
252
253 #writeLines( c("<h2>Venn Plot</h2>"), file.conn)
254 #writeLines( c("<img src='venn.png'><br/><br/>"), file.conn)
255 #writeLines( c("</body></html>"), file.conn)
256 #close(file.conn)
257 #print("passe6")
258 #sink(NULL)