Mercurial > repos > sblanck > smagexp
annotate Analyse.R @ 11:f9732f6bf218 draft
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
| author | sblanck |
|---|---|
| date | Thu, 11 May 2017 10:20:35 -0400 |
| parents | 56267e3293b2 |
| children | 7c49b4ec8cff |
| rev | line source |
|---|---|
| 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"), | |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
23 make_option("--htmltemplate",type="character",default=NULL,help="html template)"), |
|
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
24 make_option("--tooldirectory",type="character",default=NULL,help="tool directory)") |
| 6 | 25 |
| 26 | |
| 27 ); | |
| 28 | |
| 29 opt_parser = OptionParser(option_list=option_list); | |
| 30 opt = parse_args(opt_parser); | |
| 31 | |
| 32 if(is.null(opt$rdatainput)){ | |
| 33 print_help(opt_parser) | |
| 34 stop("rData input required.", call.=FALSE) | |
| 35 } | |
| 36 | |
| 37 if(is.null(opt$conditions)){ | |
| 38 print_help(opt_parser) | |
| 39 stop("conditions input required.", call.=FALSE) | |
| 40 } | |
| 41 | |
| 42 | |
| 43 #loading libraries | |
| 44 suppressPackageStartupMessages(require(GEOquery)) | |
| 45 | |
| 46 suppressPackageStartupMessages(require(Biobase)) | |
| 47 suppressPackageStartupMessages(require(GEOquery)) | |
| 48 suppressPackageStartupMessages(require(GEOmetadb)) | |
| 49 suppressPackageStartupMessages(require(limma)) | |
| 50 suppressPackageStartupMessages(require(jsonlite)) | |
| 51 suppressPackageStartupMessages(require(affy)) | |
| 52 suppressPackageStartupMessages(require(dplyr)) | |
| 2 | 53 |
| 6 | 54 load(opt$rdatainput) |
| 55 targetFile=opt$conditions | |
| 56 condition1Name=opt$selectcondition1 | |
| 57 #condition1=opt$condition1 | |
| 58 condition2Name=opt$selectcondition2 | |
| 59 #condition2=opt$condition2 | |
| 60 nbresult=opt$nbresult | |
| 61 result_export_eset=opt$rdataoutput | |
| 62 result=opt$htmloutput | |
| 63 result.path=opt$htmloutputpath | |
| 64 result.tabular=opt$tabularoutput | |
| 65 result.template=opt$htmltemplate | |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
66 tooldirectory=opt$tooldirectory |
| 2 | 67 |
| 68 #file.copy(targetFile,"./targetFile.txt") | |
| 69 | |
| 6 | 70 targets <- read.table(targetFile,sep="\t",stringsAsFactors=FALSE) |
| 2 | 71 |
| 6 | 72 #condition1_tmp <- strsplit(condition1,",") |
| 73 condition1 <-targets[which(targets$V2==condition1Name),1] | |
| 74 | |
| 75 #condition2_tmp <- strsplit(condition2,",") | |
| 76 #condition2<-unlist(condition2_tmp) | |
| 77 condition2 <-targets[which(targets$V2==condition2Name),1] | |
| 2 | 78 |
| 79 conditions=c(condition1,condition2) | |
| 80 | |
| 81 #nbresult=1000 | |
| 82 dir.create(result.path, showWarnings = TRUE, recursive = FALSE) | |
| 83 | |
| 6 | 84 eset=eset[,which(rownames(eset@phenoData@data) %in% conditions)] |
| 2 | 85 |
| 6 | 86 rownames(eset@phenoData@data) |
| 87 which(rownames(eset@phenoData@data) %in% conditions) | |
| 2 | 88 #condition1Name=make.names(condition1Name) |
| 89 #condition2Name=make.names(condition2Name) | |
| 90 #condition1Name=gsub("_","",condition1Name) | |
| 91 #condition2Name=gsub("_","",condition2Name) | |
| 92 #condition1Name | |
| 93 #condition2Name | |
| 94 | |
| 95 eset@phenoData@data$source_name_ch1="" | |
| 96 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition1)]=condition1Name | |
| 97 eset@phenoData@data$source_name_ch1[which(rownames(eset@phenoData@data) %in% condition2)]=condition2Name | |
| 98 #condition1Name | |
| 99 #condition2Name | |
| 100 | |
| 101 condNames=paste0("G",as.numeric(as.character(pData(eset)["source_name_ch1"][,1])!=condition1Name)) | |
| 102 #condNames=make.names(targets[,2]) | |
| 103 #condNames=gsub("_","",condNames) | |
| 104 f <- as.factor(condNames) | |
| 105 #eset$description <- factors | |
| 106 design <- model.matrix(~ 0+f) | |
| 107 | |
| 108 colnames(design) <- levels(f) | |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
109 #colnames(design) |
| 2 | 110 fit <- lmFit(eset, design) |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
111 |
| 2 | 112 #cont.matrix <- makeContrasts(C1=paste0(condition1Name,"-",condition2Name), levels=design) |
| 113 cont.matrix <- makeContrasts(G0-G1, levels=design) | |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
114 #cont.matrix |
| 2 | 115 fit2 <- contrasts.fit(fit, cont.matrix) |
| 116 fit2 <- eBayes(fit2) | |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
117 |
| 2 | 118 tT <- topTable(fit2, adjust="fdr", sort.by="B", number=nbresult) |
| 119 | |
| 120 #head(exprs(eset)) | |
| 121 | |
| 122 gpl <- annotation(eset) | |
| 123 if (substr(x = gpl,1,3)!="GPL"){ | |
|
10
56267e3293b2
planemo upload commit 0ad899de37d891f387912b7abd63506605a25976-dirty
sblanck
parents:
7
diff
changeset
|
124 #if the annotation info does not start with "GPL" we retrieve the corresponding GPL annotation |
|
11
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
125 print(getwd()) |
|
f9732f6bf218
planemo upload commit 3cad4bb49300f86625d62d40ea7e2ea2ba40e314-dirty
sblanck
parents:
10
diff
changeset
|
126 mapping=read.csv(paste0(tooldirectory,"/gplToBioc.csv"),stringsAsFactors=FALSE) |
| 2 | 127 gpl=mapping[which(mapping$bioc_package==annotation(eset)),]$gpl |
| 128 gpl=gpl[1] | |
| 129 | |
| 130 annotation(eset)=gpl | |
| 131 | |
| 132 platf <- getGEO(gpl, AnnotGPL=TRUE) | |
| 133 ncbifd <- data.frame(attr(dataTable(platf), "table")) | |
| 134 | |
| 135 fData(eset)["ID"]=row.names(fData(eset)) | |
| 136 fData(eset)=merge(x=fData(eset),y=ncbifd,all.x = TRUE, by = "ID") | |
| 137 colnames(fData(eset))[4]="ENTREZ_GENE_ID" | |
| 138 row.names(fData(eset))=fData(eset)[,"ID"] | |
| 139 | |
| 140 tT <- add_rownames(tT, "ID") | |
| 141 | |
| 142 } else { | |
| 143 | |
| 144 gpl <- annotation(eset) | |
| 145 platf <- getGEO(gpl, AnnotGPL=TRUE) | |
| 146 ncbifd <- data.frame(attr(dataTable(platf), "table")) | |
| 147 | |
| 148 if (!("ID" %in% colnames(tT))){ | |
| 149 tT <- add_rownames(tT, "ID")} | |
| 150 | |
| 151 } | |
| 152 | |
| 153 tT <- merge(tT, ncbifd, by="ID") | |
| 154 tT <- tT[order(tT$P.Value), ] | |
| 155 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")) | |
| 156 tT<-format(tT, digits=2, nsmall=2) | |
| 157 head(tT) | |
| 158 colnames(tT)=gsub(pattern = "\\.",replacement = "_",colnames(tT)) | |
| 159 matrixtT=as.matrix(tT) | |
| 160 datajson=toJSON(matrixtT,pretty = TRUE) | |
| 161 | |
| 162 htmlfile=readChar(result.template, file.info(result.template)$size) | |
| 163 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE) | |
| 164 dir.create(result.path, showWarnings = TRUE, recursive = FALSE) | |
| 165 | |
| 166 boxplot="boxplot.png" | |
| 167 png(boxplot,width=800,height = 400) | |
| 168 par(mar=c(7,5,1,1)) | |
| 169 boxplot(exprs(eset),las=2,outline=FALSE) | |
| 170 dev.off() | |
| 171 htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOT###",replacement = boxplot, fixed = TRUE) | |
| 172 file.copy(boxplot,result.path) | |
| 173 | |
| 174 histopvalue="histopvalue.png" | |
| 175 | |
| 176 png(histopvalue,width=800,height = 400) | |
| 177 par(mfrow=c(1,2)) | |
| 178 hist(fit2$F.p.value,nclass=100,main="Histogram of p-values", xlab="p-values",ylab="frequency") | |
| 179 volcanoplot(fit2,coef=1,highlight=10,main="Volcano plot") | |
| 180 htmlfile=gsub(x=htmlfile,pattern = "###HIST###",replacement = histopvalue, fixed = TRUE) | |
| 181 dev.off() | |
| 182 file.copy(histopvalue,result.path) | |
| 183 | |
| 184 #write.table(tolower(c(condition1Name,condition2Name)),quote = FALSE,col.names = FALSE, row.names=FALSE,file=result_export_conditions) | |
| 185 saveConditions=c(condition1Name,condition2Name) | |
| 186 save(eset,saveConditions,file=result_export_eset) | |
| 187 write.table(x=tT[,-1],file=result.tabular,quote=FALSE,row.names=FALSE,col.names=TRUE,sep="\t") | |
| 188 write(htmlfile,result) | |
| 189 |