comparison GEOQuery.R @ 2:93451f832736 draft

Uploaded
author sblanck
date Tue, 21 Mar 2017 10:28:47 -0400
parents
children 2c9e44ff68dc
comparison
equal deleted inserted replaced
1:f8a2f1fec8ef 2:93451f832736
1
2 library(GEOquery)
3
4
5 cargs<-commandArgs()
6 cargs<-cargs[(which(cargs=="--args")+1):length(cargs)]
7
8 GEOQueryID<-cargs[[1]]
9 GEOQueryData<-cargs[[2]]
10 GEOQueryRData<-cargs[[3]]
11 conditionFile<-cargs[[4]]
12 transformation<-cargs[[5]]
13
14 data1 = getGEO(GEOQueryID)
15 eset=data1[[1]]
16
17 #check if datas are in log2 space
18 normalization<-function(data){
19 ex <- exprs(data)
20 qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
21 LogC <- (qx[5] > 100) ||
22 (qx[6]-qx[1] > 50 && qx[2] > 0) ||
23 (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
24 if (LogC) { ex[which(ex <= 0)] <- NaN
25 return (log2(ex)) } else {
26 return (ex)
27 }
28 }
29
30 if (transformation=="auto"){
31 exprs(eset)=normalization(eset)
32 } else if (transformation=="yes"){
33 exprs(eset)=log2(exprs(eset))
34 }
35
36 matrixData=exprs(eset)
37 write.table(matrixData,col.names=NA,row.names=TRUE,sep="\t",file=GEOQueryData)
38
39 #Construcion of condition file
40 #if there is data in "source_name_ch1" field, we keep this data as a condition
41 #else we keep the "description" field data.
42 if (length(unique(tolower(pData(data1[[1]])["source_name_ch1"][,1])))>1)
43 {
44 conditions=pData(data1[[1]])["source_name_ch1"]
45 description=paste0(as.vector(pData(data1[[1]])["geo_accession"][,1]), " ",as.vector(pData(data1[[1]])["title"][,1]), " ", as.vector(conditions[,1]))
46 } else
47 {
48 conditions=pData(data1[[1]])["description"]
49 description=paste0(as.vector(pData(data1[[1]])["geo_accession"][,1]), " ",as.vector(pData(data1[[1]])["title"][,1]), " ", as.vector(conditions[,1]))
50 }
51
52 conditions[,1]=tolower(conditions[,1])
53 pData(eset)["source_name_ch1"]=conditions
54
55 write.table(cbind(conditions,description),quote = FALSE,col.names = FALSE, row.names=TRUE,file=conditionFile,sep="\t")
56 save(eset,conditions,file=GEOQueryRData)