annotate deseq-hts_1.0/src/difftest_deseq.R @ 4:fbec5f31a035 draft

Deleted selected files
author vipints
date Wed, 27 Jun 2012 15:37:35 -0400
parents 94a108763d9e
children e27b4f7811c2
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
1 library( DESeq )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
2
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
3 ### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
4 args <- commandArgs()
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
5 INFILE<-args[4]
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
6 OUTFILE<-args[5]
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
7
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
8 INFILE_COUNTS=c(paste(INFILE, "_COUNTS.tab", sep=""))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
9 INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
10
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
11 ### read count data from file
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
12 countsTable <- read.delim( INFILE_COUNTS, header=TRUE, stringsAsFactors=TRUE )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
13 condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
14
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
15 ### use gene IDs as row names
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
16 rownames( countsTable ) <- countsTable$gene
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
17 countsTable <- countsTable[ , -1 ]
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
18 head( countsTable )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
19
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
20 conds <- factor( condsTable[ , 2] )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
21 #head( countsTable )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
22
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
23 cds <- newCountDataSet( round(countsTable), conds )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
24 #head( counts(cds) )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
25
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
26 cds <- estimateSizeFactors( cds )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
27 #sizeFactors( cds )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
28
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
29 ### estimate variance function, use blind only, if no replicates are provided
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
30 if (length(levels(conds)) < length(conds))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
31 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
32 cds <- estimateDispersions( cds )
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
33 } else {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
34 writeLines("\nYou did not enter any replicates! - The results may be less valuable without replicates!\n")
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
35 cds <- estimateDispersions( cds, method='blind', sharingMode='fit-only')
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
36 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
37 experiments <- levels(conds)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
38
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
39 res<-c()
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
40 table_col_names<-c()
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
41 for (i in 1:(length(experiments)-1))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
42 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
43 for( j in (i+1):(length(experiments)))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
44 {
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
45 print(c(i,j))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
46 tempres <- nbinomTest(cds,experiments[i],experiments[j])
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
47 res = cbind(res,tempres[,7])
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
48 #res = cbind(res,tempres[,8])
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
49 table_col_names = cbind(table_col_names,paste('cond_', experiments[i], '_vs._cond_', experiments[j], sep=''))
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
50 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
51 }
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
52
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
53 DiffTable<-res
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
54 rownames(DiffTable)<-rownames(countsTable)
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
55 colnames(DiffTable)<-table_col_names
94a108763d9e deseq-hts version 1.0 wraps the DESeq 1.6.0
vipints
parents:
diff changeset
56 write.table(DiffTable, file = OUTFILE, quote = FALSE, sep ="\t", eol ="\n", na = "1.000", dec = ".", row.names = TRUE,col.names =TRUE)