Mercurial > repos > sblanck > smagexp
diff AffyQCnormalization.R @ 0:1024245abc70 draft
planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit 5974f806f344dbcc384b931492d7f023bfbbe03b
author | sblanck |
---|---|
date | Thu, 22 Feb 2018 08:38:22 -0500 |
parents | |
children | e4e6e583b8d9 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/AffyQCnormalization.R Thu Feb 22 08:38:22 2018 -0500 @@ -0,0 +1,128 @@ +#!/usr/bin/env Rscript +# setup R error handling to go to stderr +options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) + +# we need that to not crash galaxy with an UTF8 error on German LC settings. +loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") + +library("optparse") + +##### Read options +option_list=list( + make_option("--input",type="character",default="NULL",help="rdata object containing eset object"), + make_option("--normalization",type="character",default=NULL,help="normalization method"), + make_option("--nbresult",type="character",default=NULL,help="number of result displayed results"), + make_option("--rdataoutput",type="character",default="NULL",help="output rdata object containing eset object"), + make_option("--htmloutput",type="character",default=NULL,help="Output html report"), + make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"), + make_option("--htmltemplate",type="character",default=NULL,help="html template)") +); + +opt_parser = OptionParser(option_list=option_list); +opt = parse_args(opt_parser); + +if(is.null(opt$input)){ + print_help(opt_parser) + stop("input required.", call.=FALSE) +} + +#loading libraries + +suppressPackageStartupMessages(require(Biobase)) +suppressPackageStartupMessages(require(GEOquery)) +suppressPackageStartupMessages(require(GEOmetadb)) +suppressPackageStartupMessages(require(limma)) +suppressPackageStartupMessages(require(jsonlite)) +suppressPackageStartupMessages(require(affy)) +suppressPackageStartupMessages(require(affyPLM)) +suppressPackageStartupMessages(require(dplyr)) + +listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) ) + +celList=vector() +celFileNameList=vector() + +for (i in 1:length(listInput)) +{ + inputFileInfo <- unlist( strsplit( listInput[i], ';' ) ) + celList=c(celList,inputFileInfo[1]) + celFileNameList=c(celFileNameList,inputFileInfo[2]) +} + + +normalization=opt$normalization +result_export_eset=opt$rdataoutput +result=opt$htmloutput +result.path=opt$htmloutputpath +result.template=opt$htmltemplate + +dir.create(result.path, showWarnings = TRUE, recursive = TRUE) +for(i in 1:length(celList)) +{ + file.copy(celList[i],paste0("./",celFileNameList[i])) +} + +data <- ReadAffy(filenames=celFileNameList, celfile.path=".") +htmlfile=readChar(result.template, file.info(result.template)$size) + +boxplot="boxplot.png" +png(boxplot,width=800,height = 400) +par(mar=c(7,5,1,1)) +boxplot(data,las=2,outline=FALSE) +dev.off() +htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOT###",replacement = boxplot, fixed = TRUE) +file.copy(boxplot,result.path) + +images="images.png" +nblines=length(celList)%/%4 + as.numeric((length(celList)%%4)!=0) +png(images,width=800,height = 200*nblines) +par(mfrow=c(nblines,4)) +image(data) +dev.off() +htmlfile=gsub(x=htmlfile,pattern = "###IMAGES###",replacement = images, fixed = TRUE) +file.copy(images,result.path) + + +plotMA="plotMA.png" +nblines=length(celList)%/%3 + as.numeric((length(celList)%%3)!=0) +png(plotMA,width=800,height =300*nblines ) +par(mfrow=c(nblines,3)) +MAplot(data) +dev.off() +htmlfile=gsub(x=htmlfile,pattern = "###PLOTMA###",replacement = plotMA, fixed = TRUE) +file.copy(plotMA,result.path) + + +if (normalization == "rma") { + eset <- rma(data) +} else if (normalization == "quantile") { + eset = rma(data,background = FALSE,normalize = TRUE) +} else if (normalization == "background"){ + eset = rma(data,background = TRUE ,normalize = FALSE) +} else if (normalization == "log2") { + eset = rma(data,background = FALSE ,normalize = FALSE) +} + + +boxplotnorm="boxplotnorm.png" +png(boxplotnorm,width=800,height = 400) +par(mar=c(7,5,1,1)) +boxplot(data.frame(exprs(eset)),las=2,outline=FALSE) +dev.off() +htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOTNORM###",replacement = boxplotnorm, fixed = TRUE) +file.copy(boxplotnorm,result.path) + +plotMAnorm="plotMAnorm.png" +nblines=length(celList)%/%3 + as.numeric((length(celList)%%3)!=0) +png(plotMAnorm,width=800,height =300*nblines ) +par(mfrow=c(nblines,3)) +#for (i in 1:length(celList)){ +MAplot(eset) +#} + +dev.off() +htmlfile=gsub(x=htmlfile,pattern = "###PLOTMANORM###",replacement = plotMAnorm, fixed = TRUE) +file.copy(plotMAnorm,result.path) +save(eset,file=result_export_eset) +write(htmlfile,result) +