annotate MetaRNASeq.R @ 21:f953fdee364c draft

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