|
6
|
1 #!/usr/bin/env Rscript
|
|
|
2 # setup R error handling to go to stderr
|
|
|
3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
|
|
|
4
|
|
|
5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
|
|
|
6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
|
|
|
7
|
|
|
8 library("optparse")
|
|
|
9
|
|
|
10 ##### Read options
|
|
|
11 option_list=list(
|
|
|
12 make_option("--rdatainput",type="character",default="NULL",help="rdata object containing eset object"),
|
|
|
13 make_option("--conditions",type="character",default=NULL,help="Text file summarizing conditions of the experiment (required)"),
|
|
|
14 make_option("--selectcondition1",type="character",default=NULL,help="log2 transformation"),
|
|
|
15 # make_option("--condition1",type="character",default=NULL,help="A table containing the expression data"),
|
|
|
16 make_option("--selectcondition2",type="character",default="NULL",help="rdata object containing eset object"),
|
|
|
17 # make_option("--condition2",type="character",default=NULL,help="Text file summarizing conditions of the experiment (required)"),
|
|
|
18 make_option("--nbresult",type="character",default=NULL,help="number of result displayed results"),
|
|
|
19 make_option("--rdataoutput",type="character",default="NULL",help="output rdata object containing eset object"),
|
|
|
20 make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
|
|
|
21 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
|
|
|
22 make_option("--tabularoutput",type="character",default=NULL,help="Output text file"),
|
|
|
23 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
|
|
|
24
|
|
|
25
|
|
|
26 );
|
|
|
27
|
|
|
28 opt_parser = OptionParser(option_list=option_list);
|
|
|
29 opt = parse_args(opt_parser);
|
|
|
30
|
|
|
31 if(is.null(opt$rdatainput)){
|
|
|
32 print_help(opt_parser)
|
|
|
33 stop("rData input required.", call.=FALSE)
|
|
|
34 }
|
|
|
35
|
|
|
36 if(is.null(opt$conditions)){
|
|
|
37 print_help(opt_parser)
|
|
|
38 stop("conditions input required.", call.=FALSE)
|
|
|
39 }
|
|
|
40
|
|
|
41
|
|
|
42 #loading libraries
|
|
|
43 suppressPackageStartupMessages(require(GEOquery))
|
|
|
44
|
|
|
45 suppressPackageStartupMessages(require(Biobase))
|
|
|
46 suppressPackageStartupMessages(require(GEOquery))
|
|
|
47 suppressPackageStartupMessages(require(GEOmetadb))
|
|
|
48 suppressPackageStartupMessages(require(limma))
|
|
|
49 suppressPackageStartupMessages(require(jsonlite))
|
|
|
50 suppressPackageStartupMessages(require(affy))
|
|
|
51 suppressPackageStartupMessages(require(dplyr))
|
|
2
|
52
|
|
|
53 cargs<-commandArgs()
|
|
|
54 cargs<-cargs[(which(cargs=="--args")+1):length(cargs)]
|
|
|
55 nbargs=length(cargs)
|
|
|
56
|
|
6
|
57 load(opt$rdatainput)
|
|
|
58 targetFile=opt$conditions
|
|
|
59 condition1Name=opt$selectcondition1
|
|
|
60 #condition1=opt$condition1
|
|
|
61 condition2Name=opt$selectcondition2
|
|
|
62 #condition2=opt$condition2
|
|
|
63 nbresult=opt$nbresult
|
|
|
64 result_export_eset=opt$rdataoutput
|
|
|
65 result=opt$htmloutput
|
|
|
66 result.path=opt$htmloutputpath
|
|
|
67 result.tabular=opt$tabularoutput
|
|
|
68 result.template=opt$htmltemplate
|
|
2
|
69
|
|
|
70 #file.copy(targetFile,"./targetFile.txt")
|
|
|
71
|
|
6
|
72 targets <- read.table(targetFile,sep="\t",stringsAsFactors=FALSE)
|
|
2
|
73
|
|
6
|
74 #condition1_tmp <- strsplit(condition1,",")
|
|
|
75 condition1 <-targets[which(targets$V2==condition1Name),1]
|
|
|
76
|
|
|
77 #condition2_tmp <- strsplit(condition2,",")
|
|
|
78 #condition2<-unlist(condition2_tmp)
|
|
|
79 condition2 <-targets[which(targets$V2==condition2Name),1]
|
|
2
|
80
|
|
|
81 conditions=c(condition1,condition2)
|
|
|
82
|
|
|
83 #nbresult=1000
|
|
|
84 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
85
|
|
6
|
86 eset=eset[,which(rownames(eset@phenoData@data) %in% conditions)]
|
|
2
|
87
|
|
6
|
88 rownames(eset@phenoData@data)
|
|
|
89 which(rownames(eset@phenoData@data) %in% conditions)
|
|
2
|
90 #condition1Name=make.names(condition1Name)
|
|
|
91 #condition2Name=make.names(condition2Name)
|
|
|
92 #condition1Name=gsub("_","",condition1Name)
|
|
|
93 #condition2Name=gsub("_","",condition2Name)
|
|
|
94 #condition1Name
|
|
|
95 #condition2Name
|
|
|
96
|
|
|
97 eset@phenoData@data$source_name_ch1=""
|
|
|
98 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition1)]=condition1Name
|
|
|
99 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition2)]=condition2Name
|
|
|
100 #condition1Name
|
|
|
101 #condition2Name
|
|
|
102
|
|
|
103 condNames=paste0("G",as.numeric(as.character(pData(eset)["source_name_ch1"][,1])!=condition1Name))
|
|
|
104 #condNames=make.names(targets[,2])
|
|
|
105 #condNames=gsub("_","",condNames)
|
|
|
106 f <- as.factor(condNames)
|
|
|
107 #eset$description <- factors
|
|
|
108 design <- model.matrix(~ 0+f)
|
|
|
109 design
|
|
|
110
|
|
|
111 colnames(design) <- levels(f)
|
|
|
112 colnames(design)
|
|
|
113 fit <- lmFit(eset, design)
|
|
|
114 fit
|
|
|
115 #cont.matrix <- makeContrasts(C1=paste0(condition1Name,"-",condition2Name), levels=design)
|
|
|
116 cont.matrix <- makeContrasts(G0-G1, levels=design)
|
|
|
117 cont.matrix
|
|
|
118 fit2 <- contrasts.fit(fit, cont.matrix)
|
|
|
119 fit2 <- eBayes(fit2)
|
|
|
120 fit2
|
|
|
121 tT <- topTable(fit2, adjust="fdr", sort.by="B", number=nbresult)
|
|
|
122
|
|
|
123 #head(exprs(eset))
|
|
|
124
|
|
|
125 gpl <- annotation(eset)
|
|
|
126 if (substr(x = gpl,1,3)!="GPL"){
|
|
|
127 #if the annotation info does not start with "GPL" we retrieve the correspondin GPL annotation
|
|
|
128 mapping=read.csv("/galaxy-tools/transcriptomics/db/gplToBioc.csv",stringsAsFactors=FALSE)
|
|
|
129 gpl=mapping[which(mapping$bioc_package==annotation(eset)),]$gpl
|
|
|
130 gpl=gpl[1]
|
|
|
131
|
|
|
132 annotation(eset)=gpl
|
|
|
133
|
|
|
134 platf <- getGEO(gpl, AnnotGPL=TRUE)
|
|
|
135 ncbifd <- data.frame(attr(dataTable(platf), "table"))
|
|
|
136
|
|
|
137 fData(eset)["ID"]=row.names(fData(eset))
|
|
|
138 fData(eset)=merge(x=fData(eset),y=ncbifd,all.x = TRUE, by = "ID")
|
|
|
139 colnames(fData(eset))[4]="ENTREZ_GENE_ID"
|
|
|
140 row.names(fData(eset))=fData(eset)[,"ID"]
|
|
|
141
|
|
|
142 tT <- add_rownames(tT, "ID")
|
|
|
143
|
|
|
144 } else {
|
|
|
145
|
|
|
146 gpl <- annotation(eset)
|
|
|
147 platf <- getGEO(gpl, AnnotGPL=TRUE)
|
|
|
148 ncbifd <- data.frame(attr(dataTable(platf), "table"))
|
|
|
149
|
|
|
150 if (!("ID" %in% colnames(tT))){
|
|
|
151 tT <- add_rownames(tT, "ID")}
|
|
|
152
|
|
|
153 }
|
|
|
154
|
|
|
155 tT <- merge(tT, ncbifd, by="ID")
|
|
|
156 tT <- tT[order(tT$P.Value), ]
|
|
|
157 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"))
|
|
|
158 tT<-format(tT, digits=2, nsmall=2)
|
|
|
159 head(tT)
|
|
|
160 colnames(tT)=gsub(pattern = "\\.",replacement = "_",colnames(tT))
|
|
|
161 matrixtT=as.matrix(tT)
|
|
|
162 datajson=toJSON(matrixtT,pretty = TRUE)
|
|
|
163
|
|
|
164 htmlfile=readChar(result.template, file.info(result.template)$size)
|
|
|
165 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
|
|
|
166 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
|
|
|
167
|
|
|
168 boxplot="boxplot.png"
|
|
|
169 png(boxplot,width=800,height = 400)
|
|
|
170 par(mar=c(7,5,1,1))
|
|
|
171 boxplot(exprs(eset),las=2,outline=FALSE)
|
|
|
172 dev.off()
|
|
|
173 htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOT###",replacement = boxplot, fixed = TRUE)
|
|
|
174 file.copy(boxplot,result.path)
|
|
|
175
|
|
|
176 histopvalue="histopvalue.png"
|
|
|
177
|
|
|
178 png(histopvalue,width=800,height = 400)
|
|
|
179 par(mfrow=c(1,2))
|
|
|
180 hist(fit2$F.p.value,nclass=100,main="Histogram of p-values", xlab="p-values",ylab="frequency")
|
|
|
181 volcanoplot(fit2,coef=1,highlight=10,main="Volcano plot")
|
|
|
182 htmlfile=gsub(x=htmlfile,pattern = "###HIST###",replacement = histopvalue, fixed = TRUE)
|
|
|
183 dev.off()
|
|
|
184 file.copy(histopvalue,result.path)
|
|
|
185
|
|
|
186 #write.table(tolower(c(condition1Name,condition2Name)),quote = FALSE,col.names = FALSE, row.names=FALSE,file=result_export_conditions)
|
|
|
187 saveConditions=c(condition1Name,condition2Name)
|
|
|
188 save(eset,saveConditions,file=result_export_eset)
|
|
|
189 write.table(x=tT[,-1],file=result.tabular,quote=FALSE,row.names=FALSE,col.names=TRUE,sep="\t")
|
|
|
190 write(htmlfile,result)
|
|
|
191
|
