Mercurial > repos > sblanck > smagexp

comparison AffyQCnormalization.R @ 0:1024245abc70 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
-1:000000000000 0:1024245abc70
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("--input",type="character",default="NULL",help="rdata object containing eset object"),
13 make_option("--normalization",type="character",default=NULL,help="normalization method"),
14 make_option("--nbresult",type="character",default=NULL,help="number of result displayed results"),
15 make_option("--rdataoutput",type="character",default="NULL",help="output rdata object containing eset object"),
16 make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
17 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
18 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
19 );
20
21 opt_parser = OptionParser(option_list=option_list);
22 opt = parse_args(opt_parser);
23
24 if(is.null(opt$input)){
25 print_help(opt_parser)
26 stop("input required.", call.=FALSE)
27 }
28
29 #loading libraries
30
31 suppressPackageStartupMessages(require(Biobase))
32 suppressPackageStartupMessages(require(GEOquery))
33 suppressPackageStartupMessages(require(GEOmetadb))
34 suppressPackageStartupMessages(require(limma))
35 suppressPackageStartupMessages(require(jsonlite))
36 suppressPackageStartupMessages(require(affy))
37 suppressPackageStartupMessages(require(affyPLM))
38 suppressPackageStartupMessages(require(dplyr))
39
40 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
41
42 celList=vector()
43 celFileNameList=vector()
44
45 for (i in 1:length(listInput))
46 {
47 inputFileInfo <- unlist( strsplit( listInput[i], ';' ) )
48 celList=c(celList,inputFileInfo[1])
49 celFileNameList=c(celFileNameList,inputFileInfo[2])
50 }
51
52
53 normalization=opt$normalization
54 result_export_eset=opt$rdataoutput
55 result=opt$htmloutput
56 result.path=opt$htmloutputpath
57 result.template=opt$htmltemplate
58
59 dir.create(result.path, showWarnings = TRUE, recursive = TRUE)
60 for(i in 1:length(celList))
61 {
62 file.copy(celList[i],paste0("./",celFileNameList[i]))
63 }
64
65 data <- ReadAffy(filenames=celFileNameList, celfile.path=".")
66 htmlfile=readChar(result.template, file.info(result.template)$size)
67
68 boxplot="boxplot.png"
69 png(boxplot,width=800,height = 400)
70 par(mar=c(7,5,1,1))
71 boxplot(data,las=2,outline=FALSE)
72 dev.off()
73 htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOT###",replacement = boxplot, fixed = TRUE)
74 file.copy(boxplot,result.path)
75
76 images="images.png"
77 nblines=length(celList)%/%4 + as.numeric((length(celList)%%4)!=0)
78 png(images,width=800,height = 200*nblines)
79 par(mfrow=c(nblines,4))
80 image(data)
81 dev.off()
82 htmlfile=gsub(x=htmlfile,pattern = "###IMAGES###",replacement = images, fixed = TRUE)
83 file.copy(images,result.path)
84
85
86 plotMA="plotMA.png"
87 nblines=length(celList)%/%3 + as.numeric((length(celList)%%3)!=0)
88 png(plotMA,width=800,height =300*nblines )
89 par(mfrow=c(nblines,3))
90 MAplot(data)
91 dev.off()
92 htmlfile=gsub(x=htmlfile,pattern = "###PLOTMA###",replacement = plotMA, fixed = TRUE)
93 file.copy(plotMA,result.path)
94
95
96 if (normalization == "rma") {
97 eset <- rma(data)
98 } else if (normalization == "quantile") {
99 eset = rma(data,background = FALSE,normalize = TRUE)
100 } else if (normalization == "background"){
101 eset = rma(data,background = TRUE ,normalize = FALSE)
102 } else if (normalization == "log2") {
103 eset = rma(data,background = FALSE ,normalize = FALSE)
104 }
105
106
107 boxplotnorm="boxplotnorm.png"
108 png(boxplotnorm,width=800,height = 400)
109 par(mar=c(7,5,1,1))
110 boxplot(data.frame(exprs(eset)),las=2,outline=FALSE)
111 dev.off()
112 htmlfile=gsub(x=htmlfile,pattern = "###BOXPLOTNORM###",replacement = boxplotnorm, fixed = TRUE)
113 file.copy(boxplotnorm,result.path)
114
115 plotMAnorm="plotMAnorm.png"
116 nblines=length(celList)%/%3 + as.numeric((length(celList)%%3)!=0)
117 png(plotMAnorm,width=800,height =300*nblines )
118 par(mfrow=c(nblines,3))
119 #for (i in 1:length(celList)){
120 MAplot(eset)
121 #}
122
123 dev.off()
124 htmlfile=gsub(x=htmlfile,pattern = "###PLOTMANORM###",replacement = plotMAnorm, fixed = TRUE)
125 file.copy(plotMAnorm,result.path)
126 save(eset,file=result_export_eset)
127 write(htmlfile,result)
128