Mercurial > repos > vandelj > giant_gsea_format

comparison src/VolcanoPlotsScript.R @ 0:3022feec50fe draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/juliechevalier/GIANT/tree/master commit cb276a594444c8f32e9819fefde3a21f121d35df"
author vandelj
date Fri, 26 Jun 2020 09:36:46 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:3022feec50fe
1 # R script to plot volcanos through Galaxy based GIANT tool
2 # written by Jimmy Vandel
3 #
4 #
5 initial.options <- commandArgs(trailingOnly = FALSE)
6 file.arg.name <- "--file="
7 script.name <- sub(file.arg.name, "", initial.options[grep(file.arg.name, initial.options)])
8 script.basename <- dirname(script.name)
9 source(file.path(script.basename, "utils.R"))
10 source(file.path(script.basename, "getopt.R"))
11
12 #addComment("Welcome R!")
13
14 # setup R error handling to go to stderr
15 options( show.error.messages=F, error = function () { cat(geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
16
17 # we need that to not crash galaxy with an UTF8 error on German LC settings.
18 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
19 loc <- Sys.setlocale("LC_NUMERIC", "C")
20
21 #get starting time
22 start.time <- Sys.time()
23
24 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
25 args <- commandArgs()
26
27 # get options, using the spec as defined by the enclosed list.
28 # we read the options from the default: commandArgs(TRUE).
29 spec <- matrix(c(
30 "statisticsFile", "i", 1, "character",
31 "volcanoName" , "n", 1, "character",
32 "pvalColumnName" , "p", 1, "character",
33 "fdrColumnName" , "m", 1, "character",
34 "fcColumnName" , "c", 1, "character",
35 "fcKind","d", 1, "character",
36 "fdrThreshold","s", 1, "double",
37 "fcThreshold","e", 1, "double",
38 "organismID","x",1,"character",
39 "rowNameType","y",1,"character",
40 "log", "l", 1, "character",
41 "outputFile" , "o", 1, "character",
42 "format", "f", 1, "character",
43 "quiet", "q", 0, "logical"),
44 byrow=TRUE, ncol=4)
45 opt <- getopt(spec)
46
47 # enforce the following required arguments
48 if (is.null(opt$log)) {
49 addComment("[ERROR]'log file' is required\n")
50 q( "no", 1, F )
51 }
52 addComment("[INFO]Start of R script",T,opt$log,display=FALSE)
53 if (is.null(opt$statisticsFile)) {
54 addComment("[ERROR]'statisticsFile' is required",T,opt$log)
55 q( "no", 1, F )
56 }
57 if (length(opt$pvalColumnName)==0 || length(opt$fdrColumnName)==0 || length(opt$fcColumnName)==0) {
58 addComment("[ERROR]no selected columns",T,opt$log)
59 q( "no", 1, F )
60 }
61 if (length(opt$pvalColumnName)!=length(opt$fcColumnName) || length(opt$pvalColumnName)!=length(opt$fdrColumnName)) {
62 addComment("[ERROR]different number of selected columns between p.val, adj-p.val and FC ",T,opt$log)
63 q( "no", 1, F )
64 }
65 if (is.null(opt$fcKind)) {
66 addComment("[ERROR]'fcKind' is required",T,opt$log)
67 q( "no", 1, F )
68 }
69 if (is.null(opt$fdrThreshold)) {
70 addComment("[ERROR]'FDR threshold' is required",T,opt$log)
71 q( "no", 1, F )
72 }
73 if (is.null(opt$fcThreshold)) {
74 addComment("[ERROR]'FC threshold' is required",T,opt$log)
75 q( "no", 1, F )
76 }
77 if (is.null(opt$outputFile)) {
78 addComment("[ERROR]'output file' is required",T,opt$log)
79 q( "no", 1, F )
80 }
81 if (is.null(opt$format)) {
82 addComment("[ERROR]'output format' is required",T,opt$log)
83 q( "no", 1, F )
84 }
85
86 #demande si le script sera bavard
87 verbose <- if (is.null(opt$quiet)) {
88 TRUE
89 }else{
90 FALSE
91 }
92
93 #paramètres internes
94 addComment("[INFO]Parameters checked test mode !",T,opt$log,display=FALSE)
95
96 addComment(c("[INFO]Working directory: ",getwd()),TRUE,opt$log,display=FALSE)
97 addComment(c("[INFO]Command line: ",args),TRUE,opt$log,display=FALSE)
98
99 #directory for plots
100 dir.create(file.path(getwd(), "plotDir"))
101 dir.create(file.path(getwd(), "plotLyDir"))
102
103 #charge des packages silencieusement
104 suppressPackageStartupMessages({
105 library("methods")
106 library("biomaRt")
107 library("ggplot2")
108 library("plotly")
109 library("stringr")
110 })
111
112 #define some usefull variable
113 nbVolcanosToPlot=length(opt$pvalColumnName)
114
115 #load input file
116 statDataMatrix=read.csv(file=file.path(getwd(), opt$statisticsFile),header=F,sep="\t",colClasses="character")
117 #remove first colum to convert it as rownames
118 rownames(statDataMatrix)=statDataMatrix[,1]
119 statDataMatrix=statDataMatrix[,-1]
120
121 #identify lines without adjusted p-value info (should contain the same content as rownames) and replace them with NA values
122 FDRinfo=rep(TRUE,nbVolcanosToPlot)
123 for(iVolcano in 1:nbVolcanosToPlot){
124 #input parameter should be None when adjusted p-val are not available
125 if(opt$fdrColumnName[iVolcano]=="None"){
126 #content of the corresponding column should also be the same as rownames
127 if(!all(statDataMatrix[,(iVolcano-1)*3+2]==rownames(statDataMatrix))){
128 addComment(c("[ERROR]It seems that input stat matrix contains adjusted p-values for volcano",iVolcano,"whereas input parameter indicates that not."),T,opt$log)
129 q( "no", 1, F )
130 }
131 FDRinfo[iVolcano]=FALSE
132 statDataMatrix[,(iVolcano-1)*3+2]=NA
133 }
134 }
135
136 if(is.data.frame(statDataMatrix)){
137 statDataMatrix=data.matrix(statDataMatrix)
138 }else{
139 statDataMatrix=data.matrix(as.numeric(statDataMatrix))
140 }
141
142 #check if available column number match with volcano requested number
143 if(ncol(statDataMatrix)!=3*nbVolcanosToPlot){
144 addComment("[ERROR]Input file column number is different from requested volcano number",T,opt$log)
145 q( "no", 1, F )
146 }
147
148 #build global dataFrame with data and fill with p.val and log2(FC) and FDR
149 dataFrame=data.frame(row.names = rownames(statDataMatrix))
150 #start with p-value
151 dataFrame$p.value=statDataMatrix[,seq(1,nbVolcanosToPlot*3,3),drop=FALSE]
152 #compute FDR if needed or just get available info
153 dataFrame$adj_p.value=dataFrame$p.value
154 for(iVolcano in 1:nbVolcanosToPlot){
155 #adjusted p-value are already computed
156 if(FDRinfo[iVolcano]){
157 dataFrame$adj_p.value[,iVolcano]=statDataMatrix[,(iVolcano-1)*3+2,drop=FALSE]
158 }else{
159 #adjusted p-value should be computed based on p-val using FDR
160 dataFrame$adj_p.value[,iVolcano]=p.adjust(dataFrame$p.value[,iVolcano,drop=FALSE],"fdr")
161 addComment(c("[INFO]Adjusted p-values are not available in input for volcano",iVolcano,", FDR approach will be used on available raw p-values"),T,opt$log)
162 }
163 }
164 if(opt$fcKind=="FC"){
165 #we should transform as Log2FC
166 dataFrame$coefficients=log2(statDataMatrix[,seq(3,nbVolcanosToPlot*3,3),drop=FALSE])
167 addComment(c("[INFO]FC are converted in log2(FC) for plotting"),T,opt$log)
168 }else{
169 dataFrame$coefficients=statDataMatrix[,seq(3,nbVolcanosToPlot*3,3),drop=FALSE]
170 }
171
172 addComment(c("[INFO]Input data available for",nbVolcanosToPlot,"volcano(s) with",nrow(statDataMatrix),"rows"),T,opt$log)
173
174
175 #plot VOLCANOs
176 volcanoPerPage=1
177 logFCthreshold=log2(opt$fcThreshold)
178 iToPlot=1
179 plotVector=list()
180 volcanoNameList=c()
181 for (iVolcano in 1:nbVolcanosToPlot){
182
183 if(nchar(opt$volcanoName[iVolcano])>0){
184 curentVolcanoName=opt$volcanoName[iVolcano]
185 }else{
186 curentVolcanoName=paste(iVolcano,opt$pvalColumnName[iVolcano],sep="_")
187 }
188
189 #keep only rows without NA for p-val, adjusted p-val and coeff
190 pValToPlot=dataFrame$p.value[,iVolcano]
191 fdrToPlot=dataFrame$adj_p.value[,iVolcano]
192 coeffToPlot=dataFrame$coefficients[,iVolcano]
193
194 rowToRemove=unique(c(which(is.na(pValToPlot)),which(is.na(fdrToPlot)),which(is.na(coeffToPlot))))
195 if(length(rowToRemove)>0){
196 pValToPlot=pValToPlot[-rowToRemove]
197 fdrToPlot=fdrToPlot[-rowToRemove]
198 coeffToPlot=coeffToPlot[-rowToRemove]
199 }
200 addComment(c("[INFO]For",curentVolcanoName,"volcano,",length(rowToRemove),"rows are discarded due to NA values,",length(pValToPlot),"remaining rows."),T,opt$log)
201
202 #save volcano name
203 volcanoNameList=c(volcanoNameList,curentVolcanoName)
204
205 #remove characters possibly troubling
206 volcanoFileName=iVolcano
207
208 #define the log10(p-val) threshold corresponding to FDR threshold fixed by user
209 probeWithLowFDR=-log10(pValToPlot[which(fdrToPlot<=opt$fdrThreshold)])
210 pvalThresholdFDR=NULL
211 if(length(probeWithLowFDR)>0)pvalThresholdFDR=min(probeWithLowFDR)
212
213 #get significant points over FC and FDR thresholds
214 significativePoints=intersect(which(abs(coeffToPlot)>=logFCthreshold),which(fdrToPlot<=opt$fdrThreshold))
215
216 #to reduce size of html plot, we keep 20000 points maximum sampled amongst genes with pval>=33%(pval) and abs(log2(FC))<=66%(abs(log2(FC)))
217 htmlPointsToRemove=intersect(which(abs(coeffToPlot)<=quantile(abs(coeffToPlot),c(0.66))),which(pValToPlot>=quantile(abs(pValToPlot),c(0.33))))
218 if(length(htmlPointsToRemove)>20000){
219 htmlPointsToRemove=setdiff(htmlPointsToRemove,sample(htmlPointsToRemove,20000))
220 }else{
221 htmlPointsToRemove=c()
222 }
223
224 xMinLimPlot=min(coeffToPlot)-0.2
225 xMaxLimPlot=max(coeffToPlot)+0.2
226 yMaxLimPlot= max(-log10(pValToPlot))+0.2
227
228 if(length(significativePoints)>0){
229 dataSignifToPlot=data.frame(pval=-log10(pValToPlot[significativePoints]),FC=coeffToPlot[significativePoints],description=paste(names(coeffToPlot[significativePoints]),"\n","FC: " , round(2^coeffToPlot[significativePoints],2) , " | Adjusted p-val: ",prettyNum(fdrToPlot[significativePoints],digits=4), sep=""))
230 #to test if remains any normal points to draw
231 if(length(significativePoints)<length(pValToPlot)){
232 dataToPlot=data.frame(pval=-log10(pValToPlot[-significativePoints]),FC=coeffToPlot[-significativePoints],description=paste("FC: " , round(2^coeffToPlot[-significativePoints],2) , " | Adjusted p-val: ",prettyNum(fdrToPlot[-significativePoints],digits=4), sep=""))
233 }else{
234 dataToPlot=data.frame(pval=0,FC=0,description="null")
235 }
236 }else{
237 dataToPlot=data.frame(pval=-log10(pValToPlot),FC=coeffToPlot,description=paste("FC: " , round(2^coeffToPlot,2) , " | Adjusted p-val: ",prettyNum(fdrToPlot,digits=4), sep=""))
238 }
239
240 ##traditional plot
241
242 p <- ggplot(data=dataToPlot, aes(x=FC, y=pval)) + geom_point() +
243 theme_bw() + ggtitle(curentVolcanoName) + ylab(label="-Log10(p-val)") + xlab(label="Log2 Fold Change") +
244 theme(panel.border=element_blank(),plot.title = element_text(hjust = 0.5),legend.position="none")
245 if(logFCthreshold!=0) p <- p + geom_vline(xintercept=-logFCthreshold, color="salmon",linetype="dotted", size=1) + geom_vline(xintercept=logFCthreshold, color="salmon",linetype="dotted", size=1) + geom_text(data.frame(text=c(paste(c("log2(1/FC=",opt$fcThreshold,")"),collapse=""),paste(c("log2(FC=",opt$fcThreshold,")"),collapse="")),x=c(-logFCthreshold,logFCthreshold),y=c(0,0)),mapping=aes(x=x, y=y, label=text), size=4, angle=90, vjust=-0.4, hjust=0, color="salmon")
246 if(!is.null(pvalThresholdFDR)) p <- p + geom_hline(yintercept=pvalThresholdFDR, color="skyblue1",linetype="dotted", size=0.5) + geom_text(data.frame(text=c(paste(c("Adjusted pval limit(",opt$fdrThreshold,")"),collapse="")),x=c(xMinLimPlot),y=c(pvalThresholdFDR)),mapping=aes(x=x, y=y, label=text), size=4, vjust=0, hjust=0, color="skyblue3")
247 if(length(significativePoints)>0)p <- p + geom_point(data=dataSignifToPlot,aes(colour=description))
248
249 ##interactive plot
250
251 if(length(htmlPointsToRemove)>0){
252 pointToRemove=union(htmlPointsToRemove,significativePoints)
253 #to test if it remains any normal points to draw
254 if(length(pointToRemove)<length(pValToPlot)){
255 dataToPlot=data.frame(pval=-log10(pValToPlot[-pointToRemove]),FC=coeffToPlot[-pointToRemove],description=paste("FC: " , round(2^coeffToPlot[-pointToRemove],2) , " | Adjusted p-val: ", prettyNum(fdrToPlot[-pointToRemove],digits=4), sep=""))
256 }else{
257 dataToPlot=data.frame(pval=0,FC=0,description="null")
258 }
259 }
260
261 if((nrow(dataToPlot)+length(significativePoints))>40000)addComment(c("[WARNING]For",curentVolcanoName,"volcano, numerous points to plot(",nrow(dataToPlot)+nrow(dataSignifToPlot),"), resulting volcano could be heavy, using more stringent thresholds could be helpful."),T,opt$log)
262
263 phtml <- plot_ly(data=dataToPlot, x=~FC, y=~pval,type="scatter", mode="markers",showlegend = FALSE, marker = list(color="gray",opacity=0.5), text=~description, hoverinfo="text") %>%
264 layout(title = curentVolcanoName[iVolcano],xaxis=list(title="Log2 Fold Change",showgrid=TRUE, zeroline=FALSE),yaxis=list(title="-Log10(p-val)", showgrid=TRUE, zeroline=FALSE))
265 if(length(significativePoints)>0) phtml=add_markers(phtml,data=dataSignifToPlot, x=~FC, y=~pval, mode="markers" , marker=list( color=log10(abs(dataSignifToPlot$FC)*dataSignifToPlot$pval),colorscale='Rainbow'), text=~description, hoverinfo="text", inherit = FALSE) %>% hide_colorbar()
266 if(logFCthreshold!=0){
267 phtml=add_trace(phtml,x=c(-logFCthreshold,-logFCthreshold), y=c(0,yMaxLimPlot), type="scatter", mode = "lines", line=list(color="coral",dash="dash"), hoverinfo='none', showlegend = FALSE,inherit = FALSE)
268 phtml=add_annotations(phtml,x=-logFCthreshold,y=0,xref = "x",yref = "y",text = paste(c("log2(1/FC=",opt$fcThreshold,")"),collapse=""),xanchor = 'right',showarrow = F,textangle=270,font=list(color="coral"))
269 phtml=add_trace(phtml,x=c(logFCthreshold,logFCthreshold), y=c(0, yMaxLimPlot), type="scatter", mode = "lines", line=list(color="coral",dash="dash"), hoverinfo='none', showlegend = FALSE,inherit = FALSE)
270 phtml=add_annotations(phtml,x=logFCthreshold,y=0,xref = "x",yref = "y",text = paste(c("log2(FC=",opt$fcThreshold,")"),collapse=""),xanchor = 'right',showarrow = F,textangle=270,font=list(color="coral"))
271 }
272 if(!is.null(pvalThresholdFDR)){
273 phtml=add_trace(phtml,x=c(xMinLimPlot,xMaxLimPlot), y=c(pvalThresholdFDR,pvalThresholdFDR), type="scatter", mode = "lines", line=list(color="cornflowerblue",dash="dash"), hoverinfo='none', showlegend = FALSE,inherit = FALSE)
274 phtml=add_annotations(phtml,x=xMinLimPlot,y=pvalThresholdFDR+0.1,xref = "x",yref = "y",text = paste(c("Adjusted pval limit(",opt$fdrThreshold,")"),collapse=""),xanchor = 'left',showarrow = F,font=list(color="cornflowerblue"))
275 }
276 plotVector[[length(plotVector)+1]]=p
277
278 #save plotly files
279 pp <- ggplotly(phtml)
280 htmlwidgets::saveWidget(as_widget(pp), paste(c(file.path(getwd(), "plotLyDir"),"/Volcanos_",volcanoFileName,".html"),collapse=""),selfcontained = F)
281
282
283 if(iVolcano==nbVolcanosToPlot || length(plotVector)==volcanoPerPage){
284 #plot and close the actual plot
285 if(opt$format=="pdf"){
286 pdf(paste(c("./plotDir/Volcanos_",volcanoFileName,".pdf"),collapse=""))}else{
287 png(paste(c("./plotDir/Volcanos_",volcanoFileName,".png"),collapse=""))
288 }
289 multiplot(plotlist=plotVector,cols=1)
290 dev.off()
291 if(iVolcano<nbVolcanosToPlot){
292 #prepare for a new ploting file if necessary
293 plotVector=list()
294 iToPlot=iToPlot+1
295 }
296 }
297 }
298 remove(dataToPlot,dataSignifToPlot)
299 addComment("[INFO]Volcanos drawn",T,opt$log,T,display=FALSE)
300
301
302 #now add anotation infos about genes
303
304 rowItemInfo=NULL
305 if(!is.null(opt$rowNameType) && !is.null(opt$organismID)){
306 ##get gene information from BioMart
307 #if(!require("biomaRt")){
308 # source("https://bioconductor.org/biocLite.R")
309 # biocLite("biomaRt")
310 #}
311
312 ensembl_hs_mart <- useMart(biomart="ensembl", dataset=opt$organismID)
313 ensembl_df <- getBM(attributes=c(opt$rowNameType,"description"),mart=ensembl_hs_mart)
314 rowItemInfo=ensembl_df[which(ensembl_df[,1]!=""),2]
315 rowItemInfo=unlist(lapply(rowItemInfo,function(x)substr(unlist(strsplit(x," \\[Source"))[1],1,30)))
316 names(rowItemInfo)=ensembl_df[which(ensembl_df[,1]!=""),1]
317 }
318
319 #filter out genes with higher p-values for all comparisons
320 genesToKeep=names(which(apply(dataFrame$adj_p.value,1,function(x)length(which(x<=opt$fdrThreshold))>0)))
321 #filter out genes with lower FC for all comparisons
322 genesToKeep=intersect(genesToKeep,names(which(apply(dataFrame$coefficients,1,function(x)length(which(abs(x)>=logFCthreshold))>0))))
323
324 if(length(genesToKeep)>0){
325 dataFrameNew=data.frame(row.names=genesToKeep)
326
327 dataFrameNew$adj_p.value=matrix(dataFrame$adj_p.value[genesToKeep,,drop=FALSE],ncol=ncol(dataFrame$adj_p.value))
328 rownames(dataFrameNew$adj_p.value)=genesToKeep
329 colnames(dataFrameNew$adj_p.value)=colnames(dataFrame$p.value)
330
331 dataFrameNew$p.value=matrix(dataFrame$p.value[genesToKeep,,drop=FALSE],ncol=ncol(dataFrame$p.value))
332 rownames(dataFrameNew$p.value)=genesToKeep
333 colnames(dataFrameNew$p.value)=colnames(dataFrame$adj_p.value)
334
335 dataFrameNew$coefficients=matrix(dataFrame$coefficients[genesToKeep,,drop=FALSE],ncol=ncol(dataFrame$coefficients))
336 rownames(dataFrameNew$coefficients)=genesToKeep
337 colnames(dataFrameNew$coefficients)=colnames(dataFrame$adj_p.value)
338
339 dataFrame=dataFrameNew
340 rm(dataFrameNew)
341 }else{
342 addComment("[WARNING]No significative genes",T,opt$log,display=FALSE)
343 }
344
345 addComment("[INFO]Significant genes filtering done",T,opt$log,T,display=FALSE)
346
347
348 #plot VennDiagramm for genes below threshold between comparisons
349 #t=apply(dataFrame$adj_p.value[,1:4],2,function(x)names(which(x<=opt$threshold)))
350 #get.venn.partitions(t)
351 #vennCounts(dataFrame$adj_p.value[,1:4]<=opt$threshold)
352
353 #make a simple sort genes based only on the first comparison
354 #newOrder=order(dataFrame$adj_p.value[,1])
355 #dataFrame$adj_p.value=dataFrame$adj_p.value[newOrder,]
356
357 #alternative sorting strategy based on the mean gene rank over all comparisons
358 if(length(genesToKeep)>1){
359 currentRank=rep(0,nrow(dataFrame$adj_p.value))
360 for(iVolcano in 1:ncol(dataFrame$adj_p.value)){
361 currentRank=currentRank+rank(dataFrame$adj_p.value[,iVolcano])
362 }
363 currentRank=currentRank/ncol(dataFrame$adj_p.value)
364 newOrder=order(currentRank)
365 rownames(dataFrame)=rownames(dataFrame)[newOrder]
366
367 dataFrame$adj_p.value=matrix(dataFrame$adj_p.value[newOrder,],ncol=ncol(dataFrame$adj_p.value))
368 rownames(dataFrame$adj_p.value)=rownames(dataFrame$p.value)[newOrder]
369 colnames(dataFrame$adj_p.value)=colnames(dataFrame$p.value)
370
371 dataFrame$p.value=matrix(dataFrame$p.value[newOrder,],ncol=ncol(dataFrame$p.value))
372 rownames(dataFrame$p.value)=rownames(dataFrame$adj_p.value)
373 colnames(dataFrame$p.value)=colnames(dataFrame$adj_p.value)
374
375 dataFrame$coefficients=matrix(dataFrame$coefficients[newOrder,],ncol=ncol(dataFrame$coefficients))
376 rownames(dataFrame$coefficients)=rownames(dataFrame$adj_p.value)
377 colnames(dataFrame$coefficients)=colnames(dataFrame$adj_p.value)
378 }
379
380 #formating output matrix depending on genes to keep
381 if(length(genesToKeep)==0){
382 outputData=matrix(0,ncol=ncol(dataFrame$adj_p.value)*4+2,nrow=3)
383 outputData[1,]=c("X","X",rep(volcanoNameList,each=4))
384 outputData[2,]=c("X","X",rep(c("p-val","Adjusted.p-val","FC","log2(FC)"),ncol(dataFrame$adj_p.value)))
385 outputData[,1]=c("Volcano","Gene","noGene")
386 outputData[,2]=c("Comparison","Info","noInfo")
387 }else{
388 if(length(genesToKeep)==1){
389 outputData=matrix(0,ncol=ncol(dataFrame$adj_p.value)*4+2,nrow=3)
390 outputData[1,]=c("X","X",rep(volcanoNameList,each=4))
391 outputData[2,]=c("X","X",rep(c("p-val","Adjusted.p-val","FC","log2(FC)"),ncol(dataFrame$adj_p.value)))
392 outputData[,1]=c("Volcano","Gene",genesToKeep)
393 outputData[,2]=c("Comparison","Info","na")
394 if(!is.null(rowItemInfo))outputData[3,2]=rowItemInfo[genesToKeep]
395 outputData[3,seq(3,ncol(outputData),4)]=prettyNum(dataFrame$p.value,digits=4)
396 outputData[3,seq(4,ncol(outputData),4)]=prettyNum(dataFrame$adj_p.value,digits=4)
397 outputData[3,seq(5,ncol(outputData),4)]=prettyNum(2^dataFrame$coefficients,digits=4)
398 outputData[3,seq(6,ncol(outputData),4)]=prettyNum(dataFrame$coefficients,digits=4)
399 }else{
400 #format matrix to be correctly read by galaxy (move headers in first column and row)
401 outputData=matrix(0,ncol=ncol(dataFrame$adj_p.value)*4+2,nrow=nrow(dataFrame$adj_p.value)+2)
402 outputData[1,]=c("X","X",rep(volcanoNameList,each=4))
403 outputData[2,]=c("X","X",rep(c("p-val","Adjusted.p-val","FC","log2(FC)"),ncol(dataFrame$adj_p.value)))
404 outputData[,1]=c("Volcano","Gene",rownames(dataFrame$adj_p.value))
405 outputData[,2]=c("Comparison","Info",rep("na",nrow(dataFrame$adj_p.value)))
406 if(!is.null(rowItemInfo))outputData[3:nrow(outputData),2]=rowItemInfo[rownames(dataFrame$adj_p.value)]
407 outputData[3:nrow(outputData),seq(3,ncol(outputData),4)]=prettyNum(dataFrame$p.value,digits=4)
408 outputData[3:nrow(outputData),seq(4,ncol(outputData),4)]=prettyNum(dataFrame$adj_p.value,digits=4)
409 outputData[3:nrow(outputData),seq(5,ncol(outputData),4)]=prettyNum(2^dataFrame$coefficients,digits=4)
410 outputData[3:nrow(outputData),seq(6,ncol(outputData),4)]=prettyNum(dataFrame$coefficients,digits=4)
411 }
412 }
413 addComment("[INFO]Formated output",T,opt$log,display=FALSE)
414
415 #write output results
416 write.table(outputData,file=opt$outputFile,quote=FALSE,sep="\t",col.names = F,row.names = F)
417
418
419 end.time <- Sys.time()
420 addComment(c("[INFO]Total execution time for R script:",as.numeric(end.time - start.time,units="mins"),"mins"),T,opt$log,display=FALSE)
421
422 addComment("[INFO]End of R script",T,opt$log,display=FALSE)
423
424 printSessionInfo(opt$log)
425
426 #sessionInfo()