|
2
|
1 library(Biobase)
|
|
|
2 library(GEOquery)
|
|
|
3 library(GEOmetadb)
|
|
|
4 library(limma)
|
|
|
5 library(jsonlite)
|
|
|
6 library(affy)
|
|
|
7 library(dplyr)
|
|
|
8
|
|
|
9 cargs<-commandArgs()
|
|
|
10 cargs<-cargs[(which(cargs=="--args")+1):length(cargs)]
|
|
|
11 nbargs=length(cargs)
|
|
|
12
|
|
|
13 load(cargs[[nbargs-13]])
|
|
|
14 targetFile=cargs[[nbargs-12]]
|
|
|
15 condition1Name=cargs[[nbargs-11]]
|
|
|
16 condition1=cargs[[nbargs-10]]
|
|
|
17 condition2Name=cargs[[nbargs-9]]
|
|
|
18 condition2=cargs[[nbargs-8]]
|
|
|
19 nbresult=cargs[[nbargs-7]]
|
|
|
20 result_export_eset=cargs[[nbargs-6]]
|
|
|
21 result=cargs[[nbargs-5]]
|
|
|
22 result.path=cargs[[nbargs-4]]
|
|
|
23 result.tabular=cargs[[nbargs-3]]
|
|
|
24 result.template=cargs[[nbargs-2]]
|
|
|
25
|
|
|
26 #file.copy(targetFile,"./targetFile.txt")
|
|
|
27
|
|
|
28 condition1_tmp <- strsplit(condition1,",")
|
|
|
29 condition1 <-unlist(condition1_tmp)
|
|
|
30
|
|
|
31 condition2_tmp <- strsplit(condition2,",")
|
|
|
32 condition2<-unlist(condition2_tmp)
|
|
|
33
|
|
|
34 conditions=c(condition1,condition2)
|
|
|
35
|
|
|
36 #nbresult=1000
|
|
|
37 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
38
|
|
|
39 targets <- read.table(targetFile,sep="\t")
|
|
|
40
|
|
|
41 eset=eset[,which(rownames(eset@phenoData@data) %in% conditions)]
|
|
|
42 #condition1Name=make.names(condition1Name)
|
|
|
43 #condition2Name=make.names(condition2Name)
|
|
|
44 #condition1Name=gsub("_","",condition1Name)
|
|
|
45 #condition2Name=gsub("_","",condition2Name)
|
|
|
46 #condition1Name
|
|
|
47 #condition2Name
|
|
|
48
|
|
|
49
|
|
|
50 eset@phenoData@data$source_name_ch1=""
|
|
|
51 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition1)]=condition1Name
|
|
|
52 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition2)]=condition2Name
|
|
|
53 #condition1Name
|
|
|
54 #condition2Name
|
|
|
55
|
|
|
56 condNames=paste0("G",as.numeric(as.character(pData(eset)["source_name_ch1"][,1])!=condition1Name))
|
|
|
57 #condNames=make.names(targets[,2])
|
|
|
58 #condNames=gsub("_","",condNames)
|
|
|
59
|
|
|
60 f <- as.factor(condNames)
|
|
|
61 #eset$description <- factors
|
|
|
62 design <- model.matrix(~ 0+f)
|
|
|
63 design
|
|
|
64
|
|
|
65 colnames(design) <- levels(f)
|
|
|
66 colnames(design)
|
|
|
67 fit <- lmFit(eset, design)
|
|
|
68 fit
|
|
|
69 #cont.matrix <- makeContrasts(C1=paste0(condition1Name,"-",condition2Name), levels=design)
|
|
|
70 cont.matrix <- makeContrasts(G0-G1, levels=design)
|
|
|
71 cont.matrix
|
|
|
72 fit2 <- contrasts.fit(fit, cont.matrix)
|
|
|
73 fit2 <- eBayes(fit2)
|
|
|
74 fit2
|
|
|
75 tT <- topTable(fit2, adjust="fdr", sort.by="B", number=nbresult)
|
|
|
76
|
|
|
77 #head(exprs(eset))
|
|
|
78
|
|
|
79 gpl <- annotation(eset)
|
|
|
80 if (substr(x = gpl,1,3)!="GPL"){
|
|
|
81 #if the annotation info does not start with "GPL" we retrieve the correspondin GPL annotation
|
|
|
82 mapping=read.csv("/galaxy-tools/transcriptomics/db/gplToBioc.csv",stringsAsFactors=FALSE)
|
|
|
83 gpl=mapping[which(mapping$bioc_package==annotation(eset)),]$gpl
|
|
|
84 gpl=gpl[1]
|
|
|
85
|
|
|
86 annotation(eset)=gpl
|
|
|
87
|
|
|
88 platf <- getGEO(gpl, AnnotGPL=TRUE)
|
|
|
89 ncbifd <- data.frame(attr(dataTable(platf), "table"))
|
|
|
90
|
|
|
91 fData(eset)["ID"]=row.names(fData(eset))
|
|
|
92 fData(eset)=merge(x=fData(eset),y=ncbifd,all.x = TRUE, by = "ID")
|
|
|
93 colnames(fData(eset))[4]="ENTREZ_GENE_ID"
|
|
|
94 row.names(fData(eset))=fData(eset)[,"ID"]
|
|
|
95
|
|
|
96 tT <- add_rownames(tT, "ID")
|
|
|
97
|
|
|
98 } else {
|
|
|
99
|
|
|
100 gpl <- annotation(eset)
|
|
|
101 platf <- getGEO(gpl, AnnotGPL=TRUE)
|
|
|
102 ncbifd <- data.frame(attr(dataTable(platf), "table"))
|
|
|
103
|
|
|
104 if (!("ID" %in% colnames(tT))){
|
|
|
105 tT <- add_rownames(tT, "ID")}
|
|
|
106
|
|
|
107 }
|
|
|
108
|
|
|
109 tT <- merge(tT, ncbifd, by="ID")
|
|
|
110 tT <- tT[order(tT$P.Value), ]
|
|
|
111 tT <- subset(tT, select=c("Platform_SPOTID","ID","adj.P.Val","P.Value","t","B","logFC","Gene.symbol","Gene.title","Gene.ID","Chromosome.annotation","GO.Function.ID"))
|
|
|
112 tT<-format(tT, digits=2, nsmall=2)
|
|
|
113 head(tT)
|
|
|
114 colnames(tT)=gsub(pattern = "\\.",replacement = "_",colnames(tT))
|
|
|
115 matrixtT=as.matrix(tT)
|
|
|
116 datajson=toJSON(matrixtT,pretty = TRUE)
|
|
|
117
|
|
|
118 htmlfile=readChar(result.template, file.info(result.template)$size)
|
|
|
119 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
|
|
|
120 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
121
|
|
|
122 boxplot="boxplot.png"
|
|
|
123 png(boxplot,width=800,height = 400)
|
|
|
124 par(mar=c(7,5,1,1))
|
|
|
125 boxplot(exprs(eset),las=2,outline=FALSE)
|
|
|
126 dev.off()
|
|
|
127 htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOT###",replacement = boxplot, fixed = TRUE)
|
|
|
128 file.copy(boxplot,result.path)
|
|
|
129
|
|
|
130 histopvalue="histopvalue.png"
|
|
|
131
|
|
|
132 png(histopvalue,width=800,height = 400)
|
|
|
133 par(mfrow=c(1,2))
|
|
|
134 hist(fit2$F.p.value,nclass=100,main="Histogram of p-values", xlab="p-values",ylab="frequency")
|
|
|
135 volcanoplot(fit2,coef=1,highlight=10,main="Volcano plot")
|
|
|
136 htmlfile=gsub(x=htmlfile,pattern = "###HIST###",replacement = histopvalue, fixed = TRUE)
|
|
|
137 dev.off()
|
|
|
138 file.copy(histopvalue,result.path)
|
|
|
139
|
|
|
140 #write.table(tolower(c(condition1Name,condition2Name)),quote = FALSE,col.names = FALSE, row.names=FALSE,file=result_export_conditions)
|
|
|
141 saveConditions=c(condition1Name,condition2Name)
|
|
|
142 save(eset,saveConditions,file=result_export_eset)
|
|
|
143 write.table(x=tT[,-1],file=result.tabular,quote=FALSE,row.names=FALSE,col.names=TRUE,sep="\t")
|
|
|
144 write(htmlfile,result)
|
|
|
145
|
