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