Mercurial > repos > vipints > deseq_hts

comparison dexseq-hts_1.0/src/run_DEXseq.R @ 11:cec4b4fb30be draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DEXSeq version 1.6 added
author vipints <vipin@cbio.mskcc.org>
date Tue, 08 Oct 2013 08:22:45 -0400
parents
children
comparison
equal deleted inserted replaced
10:2fe512c7bfdf 11:cec4b4fb30be
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)