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