|
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",
|
|
15
|
21 'reportdir', 'r', 1, "character",
|
|
1
|
22 'factors', 'f', 1, "character",
|
|
|
23 'threads', 'p', 1, "integer",
|
|
15
|
24 'fdr', 'c', 1, "double"
|
|
1
|
25 ), byrow=TRUE, ncol=4);
|
|
|
26 opt = getopt(spec);
|
|
|
27
|
|
|
28 # if help was asked for print a friendly message
|
|
|
29 # and exit with a non-zero error code
|
|
|
30 if ( !is.null(opt$help) ) {
|
|
|
31 cat(getopt(spec, usage=TRUE));
|
|
|
32 q(status=1);
|
|
|
33 }
|
|
|
34
|
|
|
35 trim <- function (x) gsub("^\\s+|\\s+$", "", x)
|
|
|
36 opt$samples <- trim(opt$samples)
|
|
|
37 opt$factors <- trim(opt$factors)
|
|
|
38
|
|
|
39 parser <- newJSONParser()
|
|
|
40 parser$addData( opt$factors )
|
|
|
41 factorsList <- parser$getObject()
|
|
|
42
|
|
|
43 sampleTable<-data.frame()
|
|
|
44 countFiles<-c()
|
|
|
45 factorNames<-c()
|
|
|
46 primaryFactor<-""
|
|
|
47 for(factor in factorsList){
|
|
|
48 factorName<-factor[[1]]
|
|
|
49 factorNames<-append(factorNames, paste(factorName,"exon",sep=":"))
|
|
|
50 factorValuesMapList<-factor[[2]]
|
|
|
51 c = length(factorValuesMapList)
|
|
|
52 for (factorValuesMap in factorValuesMapList){
|
|
|
53 for(files in factorValuesMap){
|
|
|
54 for(file in files){
|
|
|
55 if(primaryFactor == "") {
|
|
|
56 countFiles<-append(countFiles,file)
|
|
|
57 }
|
|
|
58 sampleTable[basename(file),factorName]<-paste(c,names(factorValuesMap),sep="_")
|
|
|
59 }
|
|
|
60 }
|
|
|
61 c = c-1
|
|
|
62 }
|
|
|
63 if(primaryFactor == ""){
|
|
|
64 primaryFactor <- factorName
|
|
|
65 }
|
|
|
66 }
|
|
|
67
|
|
|
68 factorNames<-append(factorNames,"exon")
|
|
|
69 factorNames<-append(factorNames,"sample")
|
|
|
70 factorNames<-rev(factorNames)
|
|
|
71 formulaFullModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ "))
|
|
|
72 factorNames <- head(factorNames,-1)
|
|
|
73 formulaReducedModel <- as.formula(paste("", paste(factorNames, collapse=" + "), sep=" ~ "))
|
|
|
74
|
|
|
75 sampleTable
|
|
|
76 formulaFullModel
|
|
|
77 formulaReducedModel
|
|
|
78 primaryFactor
|
|
|
79 countFiles
|
|
15
|
80 opt$reportdir
|
|
|
81 opt$threads
|
|
1
|
82 getwd()
|
|
|
83
|
|
|
84 dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable, design= formulaFullModel, flattenedfile=opt$gtf)
|
|
|
85
|
|
|
86 colData(dxd)
|
|
|
87 dxd <- estimateSizeFactors(dxd)
|
|
22
|
88 print("Estimated size factors")
|
|
1
|
89 sizeFactors(dxd)
|
|
|
90 BPPARAM=MulticoreParam(workers=opt$threads)
|
|
|
91 dxd <- estimateDispersions(dxd, formula=formulaFullModel, BPPARAM=BPPARAM)
|
|
22
|
92 print("Estimated dispersions")
|
|
|
93 dxd <- testForDEU(dxd, fullModel=formulaFullModel, BPPARAM=BPPARAM)
|
|
|
94 print("tested for DEU")
|
|
1
|
95 dxd <- estimateExonFoldChanges(dxd, fitExpToVar=primaryFactor, BPPARAM=BPPARAM)
|
|
22
|
96 print("Estimated fold changes")
|
|
1
|
97 res <- DEXSeqResults(dxd)
|
|
22
|
98 print("Results")
|
|
1
|
99 table(res$padj <= opt$fdr)
|
|
|
100 resSorted <- res[order(res$padj),]
|
|
|
101 head(resSorted)
|
|
|
102
|
|
15
|
103 write.table(as.data.frame(resSorted), file = opt$outfile, sep="\t", quote = FALSE, col.names = FALSE)
|
|
22
|
104 print("Written Results")
|
|
1
|
105
|
|
15
|
106 if ( !is.null(opt$reportdir) ) {
|
|
|
107 save(dxd, resSorted, file = file.path(opt$reportdir,"DEXSeq_analysis.RData"))
|
|
1
|
108 save.image()
|
|
22
|
109 DEXSeqHTML(res, path=opt$reportdir, FDR=opt$fdr, color=c("#B7FEA0", "#FF8F43", "#637EE9", "#FF0000", "#F1E7A1", "#C3EEE7","#CEAEFF", "#EDC3C5", "#AAA8AA"))
|
|
15
|
110 unlink(file.path(opt$reportdir,"DEXSeq_analysis.RData"))
|
|
1
|
111 }
|
|
|
112 sessionInfo()
|
