|
11
|
1 ### load DEXSeq package
|
|
|
2 suppressMessages(require("DEXSeq"))
|
|
|
3
|
|
|
4 ### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
|
|
|
5 args <- commandArgs()
|
|
|
6 INFILE<-args[4]
|
|
|
7 EXTRAPATH<-args[5]
|
|
|
8 annodb<-args[6]
|
|
|
9 OUTFILE<-args[7]
|
|
|
10
|
|
|
11 INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
|
|
|
12
|
|
|
13 ### read count data from file
|
|
|
14 condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE, row.names=1 )
|
|
|
15 #head(condsTable)
|
|
|
16
|
|
|
17 conditions<-factor( condsTable[ , 1] )
|
|
|
18 #print(conditions)
|
|
|
19
|
|
|
20 ## unique condition to define the pair of tests
|
|
|
21 uniq_conds <- unique(conditions)
|
|
|
22 #print(uniq_conds)
|
|
|
23
|
|
|
24 ## all possible pairs of conditions
|
|
|
25 pw_tests <- list()
|
|
|
26 for(i in 1:(length(uniq_conds)-1)) {
|
|
|
27 for(j in (i+1):length(uniq_conds)) {
|
|
|
28 pw_tests[[length(pw_tests)+1]] <- c(uniq_conds[i], uniq_conds[j])
|
|
|
29 }
|
|
|
30 }
|
|
|
31 #print(pw_tests)
|
|
|
32
|
|
|
33 tab <- NULL
|
|
|
34 ## testing all possible pairs of conditions
|
|
|
35 for(i in 1:length(pw_tests)) {
|
|
|
36 ## header name
|
|
|
37 test_pair_name <- c(paste(pw_tests[[i]][1], "__vs__", pw_tests[[i]][2], sep=""))
|
|
|
38 print(test_pair_name)
|
|
|
39
|
|
|
40 sub.data <- subset(condsTable, (conditions %in% c(pw_tests[[i]][1],pw_tests[[i]][2])))
|
|
|
41 sub.data[[1]]<-as.factor(sub.data[[1]])
|
|
|
42
|
|
|
43 ecs = read.HTSeqCounts(countfiles=file.path(EXTRAPATH, row.names(sub.data)), design=sub.data, flattenedfile=annodb)
|
|
|
44
|
|
|
45 ## Normalisation
|
|
|
46 ecs <- estimateSizeFactors(ecs)
|
|
|
47 print(sizeFactors(ecs))
|
|
|
48
|
|
|
49 suppressMessages(require("parallel"))
|
|
|
50 ## Dispersion estimation
|
|
|
51 ecs <- estimateDispersions(ecs, nCores=2)
|
|
|
52 ecs <- fitDispersionFunction(ecs)
|
|
|
53
|
|
|
54 ## Testing for differential exon usage
|
|
|
55 ecs <- testForDEU(ecs, nCores=2)
|
|
|
56
|
|
|
57 ## estimate the fold change
|
|
|
58 ecs <- estimatelog2FoldChanges(ecs, nCores=2)
|
|
|
59
|
|
|
60 ## Result table
|
|
|
61 resultTable <- DEUresultTable(ecs)
|
|
|
62 print(head(resultTable))
|
|
|
63
|
|
|
64 colnames(resultTable) <- paste(test_pair_name, colnames(resultTable), sep=":")
|
|
|
65 #print(colnames(resultTable))
|
|
|
66
|
|
|
67 if(is.null(tab)) {
|
|
|
68 tab<- resultTable
|
|
|
69 }
|
|
|
70 else tab<- cbind(tab, resultTable)
|
|
|
71 }
|
|
|
72 ## printing the result to out file
|
|
|
73 write.table(tab, OUTFILE, quote=F, sep="\t", row.names=F)
|
