|
9
|
1 ### load DESeq package
|
|
|
2 suppressMessages(require("DESeq"))
|
|
0
|
3
|
|
|
4 ### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
|
|
|
5 args <- commandArgs()
|
|
|
6 INFILE<-args[4]
|
|
|
7 OUTFILE<-args[5]
|
|
|
8
|
|
|
9 INFILE_COUNTS=c(paste(INFILE, "_COUNTS.tab", sep=""))
|
|
|
10 INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
|
|
|
11
|
|
|
12 ### read count data from file
|
|
|
13 countsTable <- read.delim( INFILE_COUNTS, header=TRUE, stringsAsFactors=TRUE )
|
|
|
14 condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE )
|
|
|
15
|
|
|
16 ### use gene IDs as row names
|
|
|
17 rownames( countsTable ) <- countsTable$gene
|
|
|
18 countsTable <- countsTable[ , -1 ]
|
|
|
19 head( countsTable )
|
|
|
20
|
|
|
21 conds <- factor( condsTable[ , 2] )
|
|
|
22 #head( countsTable )
|
|
|
23
|
|
|
24 cds <- newCountDataSet( round(countsTable), conds )
|
|
|
25 #head( counts(cds) )
|
|
|
26
|
|
|
27 cds <- estimateSizeFactors( cds )
|
|
|
28 #sizeFactors( cds )
|
|
|
29
|
|
|
30 ### estimate variance function, use blind only, if no replicates are provided
|
|
|
31 if (length(levels(conds)) < length(conds))
|
|
|
32 {
|
|
|
33 cds <- estimateDispersions( cds )
|
|
|
34 } else {
|
|
9
|
35 writeLines("\n***You did not enter any replicates! - The results may be less valuable without replicates!***\n")
|
|
0
|
36 cds <- estimateDispersions( cds, method='blind', sharingMode='fit-only')
|
|
|
37 }
|
|
|
38 experiments <- levels(conds)
|
|
|
39
|
|
9
|
40 res_1<-c()
|
|
|
41 res_2<-c()
|
|
|
42 res_3<-c()
|
|
|
43 res_4<-c()
|
|
|
44 res_5<-c()
|
|
|
45 res_6<-c()
|
|
|
46 res_7<-c()
|
|
|
47 res_8<-c()
|
|
0
|
48 table_col_names<-c()
|
|
9
|
49
|
|
0
|
50 for (i in 1:(length(experiments)-1))
|
|
|
51 {
|
|
|
52 for( j in (i+1):(length(experiments)))
|
|
|
53 {
|
|
|
54 print(c(i,j))
|
|
|
55 tempres <- nbinomTest(cds,experiments[i],experiments[j])
|
|
9
|
56 res_1 = cbind(res_1,tempres[,1])
|
|
|
57 res_2 = cbind(res_2,tempres[,2])
|
|
|
58 res_3 = cbind(res_3,tempres[,3])
|
|
|
59 res_4 = cbind(res_4,tempres[,4])
|
|
|
60 res_5 = cbind(res_5,tempres[,5])
|
|
|
61 res_6 = cbind(res_6,tempres[,6])
|
|
|
62 res_7 = cbind(res_7,tempres[,7])
|
|
|
63 res_8 = cbind(res_8,tempres[,8])
|
|
|
64 table_col_names = cbind(table_col_names,paste('cond_', experiments[i], '_vs._cond_', experiments[j], sep='', 'test'))
|
|
0
|
65 }
|
|
|
66 }
|
|
|
67
|
|
9
|
68 DiffTable<-cbind(res_1,res_2,res_3,res_4,res_5,res_6,res_7,res_8)
|
|
|
69 colnames(DiffTable)<-c('feature ID', 'base mean', 'base mean A', 'base mean B', 'fold change', 'log2 fold change','p value', 'adjusted p value')
|
|
0
|
70 write.table(DiffTable, file = OUTFILE, quote = FALSE, sep ="\t", eol ="\n", na = "1.000", dec = ".", row.names = TRUE,col.names =TRUE)
|
