|
18
|
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()
|
