Mercurial > repos > sblanck > smagexp

comparison ImportDataFromMatrix.R @ 8:7f74250a083d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author sblanck
date Wed, 10 May 2017 04:12:26 -0400
parents 93451f832736
children 25b828010ca9
comparison
equal deleted inserted replaced
7:f3c021bdc000 8:7f74250a083d
1 library(Biobase) 1 #!/usr/bin/env Rscript
2 library(GEOquery) 2 # setup R error handling to go to stderr
3 library(GEOmetadb) 3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
4 library(limma)
5 library(jsonlite)
6 library(affy)
7 library(dplyr)
8 library(affyPLM)
9 4
10 cargs<-commandArgs() 5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
11 cargs<-cargs[(which(cargs=="--args")+1):length(cargs)] 6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
12 nbargs=length(cargs)
13 7
14 dataFile=cargs[[nbargs-9]] 8 library("optparse")
15 normalization=cargs[[nbargs-8]] 9
16 conditionsFile=cargs[[nbargs-7]] 10 ##### Read options
17 annotation=cargs[[nbargs-6]] 11 option_list=list(
18 result_export_eset=cargs[[nbargs-5]] 12 make_option("--input",type="character",default="NULL",help="rdata object containing eset object"),
19 result=cargs[[nbargs-4]] 13 make_option("--conditions",type="character",default=NULL,help="Text file summarizing conditions of the experiment (required)"),
20 result.path=cargs[[nbargs-3]] 14 make_option("--normalization",type="character",default=NULL,help="log2 transformation"),
21 result.template=cargs[[nbargs-2]] 15 make_option("--annotations",type="character",default="NULL",help="rdata object containing eset object"),
16 make_option("--rdataoutput",type="character",default="NULL",help="output rdata object containing eset object"),
17 make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
18 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
19 make_option("--htmltemplate",type="character",default="NULL",help="html template)")
20
21
22 );
23
24
25 opt_parser = OptionParser(option_list=option_list);
26 opt = parse_args(opt_parser);
27
28 if(is.null(opt$input)){
29 print_help(opt_parser)
30 stop("input required.", call.=FALSE)
31 }
32
33 if(is.null(opt$conditions)){
34 print_help(opt_parser)
35 stop("conditions input required.", call.=FALSE)
36 }
37
38
39 #loading libraries
40 suppressPackageStartupMessages(require(GEOquery))
41
42 suppressPackageStartupMessages(require(Biobase))
43 suppressPackageStartupMessages(require(GEOquery))
44 suppressPackageStartupMessages(require(GEOmetadb))
45 suppressPackageStartupMessages(require(limma))
46 suppressPackageStartupMessages(require(jsonlite))
47 suppressPackageStartupMessages(require(affy))
48 suppressPackageStartupMessages(require(dplyr))
49 suppressPackageStartupMessages(require(affyPLM))
50
51 dataFile=opt$input
52 normalization=opt$normalization
53 conditionsFile=opt$conditions
54 annotation=opt$annotations
55 result_export_eset=opt$rdataoutput
56 result=opt$htmloutput
57 result.path=opt$htmloutputpath
58 result.template=opt$htmltemplate
22 59
23 dir.create(result.path, showWarnings = TRUE, recursive = FALSE) 60 dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
24 61
25 data=as.matrix(read.table(file = dataFile)) 62 data=as.matrix(read.table(file = dataFile,row.names=1,header=TRUE))
26 conditions=read.table(file=conditionsFile,sep = "\t",row.names=1) 63 conditions=read.table(file=conditionsFile,sep = "\t",row.names=1)
27 htmlfile=readChar(result.template, file.info(result.template)$size) 64 htmlfile=readChar(result.template, file.info(result.template)$size)
28 65
29 colnames(conditions)=c("source_name_ch1","description") 66 colnames(conditions)=c("source_name_ch1","description")
30 phenodata<-new("AnnotatedDataFrame",data=conditions) 67 phenodata<-new("AnnotatedDataFrame",data=conditions)
31 68
69 head(data)
70 conditions
71
32 eset=ExpressionSet(assayData=data,phenoData=phenodata,annotation=annotation) 72 eset=ExpressionSet(assayData=data,phenoData=phenodata,annotation=annotation)
33 73
34 if (normalization == "quantile") { 74 if (normalization == "quantile") {
35 eset <- normalize.ExpressionSet.quantiles(eset, transfn="log") 75 eset <- normalize.ExpressionSet.quantiles(eset, transfn="log2")
36 } else if (normalization == "log2") { 76 } else if (normalization == "log2") {
37 exprs(eset) = log2(exprs(eset)) 77 exprs(eset) = log2(exprs(eset))
38 } 78 }
39 79
40 boxplotnorm="boxplotnorm.png" 80 boxplotnorm="boxplotnorm.png"
47 87
48 plotMAnorm="plotMAnorm.png" 88 plotMAnorm="plotMAnorm.png"
49 nblines=length(colnames(data))%/%3 + as.numeric((length(colnames(data))%%3)!=0) 89 nblines=length(colnames(data))%/%3 + as.numeric((length(colnames(data))%%3)!=0)
50 png(plotMAnorm,width=800,height =300*nblines ) 90 png(plotMAnorm,width=800,height =300*nblines )
51 par(mfrow=c(nblines,3)) 91 par(mfrow=c(nblines,3))
52 #for (i in 1:length(colnames(data))){ 92 ##for (i in 1:length(colnames(data))){
53 MAplot(eset) 93 MAplot(eset)
54 #} 94 #}
55 95
56 dev.off() 96 dev.off()
57 htmlfile=gsub(x=htmlfile,pattern = "###PLOTMANORM###",replacement = plotMAnorm, fixed = TRUE) 97 htmlfile=gsub(x=htmlfile,pattern = "###PLOTMANORM###",replacement = plotMAnorm, fixed = TRUE)