1
|
1 ## Setup R error handling to go to stderr
|
|
2 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
|
|
3 # we need that to not crash galaxy with an UTF8 error on German LC settings.
|
|
4 Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
|
|
5
|
|
6 library("DEXSeq")
|
|
7 library('getopt')
|
|
8 library('rjson')
|
|
9
|
|
10
|
|
11 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
|
|
12 args <- commandArgs(trailingOnly = TRUE)
|
|
13
|
|
14 #get options, using the spec as defined by the enclosed list.
|
|
15 #we read the options from the default: commandArgs(TRUE).
|
|
16 spec = matrix(c(
|
|
17 'verbose', 'v', 2, "integer",
|
|
18 'help', 'h', 0, "logical",
|
|
19 'gtf', 'a', 1, "character",
|
|
20 'outfile', 'o', 1, "character",
|
|
21 'report', 'r', 1, "character",
|
|
22 'factors', 'f', 1, "character",
|
|
23 'threads', 'p', 1, "integer",
|
|
24 'fdr', 'c', 1, "double",
|
|
25 ), byrow=TRUE, ncol=4);
|
|
26 opt = getopt(spec);
|
|
27
|
|
28 setwd(opt$report)
|
|
29
|
|
30 # if help was asked for print a friendly message
|
|
31 # and exit with a non-zero error code
|
|
32 if ( !is.null(opt$help) ) {
|
|
33 cat(getopt(spec, usage=TRUE));
|
|
34 q(status=1);
|
|
35 }
|
|
36
|
|
37 trim <- function (x) gsub("^\\s+|\\s+$", "", x)
|
|
38 opt$samples <- trim(opt$samples)
|
|
39 opt$factors <- trim(opt$factors)
|
|
40
|
|
41 parser <- newJSONParser()
|
|
42 parser$addData( opt$factors )
|
|
43 factorsList <- parser$getObject()
|
|
44
|
|
45 sampleTable<-data.frame()
|
|
46 countFiles<-c()
|
|
47 factorNames<-c()
|
|
48 primaryFactor<-""
|
|
49 for(factor in factorsList){
|
|
50 factorName<-factor[[1]]
|
|
51 factorNames<-append(factorNames, paste(factorName,"exon",sep=":"))
|
|
52 factorValuesMapList<-factor[[2]]
|
|
53 c = length(factorValuesMapList)
|
|
54 for (factorValuesMap in factorValuesMapList){
|
|
55 for(files in factorValuesMap){
|
|
56 for(file in files){
|
|
57 if(primaryFactor == "") {
|
|
58 countFiles<-append(countFiles,file)
|
|
59 }
|
|
60 sampleTable[basename(file),factorName]<-paste(c,names(factorValuesMap),sep="_")
|
|
61 }
|
|
62 }
|
|
63 c = c-1
|
|
64 }
|
|
65 if(primaryFactor == ""){
|
|
66 primaryFactor <- factorName
|
|
67 }
|
|
68 }
|
|
69
|
|
70 factorNames<-append(factorNames,"exon")
|
|
71 factorNames<-append(factorNames,"sample")
|
|
72 factorNames<-rev(factorNames)
|
|
73 formulaFullModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ "))
|
|
74 factorNames <- head(factorNames,-1)
|
|
75 formulaReducedModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ "))
|
|
76
|
|
77 sampleTable
|
|
78 formulaFullModel
|
|
79 formulaReducedModel
|
|
80 primaryFactor
|
|
81 countFiles
|
|
82 opt$report
|
|
83 file.path(opt$report,"DEXSeq_analysis.RData")
|
|
84 getwd()
|
|
85
|
|
86 dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable, design= formulaFullModel, flattenedfile=opt$gtf)
|
|
87
|
|
88 colData(dxd)
|
|
89 dxd <- estimateSizeFactors(dxd)
|
|
90 sizeFactors(dxd)
|
|
91 BPPARAM=MulticoreParam(workers=opt$threads)
|
|
92 dxd <- estimateDispersions(dxd, formula=formulaFullModel, BPPARAM=BPPARAM)
|
|
93 dxd <- testForDEU(dxd, reducedModel=formulaReducedModel, fullModel=formulaFullModel, BPPARAM=BPPARAM)
|
|
94 dxd <- estimateExonFoldChanges(dxd, fitExpToVar=primaryFactor, BPPARAM=BPPARAM)
|
|
95 res <- DEXSeqResults(dxd)
|
|
96 table(res$padj <= opt$fdr)
|
|
97 resSorted <- res[order(res$padj),]
|
|
98 head(resSorted)
|
|
99
|
|
100 write.csv(as.data.frame(resSorted), file=opt$outfile)
|
|
101
|
|
102 if ( !is.null(opt$report) ) {
|
|
103 save(dxd, resSorted, file = file.path(opt$report,"DEXSeq_analysis.RData"))
|
|
104 save.image()
|
|
105 DEXSeqHTML(res, FDR=opt$fdr, color=c("#C3EEE7","#B7FEA0","#F1E7A1","#CEAEFF","#FF8F43","#EDC3C5","#AAA8AA","#FF0000","#637EE9","#FBFBFB"), BPPARAM=BPPARAM)
|
|
106 unlink(file.path(opt$report,"DEXSeq_analysis.RData"))
|
|
107 }
|
|
108 sessionInfo()
|