|
2
|
1 library(metaMA)
|
|
|
2 library(affy)
|
|
|
3 library(annaffy)
|
|
|
4 library(VennDiagram)
|
|
|
5 library(GEOquery)
|
|
|
6
|
|
|
7 cargs<-commandArgs()
|
|
|
8 cargs<-cargs[(which(cargs=="--args")+1):length(cargs)]
|
|
|
9
|
|
|
10 nbargs=length(cargs)
|
|
|
11 rdataList=list()
|
|
|
12 condition1List=list()
|
|
|
13 condition2List=list()
|
|
|
14
|
|
|
15 for (i in 1:(nbargs-5))
|
|
|
16 {
|
|
|
17 Rdata=cargs[[i]]
|
|
|
18 #condition1=cargs[[i+1]]
|
|
|
19 #condition2=cargs[[i+2]]
|
|
|
20 load(Rdata)
|
|
|
21
|
|
|
22 rdataList=c(rdataList,(eset))
|
|
|
23 #condition1List=c(condition1List,condition1)
|
|
|
24 #condition2List=c(condition2List,condition2)
|
|
|
25 condition1List=c(condition1List,saveConditions[1])
|
|
|
26 condition2List=c(condition2List,saveConditions[2])
|
|
|
27
|
|
|
28 }
|
|
|
29
|
|
|
30 #tables<-cargs[[1]]
|
|
|
31 #tech<-cargs[[2]]
|
|
|
32 result.html<-cargs[[nbargs-4]]
|
|
|
33 result.path<-cargs[[nbargs-3]]
|
|
|
34 #result.venn<-cargs[[nbargs-3]]
|
|
|
35 result.template<-cargs[[nbargs-2]]
|
|
|
36
|
|
|
37 #sink("/dev/null")
|
|
|
38 #dir.create(temp.files.path,recursive=TRUE)
|
|
|
39 #file.conn=file(diag.html,open="w")
|
|
|
40
|
|
|
41 #writeLines(c("<html><body bgcolor='lightgray'>"),file.conn)
|
|
|
42
|
|
|
43 showVenn<-function(res,file)
|
|
|
44 {
|
|
|
45 venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),
|
|
|
46 filename = NULL, col = "black",
|
|
|
47 fill = c(1:(length(res)-2)),
|
|
|
48 margin=0.05, alpha = 0.6)
|
|
|
49 jpeg(file)
|
|
|
50 grid.draw(venn.plot)
|
|
|
51 dev.off()
|
|
|
52 }
|
|
|
53
|
|
|
54
|
|
|
55 library("org.Hs.eg.db")
|
|
|
56 x <- org.Hs.egUNIGENE
|
|
|
57 mapped_genes <- mappedkeys(x)
|
|
|
58 link <- as.list(x[mapped_genes])
|
|
|
59
|
|
|
60 probe2unigene<-function(expset){
|
|
|
61 #construction of the map probe->unigene
|
|
|
62 probes=rownames(exprs(expset))
|
|
|
63 gene_id=fData(expset)[probes,"ENTREZ_GENE_ID"]
|
|
|
64 unigene=link[gene_id]
|
|
|
65 names(unigene)<-probes
|
|
|
66 probe_unigene=unigene
|
|
|
67 }
|
|
|
68
|
|
|
69 unigene2probe<-function(map)
|
|
|
70 {
|
|
|
71 suppressWarnings(x <- cbind(unlist(map), names(map)))
|
|
|
72 unigene_probe=split(x[,2], x[,1])
|
|
|
73 }
|
|
|
74
|
|
|
75 convert2metaMA<-function(listStudies,mergemeth=mean)
|
|
|
76 {
|
|
|
77 if (!(class(listStudies) %in% c("list"))) {
|
|
|
78 stop("listStudies must be a list")
|
|
|
79 }
|
|
|
80 conv_unigene=lapply(listStudies,
|
|
|
81 FUN=function(x) unigene2probe(probe2unigene(x)))
|
|
|
82
|
|
|
83 id=lapply(conv_unigene,names)
|
|
|
84 inter=Reduce(intersect,id)
|
|
|
85 if(length(inter)<=0){stop("no common genes")}
|
|
|
86 print(paste(length(inter),"genes in common"))
|
|
|
87 esets=lapply(1:length(listStudies),FUN=function(i){
|
|
|
88 l=lapply(conv_unigene[[i]][inter],
|
|
|
89 FUN=function(x) exprs(listStudies[[i]])[x,,drop=TRUE])
|
|
|
90 esetsgr=t(sapply(l,FUN=function(ll) if(is.null(dim(ll))){ll}
|
|
|
91 else{apply(ll,2,mergemeth)}))
|
|
|
92 esetsgr
|
|
|
93 })
|
|
|
94 return(list(esets=esets,conv.unigene=conv_unigene))
|
|
|
95 }
|
|
|
96
|
|
|
97 normalization<-function(data){
|
|
|
98 ex <- exprs(data)
|
|
|
99 qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
|
|
|
100 LogC <- (qx[5] > 100) ||
|
|
|
101 (qx[6]-qx[1] > 50 && qx[2] > 0) ||
|
|
|
102 (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
|
|
|
103 if (LogC) { ex[which(ex <= 0)] <- NaN
|
|
|
104 return (log2(ex)) } else {
|
|
|
105 return (ex)
|
|
|
106 }
|
|
|
107 }
|
|
|
108
|
|
|
109
|
|
|
110
|
|
|
111
|
|
|
112 filterCondition<-function(gset,condition1, condition2){
|
|
|
113 selected=c(which((tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition1)),
|
|
|
114 which(tolower(as.character(pData(gset)["source_name_ch1"][,1]))==condition2))
|
|
|
115
|
|
|
116 return(gset[,selected])
|
|
|
117 }
|
|
|
118
|
|
|
119 rdatalist <- lapply(rdataList, FUN=function(datalist) normalization(datalist))
|
|
|
120
|
|
|
121 classes=list()
|
|
|
122 filteredRdataList=list()
|
|
|
123 for (i in 1:length(rdatalist))
|
|
|
124 {
|
|
|
125 currentData=rdataList[[i]]
|
|
|
126 currentCondition1=condition1List[[i]]
|
|
|
127 currentCondition2=condition2List[[i]]
|
|
|
128 #currentData=filterCondition(currentData,currentCondition1,currentCondition2)
|
|
|
129 currentClasses=as.numeric(tolower(as.character(pData(currentData)["source_name_ch1"][,1]))==currentCondition1)
|
|
|
130 filteredRdataList=c(filteredRdataList,currentData)
|
|
|
131 classes=c(classes,list(currentClasses))
|
|
|
132 #write(file="~/galaxy-modal/classes.txt",classes)
|
|
|
133 }
|
|
|
134
|
|
|
135 #rdataList=filteredRdataList
|
|
|
136 conv=convert2metaMA(rdataList)
|
|
|
137 esets=conv$esets
|
|
|
138 conv_unigene=conv$conv.unigene
|
|
|
139
|
|
|
140 #write(file="~/galaxy-modal/esets.txt",length(esets))
|
|
|
141 #write(file="~/galaxy-modal/classes.txt",length(classes))
|
|
|
142 res=pvalcombination(esets=esets,classes=classes,moderated="limma")
|
|
|
143 resIDDIRR=IDDIRR(res$Meta,res$AllIndStudies)
|
|
|
144 length(res$Meta)
|
|
|
145 Hs.Meta=rownames(esets[[1]])[res$Meta]
|
|
|
146 origId.Meta=lapply(conv_unigene,FUN=function(vec) as.vector(unlist(vec[Hs.Meta])))
|
|
|
147
|
|
|
148 gpllist <- lapply(rdataList, FUN=function(ann) annotation(ann))
|
|
|
149 platflist <- lapply(gpllist, FUN=function(gpl) getGEO(gpl, AnnotGPL=TRUE))
|
|
|
150 ncbifdlist <- lapply(platflist, FUN=function(data) data.frame(attr(dataTable(data), "table")))
|
|
|
151 ncbifdresult=lapply(1:length(origId.Meta), FUN=function(i) ncbifdlist[[i]][which(ncbifdlist[[i]]$ID %in% origId.Meta[[i]]),])
|
|
|
152 ncbidfannot=do.call(rbind,ncbifdresult)
|
|
|
153 ncbidfannot <- subset(ncbidfannot, select=c("Platform_SPOTID","ID","Gene.symbol","Gene.title","Gene.ID","Chromosome.annotation","GO.Function.ID"))
|
|
|
154
|
|
|
155 library(jsonlite)
|
|
|
156 matrixncbidfannot=as.matrix(ncbidfannot)
|
|
|
157 datajson=toJSON(matrixncbidfannot,pretty = TRUE)
|
|
|
158 summaryjson=toJSON(as.matrix(t(resIDDIRR)),pretty = TRUE)
|
|
|
159
|
|
|
160
|
|
|
161 #vennsplit=strsplit(result.venn,split="/")[[1]]
|
|
|
162 #venn=paste0("./",vennsplit[length(vennsplit)])
|
|
|
163
|
|
|
164
|
|
|
165 vennFilename="venn.png"
|
|
|
166 vennFile=file.path(result.path,vennFilename)
|
|
|
167 htmlfile=readChar(result.template, file.info(result.template)$size)
|
|
|
168 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
|
|
|
169 htmlfile=gsub(x=htmlfile,pattern = "###SUMMARYJSON###",replacement = summaryjson, fixed = TRUE)
|
|
|
170 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
|
|
|
171 write(htmlfile,result.html)
|
|
|
172
|
|
|
173 library(VennDiagram)
|
|
|
174 flog.threshold(ERROR)
|
|
|
175
|
|
|
176 #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")
|
|
|
177 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
178
|
|
|
179 showVenn<-function(liste,file)
|
|
|
180 {
|
|
|
181 venn.plot<-venn.diagram(x = liste,
|
|
|
182 filename = vennFilename, col = "black",
|
|
|
183 fill = 1:length(liste)+1,
|
|
|
184 margin=0.05, alpha = 0.6,imagetype = "png")
|
|
|
185 # png(file);
|
|
|
186 # grid.draw(venn.plot);
|
|
|
187 # dev.off();
|
|
|
188
|
|
|
189 }
|
|
|
190
|
|
|
191 l=list()
|
|
|
192 for(i in 1:length(esets))
|
|
|
193 {
|
|
|
194 l[[paste("study",i,sep="")]]<-res[[i]]
|
|
|
195 }
|
|
|
196 l[["Meta"]]=res[[length(res)-1]]
|
|
|
197 showVenn(l,vennFile)
|
|
|
198 file.copy(vennFilename,result.path)
|
|
|
199 #l=list()
|
|
|
200 #for(i in 1:length(esets))
|
|
|
201 #{
|
|
|
202 # l[[paste("study",i,sep="")]]<-res[[i]]
|
|
|
203 #}
|
|
|
204 #l[["Meta"]]=res[[length(res)-1]]
|
|
|
205 #showVenn(res,result.venn)
|
|
|
206 #writeLines(c("<h2>Venn diagram</h2>"),file.conn)
|
|
|
207 #writeLines(c("<img src='venn.png'><br/><br/>"),file.conn)
|
|
|
208 #writeLines(c("</body></html>"),file.conn)
|
|
|
209 #close(file.conn) |
