Mercurial > repos > vipints > deseq_hts
comparison deseq-hts_2.0/src/difftest_deseq2.R @ 10:2fe512c7bfdf draft
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) |