|
2
|
1 cargs <- commandArgs()
|
|
|
2 cargs <- cargs[(which(cargs == "--args")+1):length(cargs)]
|
|
|
3 nbargs=length(cargs)
|
|
|
4 listfiles=vector()
|
|
|
5 listfilenames=vector()
|
|
|
6 for (i in seq(1,nbargs-6,2)) {
|
|
|
7 listfiles=c(listfiles,cargs[[i]])
|
|
|
8 listfilenames=c(listfilenames,cargs[[i+1]])
|
|
|
9 }
|
|
|
10 #mod<-cargs[[length(cargs) - 6]]
|
|
|
11 outputfile <- cargs[[length(cargs) - 5]]
|
|
|
12 result.html = cargs[[length(cargs) - 4]]
|
|
|
13 html.files.path=cargs[[length(cargs) - 3]]
|
|
|
14 result.template=cargs[[length(cargs) - 2]]
|
|
|
15
|
|
|
16 alpha=0.05
|
|
|
17
|
|
|
18 #print(comparison)
|
|
|
19
|
|
|
20 listData=lapply(listfiles,read.table)
|
|
|
21 orderData=lapply(listData, function(x) x[order(x[1]), ])
|
|
|
22 rawpval=lapply(orderData,function(x) x[6])
|
|
|
23 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x)))
|
|
|
24
|
|
|
25 DE=list()
|
|
|
26 DE=lapply(orderData, function(x) ifelse(x[7]<=0.05,1,0))
|
|
|
27
|
|
|
28 FC=list()
|
|
|
29 FC=lapply(orderData, function(x) x[3])
|
|
|
30
|
|
|
31 DE=as.data.frame(DE)
|
|
|
32 colnames(DE)=listfilenames
|
|
|
33 FC=as.data.frame(FC)
|
|
|
34 colnames(FC)=listfilenames
|
|
|
35 # the comparison must only have two values and the conds must
|
|
|
36 # be a vector from those values, at least one of each.
|
|
|
37
|
|
|
38 #if (length(comparison) != 2) {
|
|
|
39 # stop("Comparison type must be a tuple: ", cargs[length(cargs) - 8])
|
|
|
40 #}
|
|
|
41
|
|
|
42 sink("/dev/null")
|
|
|
43 dir.create(html.files.path, recursive=TRUE)
|
|
|
44 #library(DESeq)
|
|
|
45 #library(HTSFilter)
|
|
|
46
|
|
|
47 #DE=list()
|
|
|
48 #FC=list()
|
|
|
49 #i=1
|
|
|
50
|
|
|
51 # Open the html output file
|
|
|
52 #file.conn = file(diag.html, open="w")
|
|
|
53
|
|
|
54 #writeLines( c("<html><body>"), file.conn)
|
|
|
55
|
|
|
56 # Perform deseq analysis on each study
|
|
|
57 #for(i in 1:length(listfiles))
|
|
|
58 #{
|
|
|
59 # f=listfiles[i]
|
|
|
60 # fname=listfilenames[i]
|
|
|
61 # study_name=unlist(strsplit(fname,"[.]"))[1]
|
|
|
62 # print(paste0("study.name ",study_name))
|
|
|
63 # d <- read.table(f, sep=" ", header=TRUE, row.names=1)
|
|
|
64 # conds<-sapply(strsplit(colnames(d),"[.]"),FUN=function(x) x[1])
|
|
|
65 # if (length(unique(conds)) != 2) {
|
|
|
66 # warning(as.data.frame(strsplit(colnames(d),"[.]")))
|
|
|
67 # stop("You can only have two columns types: ", paste(conds,collapse=" "))
|
|
|
68 # }
|
|
|
69 # if (!identical(sort(comparison), sort(unique(conds)))) {
|
|
|
70 # 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])
|
|
|
71 # }
|
|
|
72 # if (length(d) != length(conds)) {
|
|
|
73 # 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.")
|
|
|
74 # }
|
|
|
75 #
|
|
|
76 # cds <- newCountDataSet(d, conds)
|
|
|
77 # cds <- estimateSizeFactors(cds)
|
|
|
78 #
|
|
|
79 # cdsBlind <- estimateDispersions( cds, method="blind" )
|
|
|
80 #
|
|
|
81 # if (length(conds) != 2) {
|
|
|
82 # cds <- estimateDispersions( cds )
|
|
|
83 # norep = FALSE
|
|
|
84 # }
|
|
|
85 #
|
|
|
86 # if (length(conds) == 2) {
|
|
|
87 # cds <- estimateDispersions( cds, method=method, sharingMode=mod, fitType="parametric" )
|
|
|
88 # norep = TRUE
|
|
|
89 # }
|
|
|
90 #
|
|
|
91 # filter<-HTSFilter(cds, plot=FALSE)
|
|
|
92 # cds.filter<-filter$filteredData
|
|
|
93 # on.index<-which(filter$on==1)
|
|
|
94 #
|
|
|
95 # res<-as.data.frame(matrix(NA,nrow=nrow(cds),ncol=ncol(cds)))
|
|
|
96 # nbT <- nbinomTest(cds.filter, comparison[1], comparison[2])
|
|
|
97 # colnames(res)<-colnames(nbT)
|
|
|
98 # res[on.index,]<-nbT
|
|
|
99 # #write.table(res[order(res$padj), ], file=outputfile, quote=FALSE, row.names=FALSE, sep="\t")
|
|
|
100 #
|
|
|
101 #
|
|
|
102 # temp.pval.plot = file.path( html.files.path, paste("PvalHist",i,".png",sep=""))
|
|
|
103 # png( temp.pval.plot, width=500, height=500 )
|
|
|
104 # hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="")
|
|
|
105 # dev.off()
|
|
|
106 #
|
|
|
107 # writeLines( c("<h2>P-value histogram for ",study_name,"</h2>"), file.conn)
|
|
|
108 # writeLines( c("<img src='PvalHist",i,".png'><br/><br/>"), file.conn)
|
|
|
109 #
|
|
|
110 # #on enregistre la p-value
|
|
|
111 # rawpval[[study_name]]<-res$pval
|
|
|
112 # DE[[study_name]]<-ifelse(res$padj<=alpha,1,0)
|
|
|
113 # FC[[study_name]]<-res$log2FoldChange
|
|
|
114 #
|
|
|
115 # i=i+1
|
|
|
116 #}
|
|
|
117
|
|
|
118
|
|
|
119 # combinations
|
|
|
120 library(metaRNASeq)
|
|
|
121 fishcomb<-fishercomb(rawpval, BHth=alpha)
|
|
|
122 warning(length(rawpval))
|
|
|
123 invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha)
|
|
|
124 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0)
|
|
|
125 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0)
|
|
|
126
|
|
|
127 signsFC<-mapply(FC,FUN=function(x) sign(x))
|
|
|
128 sumsigns<-apply(signsFC,1,sum)
|
|
|
129 commonsgnFC<-ifelse(abs(sumsigns)==dim(signsFC)[2],sign(sumsigns),0)
|
|
|
130
|
|
|
131 DEresults <- data.frame(DE=DE,"DE.fishercomb"=ifelse(fishcomb$adjpval<=alpha,1,0),"DE.invnorm"=ifelse(invnormcomb$adjpval<=alpha,1,0))
|
|
|
132
|
|
|
133 unionDE <- unique(c(fishcomb$DEindices,invnormcomb$DEindices))
|
|
|
134 FC.selecDE <- data.frame(DEresults[unionDE,],FC[unionDE,],signFC=commonsgnFC[unionDE])
|
|
|
135 keepDE <- FC.selecDE[which(abs(FC.selecDE$signFC)==1),]
|
|
|
136
|
|
|
137 fishcomb_de <- rownames(keepDE)[which(keepDE[,"DE.fishercomb"]==1)]
|
|
|
138 invnorm_de <- rownames(keepDE)[which(keepDE[,"DE.invnorm"]==1)]
|
|
|
139 indstudy_de = list()
|
|
|
140 for (i in 1:length(listfiles)) {
|
|
|
141 currentIndstudy_de = rownames(keepDE)[which(keepDE[,i]==1)]
|
|
|
142 indstudy_de[[listfilenames[i]]]=currentIndstudy_de
|
|
|
143 }
|
|
|
144
|
|
|
145 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de)
|
|
|
146 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de)
|
|
|
147
|
|
|
148 #conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],Fishercomb=DEresults[["DE.fishercomb"]],Invnormcomb=DEresults[["DE.invnorm"]],sign=commonsgnFC)
|
|
|
149 conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC)
|
|
|
150 #write DE outputfile
|
|
|
151 write.table(conflits, outputfile,sep="\t",,row.names=FALSE)
|
|
|
152 library(VennDiagram)
|
|
|
153 DE_num=apply(DEresults, 2, FUN=function(x) which(x==1))
|
|
|
154 venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6)
|
|
|
155 temp.venn.plot = file.path( html.files.path, paste("venn.png"))
|
|
|
156 png(temp.venn.plot,width=500,height=500)
|
|
|
157 grid.draw(venn.plot)
|
|
|
158 dev.off()
|
|
|
159
|
|
|
160 library(jsonlite)
|
|
|
161 matrixConflits=as.matrix(conflits)
|
|
|
162 datajson=toJSON(matrixConflits,pretty = TRUE)
|
|
|
163 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE)
|
|
|
164 summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE)
|
|
|
165
|
|
|
166
|
|
|
167 #vennsplit=strsplit(result.venn,split="/")[[1]]
|
|
|
168 #venn=paste0("./",vennsplit[length(vennsplit)])
|
|
|
169
|
|
|
170
|
|
|
171 vennFilename="venn.png"
|
|
|
172 vennFile=file.path(html.files.path,vennFilename)
|
|
|
173 htmlfile=readChar(result.template, file.info(result.template)$size)
|
|
|
174 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
|
|
|
175 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE)
|
|
|
176 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE)
|
|
|
177 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
|
|
|
178 write(htmlfile,result.html)
|
|
|
179
|
|
|
180 #library(VennDiagram)
|
|
|
181 #flog.threshold(ERROR)
|
|
|
182 #
|
|
|
183 ##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")
|
|
|
184 #dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
185 #
|
|
|
186 #showVenn<-function(liste,file)
|
|
|
187 #{
|
|
|
188 # venn.plot<-venn.diagram(x = liste,
|
|
|
189 # filename = vennFilename, col = "black",
|
|
|
190 # fill = 1:length(liste)+1,
|
|
|
191 # margin=0.05, alpha = 0.6,imagetype = "png")
|
|
|
192 ## png(file);
|
|
|
193 ## grid.draw(venn.plot);
|
|
|
194 ## dev.off();
|
|
|
195 #
|
|
|
196 #}
|
|
|
197 #
|
|
|
198 #l=list()
|
|
|
199 #for(i in 1:length(esets))
|
|
|
200 #{
|
|
|
201 # l[[paste("study",i,sep="")]]<-res[[i]]
|
|
|
202 #}
|
|
|
203 #l[["Meta"]]=res[[length(res)-1]]
|
|
|
204 #showVenn(l,vennFile)
|
|
|
205 #file.copy(vennFilename,result.path)
|
|
|
206
|
|
|
207
|
|
|
208 #writeLines( c("<h2>Venn Plot</h2>"), file.conn)
|
|
|
209 #writeLines( c("<img src='venn.png'><br/><br/>"), file.conn)
|
|
|
210 #writeLines( c("</body></html>"), file.conn)
|
|
|
211 #close(file.conn)
|
|
|
212 #print("passe6")
|
|
|
213 #sink(NULL)
|
