Mercurial > repos > pavanvidem > dexseq
view dexseq/dexseq.R @ 0:7604d324c5aa draft
Uploaded
author | pavanvidem |
---|---|
date | Fri, 28 Aug 2015 08:37:31 -0400 |
parents | |
children |
line wrap: on
line source
## 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. Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") library("DEXSeq") library('getopt') library('rjson') options(stringAsfactors = FALSE, useFancyQuotes = FALSE) args <- commandArgs(trailingOnly = TRUE) #get options, using the spec as defined by the enclosed list. #we read the options from the default: commandArgs(TRUE). spec = matrix(c( 'verbose', 'v', 2, "integer", 'help', 'h', 0, "logical", 'gtf', 'a', 1, "character", 'outfile', 'o', 1, "character", 'report', 'r', 1, "character", 'factors', 'f', 1, "character", 'threads', 'p', 1, "integer", 'fdr', 'c', 1, "double", ), byrow=TRUE, ncol=4); opt = getopt(spec); setwd(opt$report) # if help was asked for print a friendly message # and exit with a non-zero error code if ( !is.null(opt$help) ) { cat(getopt(spec, usage=TRUE)); q(status=1); } trim <- function (x) gsub("^\\s+|\\s+$", "", x) opt$samples <- trim(opt$samples) opt$factors <- trim(opt$factors) parser <- newJSONParser() parser$addData( opt$factors ) factorsList <- parser$getObject() sampleTable<-data.frame() countFiles<-c() factorNames<-c() primaryFactor<-"" for(factor in factorsList){ factorName<-factor[[1]] factorNames<-append(factorNames, paste(factorName,"exon",sep=":")) factorValuesMapList<-factor[[2]] c = length(factorValuesMapList) for (factorValuesMap in factorValuesMapList){ for(files in factorValuesMap){ for(file in files){ if(primaryFactor == "") { countFiles<-append(countFiles,file) } sampleTable[basename(file),factorName]<-paste(c,names(factorValuesMap),sep="_") } } c = c-1 } if(primaryFactor == ""){ primaryFactor <- factorName } } factorNames<-append(factorNames,"exon") factorNames<-append(factorNames,"sample") factorNames<-rev(factorNames) formulaFullModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ ")) factorNames <- head(factorNames,-1) formulaReducedModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ ")) sampleTable formulaFullModel formulaReducedModel primaryFactor countFiles opt$report file.path(opt$report,"DEXSeq_analysis.RData") getwd() dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable, design= formulaFullModel, flattenedfile=opt$gtf) colData(dxd) dxd <- estimateSizeFactors(dxd) sizeFactors(dxd) BPPARAM=MulticoreParam(workers=opt$threads) dxd <- estimateDispersions(dxd, formula=formulaFullModel, BPPARAM=BPPARAM) dxd <- testForDEU(dxd, reducedModel=formulaReducedModel, fullModel=formulaFullModel, BPPARAM=BPPARAM) dxd <- estimateExonFoldChanges(dxd, fitExpToVar=primaryFactor, BPPARAM=BPPARAM) res <- DEXSeqResults(dxd) table(res$padj <= opt$fdr) resSorted <- res[order(res$padj),] head(resSorted) write.csv(as.data.frame(resSorted), file=opt$outfile) if ( !is.null(opt$report) ) { save(dxd, resSorted, file = file.path(opt$report,"DEXSeq_analysis.RData")) save.image() DEXSeqHTML(res, FDR=opt$fdr, color=c("#C3EEE7","#B7FEA0","#F1E7A1","#CEAEFF","#FF8F43","#EDC3C5","#AAA8AA","#FF0000","#637EE9","#FBFBFB"), BPPARAM=BPPARAM) unlink(file.path(opt$report,"DEXSeq_analysis.RData")) } sessionInfo()