Mercurial > repos > vipints > deseq_hts

comparison deseq-hts_2.0/src/difftest_deseq2.R @ 10:2fe512c7bfdf draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DESeq2 version 1.0.19 added to the repo
author vipints <vipin@cbio.mskcc.org>
date Tue, 08 Oct 2013 08:15:34 -0400
parents
children
comparison
equal deleted inserted replaced
9:e27b4f7811c2 10:2fe512c7bfdf
1 ### load DESeq package
2 suppressMessages(require("DESeq2"))
3
4 ### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
5 args <- commandArgs()
6 FITTYP<-args[4]
7 INFILE<-args[5]
8 OUTFILE<-args[6]
9
10 INFILE_COUNTS=c(paste(INFILE, "_COUNTS.tab", sep=""))
11 INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
12
13 ### read count data from file
14 countsTable <- read.delim( INFILE_COUNTS, header=TRUE, stringsAsFactors=TRUE)
15 condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE)
16
17 tagnames <- countsTable[-(1:2), 1]
18 #print(tagnames)
19
20 ## use gene IDs as row names
21 rownames( countsTable ) <- countsTable$gene
22 countsTable <- countsTable[ , -1 ]
23 #print(countsTable)
24
25 conditions<-factor( condsTable[ , 2] )
26 #print(conditions)
27
28 ## unique condition to define the pair of tests
29 uniq_conds <- unique(conditions)
30 #print(uniq_conds)
31
32 ## all possible pairs of conditions
33 pw_tests <- list()
34 for(i in 1:(length(uniq_conds)-1)) {
35 for(j in (i+1):length(uniq_conds)) {
36 pw_tests[[length(pw_tests)+1]] <- c(uniq_conds[i], uniq_conds[j])
37 }
38 }
39 #print(pw_tests)
40
41 tab <- NULL
42 ## testing all possible pairs of conditions
43 for(i in 1:length(pw_tests)) {
44 ## header name
45 test_pair_name <- c(paste(pw_tests[[i]][1], "__vs__", pw_tests[[i]][2], sep=""))
46 #print(test_pair_name)
47 ## colnames respective to the test pair
48 sub.data <- subset(condsTable, (conditions %in% c(pw_tests[[i]][1],pw_tests[[i]][2])))
49 #print(sub.data)
50 #print(sub.data[[1]]) # sample file name
51 #print(sub.data[[2]]) # condition
52 #print(sub.data[[3]]) # replicates
53 colData <- data.frame(row.names=sub.data[[1]], condition=sub.data[[2]], libType=sub.data[[3]])
54 #print(colData)
55 #print(countsTable[(sub.data[[1]])])
56 dds <- DESeqDataSetFromMatrix(countData=countsTable[(sub.data[[1]])], colData=colData, design=~condition)
57 colData(dds)$condition <- factor(colData(dds)$condition, levels=unique(sub.data[[2]]))
58 dds <- DESeq(dds, fitType=FITTYP)
59 ## concatenate the results
60 tested_pairs <- results(dds)
61 #print(typeof(tested_pairs))
62
63 colnames(tested_pairs) <- paste(test_pair_name, colnames(tested_pairs), sep=":")
64 #print(colnames(tested_pairs))
65 #print(tested_pairs)
66
67 tab_tmp <- tested_pairs[tagnames,]
68 if(is.null(tab)) {
69 tab<- as.data.frame(tab_tmp)
70 }
71 else tab<- cbind(tab, as.data.frame(tab_tmp))
72 }
73 ## TODO cbind creates a X character to the string start place.
74 colnames(tab) <- gsub("^X", "", colnames(tab))
75 ## adding gene names to the row
76 tab <- cbind(Feature=row.names(tab_tmp), tab)
77 ## priting the result
78 write.table(tab, OUTFILE, quote=F, sep="\t", row.names=F)