Mercurial > repos > sblanck > smagexp

comparison MetaMA.R @ 2:93451f832736 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author sblanck
date Tue, 21 Mar 2017 10:28:47 -0400
parents
children 328f4031e5d3
comparison
equal deleted inserted replaced
1:f8a2f1fec8ef 2:93451f832736
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)