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