Mercurial > repos > yufei-luo > s_mart

comparison SMART/DiffExpAnal/testR.R @ 18:94ab73e8a190

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Mon, 29 Apr 2013 03:20:15 -0400
parents
children
comparison
equal deleted inserted replaced
17:b0e8584489e6 18:94ab73e8a190
1 #!/usr/bin
2
3 library(DESeq)
4 library(hexbin)
5 library(latticeExtra)
6 library(gplots)
7 library(geneplotter)
8 library(Biobase)
9
10 ##In a file called test_args.R
11 args <- commandArgs()
12
13
14 fileName <- args[4]
15 colNames <- as.integer(unlist(strsplit(args[5], ",")))
16 colCond1 <- as.integer(unlist(strsplit(args[6], ",")))
17 colCond2 <- as.integer(unlist(strsplit(args[7], ",")))
18 OUTPUTCSV <- args[8]
19 OUTPUTPNG <- args[9]
20
21 if(colNames[1]!=0){
22 countsTable <- read.delim(fileName, row.names=1)
23 conditions <- c((colNames[length(colNames)]+1):ncol(countsTable))
24 } else if(colNames[1]==0){
25 countsTable <- read.delim(fileName)
26 conditions <- c(1:ncol(countsTable))
27 rownames(countsTable) <- paste( "Gene", 1:nrow(countsTable), sep="_" )}
28
29 for(i in colCond1){conditions[i] = "A"}
30 for(i in colCond2){conditions[i] = "B"}
31 conditions
32 #analysis with DESeq
33 cds <- newCountDataSet( countsTable, conditions )
34 cds <- estimateSizeFactors( cds )
35 cds <- estimateVarianceFunctions( cds )
36 result <- nbinomTest( cds, "A", "B" )
37 #stock the result dans un .tsv as output file
38 write.table(result, OUTPUTCSV, sep = " ", quote = FALSE, col.names = NA)
39
40 #figures for DE analysis
41 #pdf( OUTPUTPNG, width=4, height=4 )
42 png( filename=OUTPUTPNG, width=700, height=700 )
43 #png format is not as clear as pdf format!!!!!!!!!!!!!!!!!!!!!!!!!
44 print(xyplot(
45 log2FoldChange ~ I(baseMean),
46 result,
47 pch=16, cex=.3,
48 col=ifelse(result$padj < .1, "#FF000040","#00000040" ),
49 panel = function( x, y, col, ...) {
50 above <- (y > 5.8)
51 below <- (y < -5.8)
52 inside <- !( above | below )
53 panel.xyplot( x=x[inside], y=y[inside], col=col[inside], ...)
54 panel.arrows( x[above], 5.8, x[above], 5.95, col=col[above],length=".1", unit="native" )
55 panel.arrows( x[below], -5.8, x[below], -5.95, col=col[below],length=".1", unit="native" ) },
56 axis = function(side, ...) {
57 if( side=="left") {
58 panel.axis( side, outside=TRUE, at=seq(-14,14,by=1), labels=FALSE )
59 panel.axis( side, outside=TRUE, at=seq(-10,10,by=5), labels=TRUE )
60 }
61 if( side=="bottom") {
62 panel.axis( side, outside=TRUE, at=seq(-2,10,by=1), rot=0,
63 labels = do.call( expression,
64 lapply( seq(-2,10,by=1), function(a)
65 substitute( 10^b, list(b=a) ) ) ) )
66 } },
67 xlab = "mean", ylab = "log2 fold change",
68 scales = list(x = list( log=TRUE ),y = list( log=FALSE, limits=c( -6, 6 ) ) ) ))
69 dev.off()
70
71 #The volcano plot
72 #pdf( "vulcano_fly.pdf", width=4, height=4 )
73 #print(xyplot( -log10( pval ) ~ log2FoldChange,
74 # result,
75 # pch=20, cex=.2,
76 # col=ifelse( result$padj<.1, "#FF000050", "#00000050" ),
77 # axis = function( side, ... ) {
78 # if( side=="bottom") {
79 # panel.axis( side, outside=TRUE, at=seq(-14,14,by=1), labels=FALSE )
80 # panel.axis( side, outside=TRUE, at=seq(-10,10,by=5), labels=TRUE )
81 # }
82 # if( side=="left") {
83 # panel.axis( side, outside=TRUE, at=seq(0,25,by=1), labels=FALSE )
84 # panel.axis( side, outside=TRUE, at=seq(0,25,by=5),
85 # labels = do.call( expression,
86 # lapply( seq(0,25,by=5), function(a)
87 # substitute( 10^-b, list(b=a) ) ) ) )
88 # } },
89 # xlab = "log2 fold change", ylab = "p value",
90 # scales = list(
91 # x = list( limits=c( -6, 6 ) ),
92 # y = list( limits=c( 0, 25 ) ) ) ))
93 #dev.off()