Mercurial > repos > sblanck > smagexp
annotate MetaRNASeq.R @ 24:cb29ad7d75a5 draft
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit a0cb2703e4dbbc1fadeedc5b61acb6e1324433f7
author | sblanck |
---|---|
date | Sat, 24 Feb 2018 10:07:51 -0500 |
parents | ef7d98f9eb51 |
children | 49ce6fbe11de |
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 |
e5a94bc69bd6
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160
sblanck
parents:
2
diff
changeset
|
46 outputfile <- opt$result |
e5a94bc69bd6
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160
sblanck
parents:
2
diff
changeset
|
47 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
|
48 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
|
49 result.template=opt$htmltemplate |
2 | 50 |
51 alpha=0.05 | |
52 | |
53 #print(comparison) | |
54 | |
55 listData=lapply(listfiles,read.table) | |
56 orderData=lapply(listData, function(x) x[order(x[1]), ]) | |
57 rawpval=lapply(orderData,function(x) x[6]) | |
58 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x))) | |
59 | |
60 DE=list() | |
61 DE=lapply(orderData, function(x) ifelse(x[7]<=0.05,1,0)) | |
62 | |
63 FC=list() | |
64 FC=lapply(orderData, function(x) x[3]) | |
65 | |
66 DE=as.data.frame(DE) | |
67 colnames(DE)=listfilenames | |
68 FC=as.data.frame(FC) | |
69 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 | |
154 # combinations | |
155 library(metaRNASeq) | |
156 fishcomb<-fishercomb(rawpval, BHth=alpha) | |
157 warning(length(rawpval)) | |
158 invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha) | |
159 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0) | |
160 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0) | |
161 | |
162 signsFC<-mapply(FC,FUN=function(x) sign(x)) | |
163 sumsigns<-apply(signsFC,1,sum) | |
164 commonsgnFC<-ifelse(abs(sumsigns)==dim(signsFC)[2],sign(sumsigns),0) | |
165 | |
166 DEresults <- data.frame(DE=DE,"DE.fishercomb"=ifelse(fishcomb$adjpval<=alpha,1,0),"DE.invnorm"=ifelse(invnormcomb$adjpval<=alpha,1,0)) | |
167 | |
168 unionDE <- unique(c(fishcomb$DEindices,invnormcomb$DEindices)) | |
169 FC.selecDE <- data.frame(DEresults[unionDE,],FC[unionDE,],signFC=commonsgnFC[unionDE]) | |
170 keepDE <- FC.selecDE[which(abs(FC.selecDE$signFC)==1),] | |
171 | |
172 fishcomb_de <- rownames(keepDE)[which(keepDE[,"DE.fishercomb"]==1)] | |
173 invnorm_de <- rownames(keepDE)[which(keepDE[,"DE.invnorm"]==1)] | |
174 indstudy_de = list() | |
175 for (i in 1:length(listfiles)) { | |
176 currentIndstudy_de = rownames(keepDE)[which(keepDE[,i]==1)] | |
177 indstudy_de[[listfilenames[i]]]=currentIndstudy_de | |
178 } | |
179 | |
180 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de) | |
181 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de) | |
182 | |
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) | |
185 #write DE outputfile | |
186 write.table(conflits, outputfile,sep="\t",,row.names=FALSE) | |
187 library(VennDiagram) | |
188 DE_num=apply(DEresults, 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) | |
190 temp.venn.plot = file.path( html.files.path, paste("venn.png")) | |
191 png(temp.venn.plot,width=500,height=500) | |
192 grid.draw(venn.plot) | |
193 dev.off() | |
194 | |
195 library(jsonlite) | |
196 matrixConflits=as.matrix(conflits) | |
197 datajson=toJSON(matrixConflits,pretty = TRUE) | |
198 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE) | |
199 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 | |
206 vennFilename="venn.png" | |
207 vennFile=file.path(html.files.path,vennFilename) | |
208 htmlfile=readChar(result.template, file.info(result.template)$size) | |
209 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE) | |
210 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE) | |
211 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE) | |
212 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE) | |
213 write(htmlfile,result.html) | |
214 | |
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) |