Browse code

add tRanslatome package

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/tRanslatome@77226 bc3139a8-67e5-0310-9ffc-ced21a209358

Martin Morgan authored on 05/06/2013 23:20:36
Showing 60 changed files

1 1
new file mode 100644
... ...
@@ -0,0 +1,17 @@
1
+Package: tRanslatome
2
+Type: Package
3
+Title: Comparison between multiple levels of gene expression.
4
+Version: 0.99.1
5
+Date: 2013-05-06
6
+Author: Toma Tebaldi, Erik Dassi, Galena Kostoska
7
+Maintainer: Toma Tebaldi <tebaldi@science.unitn.it>
8
+Depends: R (>= 2.15.0), methods, limma, sigPathway, samr, anota, DESeq,
9
+        edgeR, RankProd, topGO, org.Hs.eg.db, GOSemSim, Heatplus,
10
+        gplots, plotrix
11
+Description: Detection of differentially expressed genes (DEGs) from the comparison of two biological conditions (treated vs. untreated, diseased vs. normal, mutant vs. wild-type) among different levels of gene expression (transcriptome ,translatome, proteome), using several statistical methods:  Rank Product, t-test, SAM, Limma, ANOTA, DESeq, edgeR. Possibility to plot the results with scatterplots, histograms, MA plots, standard deviation (SD) plots, coefficient of variation (CV) plots. Detection of significantly enriched post-transcriptional regulatory factors (RBPs, miRNAs, etc) and Gene Ontology terms in the lists of DEGs previously identified for the two expression levels. Comparison of GO terms enriched only in one of the levels or in both. Calculation of the semantic similarity score between the lists of enriched GO terms coming from the two expression levels. Visual examination and comparison of the enriched terms with heatmaps, radar plots and barplots. 
12
+License: LGPL
13
+LazyLoad: yes
14
+biocViews: CellBiology, GeneRegulation, Regulation, GeneExpression,
15
+        DifferentialExpression, Microarray, HighThroughputSequencing,
16
+        QualityControl, GO, MultipleComparisons, Bioinformatics
17
+Packaged: 2013-06-05 06:48:29 UTC; wolf
0 18
new file mode 100644
... ...
@@ -0,0 +1,5 @@
1
+export(newTranslatomeDataset, computeDEGs, getExprMatrix, getConditionA, getConditionB, getConditionC, getConditionD, getDataType, getConditionLabels, getLevelLabels, getDEGs)
2
+export(Histogram, Scatterplot, MAplot, SDplot, CVplot, GOEnrichment, RegulatoryEnrichment, DEGs.table, getDEGsMethod, significance.threshold, FC.threshold, label.level.DEGs, label.condition)
3
+export(Radar, Heatmap, enriched.table, label.level.enriched)
4
+export(GOComparison)
5
+export(IdentityPlot, SimilarityPlot, identity.matrix, similarity.matrix, average.similarity.scores)
0 6
\ No newline at end of file
1 7
new file mode 100644
... ...
@@ -0,0 +1,766 @@
1
+## Class declaration
2
+###############################################################################
3
+
4
+setClass("DEGs",
5
+	representation(
6
+		DEGs.table="matrix", method="character", 
7
+		significance.threshold="numeric", FC.threshold="numeric",
8
+		label.level.DEGs="character", label.condition="character"
9
+	)
10
+)
11
+
12
+
13
+## Generics declaration (getters, setters and methods)
14
+###############################################################################
15
+
16
+setGeneric("DEGs.table", signature="object", 
17
+					function(object) standardGeneric("DEGs.table"))
18
+					
19
+setGeneric("getDEGsMethod", signature="object", 
20
+					function(object) standardGeneric("getDEGsMethod"))
21
+
22
+setGeneric("significance.threshold", signature="object", 
23
+					function(object) standardGeneric("significance.threshold"))
24
+					
25
+setGeneric("FC.threshold", signature="object", 
26
+					function(object) standardGeneric("FC.threshold"))
27
+					
28
+setGeneric("label.level.DEGs", signature="object", 
29
+					function(object) standardGeneric("label.level.DEGs"))
30
+					
31
+setGeneric("label.condition", signature="object", 
32
+					function(object) standardGeneric("label.condition"))
33
+
34
+setGeneric("MAplot", signature="object", 
35
+					function(object, outputformat="on screen",track="") 
36
+					standardGeneric("MAplot"))
37
+					
38
+setGeneric("CVplot", signature="object", 
39
+					function(object, outputformat="on screen",track="") 
40
+					standardGeneric("CVplot"))
41
+					
42
+setGeneric("SDplot", signature="object", 
43
+					function(object, outputformat="on screen",track="") 
44
+					standardGeneric("SDplot"))
45
+					
46
+setGeneric("Histogram", signature="object", 
47
+					function(object, plottype="summary", outputformat="on screen") 
48
+					standardGeneric("Histogram"))
49
+					
50
+setGeneric("Scatterplot", signature="object",
51
+					function(object, outputformat="on screen",track="") 
52
+					standardGeneric("Scatterplot"))
53
+					
54
+setGeneric("GOEnrichment", signature="object", 
55
+					function(object, ontology="all", 
56
+									classOfDEGs="both", test.method="classic",
57
+									test.threshold = 0.05, mult.cor=TRUE) 
58
+					standardGeneric("GOEnrichment"))
59
+		
60
+setGeneric("RegulatoryEnrichment", signature="object", 
61
+					function(object, classOfDEGs="both", 
62
+									 significance.threshold = 0.05, mult.cor=TRUE)
63
+					standardGeneric("RegulatoryEnrichment"))
64
+
65
+
66
+## Getters implementation
67
+###############################################################################
68
+
69
+setMethod("DEGs.table", "DEGs",
70
+	function(object) {
71
+		object@DEGs.table
72
+	}
73
+)
74
+
75
+setMethod("getDEGsMethod", "DEGs",
76
+	function(object) {
77
+		object@method
78
+	}
79
+)
80
+
81
+setMethod("significance.threshold", "DEGs",
82
+	function(object) {
83
+		object@significance.threshold
84
+	}
85
+)
86
+
87
+setMethod("FC.threshold", "DEGs",
88
+	function(object) {
89
+		object@FC.threshold
90
+	}
91
+)
92
+
93
+setMethod("label.level.DEGs", "DEGs",
94
+	function(object) {
95
+		object@label.level.DEGs
96
+	}
97
+)
98
+
99
+setMethod("label.condition", "DEGs",
100
+	function(object) {
101
+		object@label.condition
102
+	}
103
+)
104
+
105
+
106
+## Methods implementation
107
+###############################################################################
108
+
109
+## Implementation of the show method
110
+setMethod("show", "DEGs",
111
+	function(object) {
112
+		print(head(object@DEGs.table))
113
+		print(object@method)
114
+		print(object@significance.threshold)
115
+		print(object@FC.threshold)
116
+		print(object@label.condition)
117
+		print(object@label.level.DEGs)
118
+	}
119
+)
120
+
121
+
122
+## Implementation of the CVplot method
123
+setMethod("CVplot", "DEGs",
124
+	function(object,outputformat="on screen",track="") {
125
+		resultMatrix <- DEGs.table(object)
126
+		a <- grep(c("FC"), colnames(resultMatrix),value=TRUE)
127
+		c <- grep(c("avg"), colnames(resultMatrix),value=TRUE)
128
+		b <- grep(c("sd"), colnames(resultMatrix),value=TRUE)
129
+
130
+		vector1 <- rowMeans(resultMatrix[,b[1:2]]) / 
131
+												abs(rowMeans(resultMatrix[,c[1:2]]))
132
+		vector2 <- resultMatrix[,a[1]] 
133
+
134
+		vector3 <- rowMeans(resultMatrix[,b[3:4]]) /
135
+												abs(rowMeans(resultMatrix[,c[3:4]]))
136
+		vector4 <- resultMatrix[,a[2]] 
137
+
138
+		xlabel <- "CV"
139
+		ylabel <- "log2 FC"
140
+
141
+		if (outputformat == "pdf") pdf(file="CVplot.pdf")
142
+		if (outputformat == "postscript") postscript(file="CVplot.ps")
143
+		if (outputformat == "jpeg") jpeg(filename="CVplot.jpeg")
144
+
145
+		par(fig=c(0,1,0,1), mar=c(4,4,1,2), mgp=c(2, 0.75, 0))	
146
+		layout(matrix(c(1,2), 2, 1, byrow = TRUE))
147
+		plot(vector1, vector2, xlab=xlabel, ylab=ylabel, 
148
+					col="grey", pch=20, frame.plot=TRUE,cex=0.8, 
149
+					main=substr(a[1], 10, nchar(a[1])-1))
150
+
151
+		list1 <- rownames(resultMatrix)[which(resultMatrix[,"level1Up"]==1 | 
152
+																					resultMatrix[,"level1Down"]==1)]
153
+		variable <- which(rownames(resultMatrix) %in% list1)
154
+		color <- "steelblue3"
155
+		points(vector1[variable],vector2[variable],col=color,pch=20,cex=0.8)
156
+		abline(h=0,lty=2,col="darkgray")
157
+		legend("topright", c("DEGs"),fill=c(color), bty="n", cex=0.8)
158
+
159
+		listm <- rownames(resultMatrix)[rownames(resultMatrix)%in%track]
160
+		if (length(listm)>0) {
161
+			variable <- which(rownames(resultMatrix) %in% track)
162
+			points(vector1[variable],vector2[variable],col="white",pch=4,cex=0.5)
163
+			text(vector1[variable],vector2[variable],labels=listm,pos=3,cex=0.7)
164
+		}
165
+
166
+		plot(vector3, vector4, xlab=xlabel, ylab=ylabel, 
167
+					col="grey", pch=20, frame.plot=TRUE,cex=0.8,
168
+					main=substr(a[2], 10, nchar(a[2])-1))
169
+		abline(h=0,lty=2,col="darkgray")
170
+
171
+		list2 <- rownames(resultMatrix)[which(resultMatrix[,"level2Up"]==1 | 
172
+																					resultMatrix[,"level2Down"]==1)]
173
+		variable <- which(rownames(resultMatrix) %in% list2)
174
+		color <- "gold1"
175
+		points(vector3[variable],vector4[variable],col=color,pch=20,cex=0.8)
176
+		abline(h=0,lty=2,col="darkgray")
177
+		legend("topright", c("DEGs"), fill=c(color), bty="n", cex=0.8)
178
+
179
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
180
+		if (length(listm) > 0) {
181
+			variable <- which(rownames(resultMatrix) %in% track)
182
+			points(vector3[variable],vector4[variable],col="white",pch=4,cex=0.5)
183
+			text(vector3[variable],vector4[variable],labels=listm,pos=3,cex=0.7)
184
+		}
185
+
186
+		# if a file device is open, close it.
187
+		if (!(outputformat == "on screen")) dev.off()
188
+	}
189
+)
190
+
191
+
192
+## Implementation of the MAplot method
193
+setMethod("MAplot", "DEGs",
194
+	function(object,outputformat="on screen",track="") {
195
+		resultMatrix <- DEGs.table(object)
196
+		a <- grep(c("FC"),colnames(resultMatrix),value=TRUE)
197
+		b <- grep(c("avg"),colnames(resultMatrix),value=TRUE)
198
+
199
+		vector1 <- rowMeans(resultMatrix[,b[1:2]])
200
+		vector2 <- resultMatrix[,a[1]] 
201
+
202
+		vector3 <- rowMeans(resultMatrix[,b[3:4]])
203
+		vector4 <- resultMatrix[,a[2]] 
204
+
205
+		xlabel="log2 avg signal (A)"
206
+		ylabel="log2 FC (M)"
207
+
208
+		if (outputformat == "pdf") pdf(file="MAplot.DEGs.pdf")
209
+		if (outputformat == "postscript") postscript(file="MAplot.DEGs.ps")
210
+		if (outputformat == "jpeg") jpeg(filename="MAplot.DEGs.jpeg")
211
+
212
+		par(fig=c(0,1,0,1), mar=c(4,4,1,2), mgp=c(2, 0.75, 0))	
213
+		layout(matrix(c(1,2), 2, 1, byrow = TRUE))
214
+		plot(vector1, vector2, xlab=xlabel, ylab=ylabel, 
215
+					col="grey", pch=20, frame.plot=TRUE, cex=0.8, 
216
+					main=substr(a[1], 10, nchar(a[1])-1))
217
+
218
+		list1 <- rownames(resultMatrix)[which(resultMatrix[,"level1Up"]==1 | 
219
+																					resultMatrix[,"level1Down"]==1)]
220
+		variable <- which(rownames(resultMatrix) %in% list1)
221
+		color <- "steelblue3"
222
+		points(vector1[variable],vector2[variable],col=color,pch=20,cex=0.8)
223
+		abline(h=0,lty=2,col="darkgray")
224
+		legend("topright", c("DEGs"), fill=c(color), bty="n", cex=0.8)
225
+
226
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
227
+		if (length(listm)>0) {
228
+			variable <- which(rownames(resultMatrix) %in% track)
229
+			points(vector1[variable],vector2[variable],col="white",pch=4,cex=0.5)
230
+			text(vector1[variable],vector2[variable],labels=listm,pos=3,cex=0.7)
231
+		}
232
+
233
+		plot(vector3, vector4, xlab=xlabel, ylab=ylabel, 
234
+					col="grey", pch=20, frame.plot=TRUE,cex=0.8,
235
+					main=substr(a[2], 10, nchar(a[2])-1))
236
+		abline(h=0,lty=2,col="darkgray")
237
+
238
+		list2 <- rownames(resultMatrix)[which(resultMatrix[,"level2Up"]==1 | 
239
+																					resultMatrix[,"level2Down"]==1)]
240
+		variable <- which(rownames(resultMatrix) %in% list2)
241
+		color <- "gold1"
242
+		points(vector3[variable],vector4[variable],col=color,pch=20,cex=0.8)
243
+		abline(h=0,lty=2,col="darkgray")
244
+		legend("topright", c("DEGs"), fill=c(color), bty="n", cex=0.8)
245
+
246
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
247
+		if (length(listm)>0) {
248
+			variable <- which(rownames(resultMatrix) %in% track)
249
+			points(vector3[variable], vector4[variable], col="white", pch=4, cex=0.5)
250
+			text(vector3[variable], vector4[variable], labels=listm, pos=3, cex=0.7)
251
+		}
252
+
253
+		# if a file device is open, close it.
254
+		if (!(outputformat == "on screen")) dev.off()
255
+	}
256
+)
257
+
258
+
259
+## Implementation of the Histogram method
260
+setMethod("Histogram", "DEGs",
261
+	function(object, plottype="summary", outputformat="on screen") {
262
+
263
+		resultMatrix <- DEGs.table(object)
264
+		a <- grep(c("FC"), colnames(resultMatrix), value=TRUE)
265
+		list1stlevelUP <- rownames(resultMatrix)[
266
+															which(resultMatrix[,"level1Up"]==1)]
267
+		list1stlevelDOWN <- rownames(resultMatrix)[
268
+															which(resultMatrix[,"level1Down"]==1)]
269
+
270
+		list2ndlevelUP <- rownames(resultMatrix)[which(resultMatrix[,"level2Up"]==1)]
271
+		list2ndlevelDOWN <- rownames(resultMatrix)[which(resultMatrix[,"level2Down"]==1)]
272
+
273
+		listUPUP <- rownames(resultMatrix)[which(resultMatrix[,"UpUp"]==1)]
274
+		listUPDOWN <- rownames(resultMatrix)[which(resultMatrix[,"UpDown"]==1)]
275
+		listDOWNUP <- rownames(resultMatrix)[which(resultMatrix[,"DownUp"]==1)]
276
+		listDOWNDOWN <- rownames(resultMatrix)[which(resultMatrix[,"DownDown"]==1)]
277
+
278
+		if (outputformat == "pdf") pdf(file="Histogram.DEGs.pdf")
279
+		if (outputformat == "postscript") postscript(file="Histogram.DEGs.ps")
280
+		if (outputformat == "jpeg") jpeg(filename="Histogram.DEGs.jpeg")
281
+			
282
+		if (plottype=="summary") {
283
+			inputforsummaryplot <- matrix(c(length(list1stlevelUP),
284
+																			length(list1stlevelDOWN), 
285
+																			length(list2ndlevelUP), 
286
+																			length(list2ndlevelDOWN),
287
+																			0, 0), nrow=2, byrow=FALSE)
288
+
289
+			par(fig=c(0, 1, 0, 1), mar=c(4, 5, 4, 2), mgp=c(2, 0.75, 0))
290
+			barplot(as.table(inputforsummaryplot), xlab="number of DEGs",
291
+							names.arg=c(substr(a[1], 10, nchar(a[1])-1), 
292
+													substr(a[2], 10, nchar(a[2])-1),""), 
293
+							col=c("grey25","grey75"), las=1, 
294
+							border=NA, horiz=TRUE, beside=FALSE) 
295
+			
296
+			legend("topright", c("up", "down"), fill=c("grey25", "grey75"), bty="n")
297
+		}
298
+
299
+		if (plottype == "detailed") {
300
+
301
+			inputfordetailedplot <- matrix(
302
+				c(length(list1stlevelUP) - (length(listUPUP) + length(listUPDOWN)),
303
+				length(list1stlevelDOWN) - (length(listDOWNDOWN) + length(listDOWNUP)),
304
+				length(list2ndlevelUP) - (length(listUPUP) + length(listDOWNUP)),
305
+				length(list2ndlevelDOWN) - (length(listDOWNDOWN) + length(listUPDOWN)),
306
+				length(listUPUP), length(listDOWNDOWN), 
307
+				length(listUPDOWN), length(listDOWNUP)),
308
+				nrow=1, byrow=FALSE)
309
+				
310
+			par(fig=c(0,1,0,1), mar=c(4,7,4,2), mgp=c(2, 0.75, 0))
311
+			barplot(as.table(inputfordetailedplot),
312
+				names.arg= c("up / -","down / -","- / up","- / down", 
313
+									   "up / up","down / down","up / down","down / up"),
314
+				las=1,border=NA,horiz=TRUE,
315
+				col=c("steelblue","steelblue4","gold","orange", 
316
+							"chartreuse","chartreuse4","red","red3"),
317
+				beside=TRUE,xlab="number of DEGs")
318
+				
319
+			mtext(paste(substr(a[1], 10, nchar(a[1])-1),
320
+									substr(a[2], 10, nchar(a[2])-1), sep="/"),
321
+						adj=0, line=0, cex=1)
322
+		}
323
+
324
+		# if a file device is open, close it.
325
+		if (!(outputformat == "on screen")) dev.off()
326
+	}
327
+)
328
+  
329
+  
330
+## Implementation of the Scatterplot method  
331
+setMethod("Scatterplot", "DEGs",
332
+	function(object, outputformat="on screen", track="") {
333
+	
334
+		resultMatrix <- DEGs.table(object)
335
+		a <- grep(c("FC"), colnames(resultMatrix), value=TRUE)
336
+		 
337
+		vector1 <- resultMatrix[,a[1]]
338
+		vector2 <- resultMatrix[,a[2]]
339
+
340
+		originalFCvalX <- 2 ^ vector1
341
+		originalFCvalY <- 2 ^ vector2
342
+
343
+		xlabel=substr(a[1], 6, nchar(a[1]))
344
+		ylabel=substr(a[2], 6, nchar(a[2]))
345
+
346
+		if (outputformat == "pdf") pdf(file="Scatterplot.DEGs.pdf")
347
+		if (outputformat == "postscript") postscript(file="Scatterplot.DEGs.ps")
348
+		if (outputformat == "jpeg") jpeg(filename="Scatterplot.DEGs.jpeg")
349
+
350
+		plot(originalFCvalX, originalFCvalY, xlab=xlabel, ylab=ylabel, 
351
+				 col="grey", pch=20, frame.plot=TRUE, cex=0.8, log="xy")
352
+
353
+		mtext(paste("Spearman all:", 
354
+					round(cor.test(vector1, vector2, method="spearman")$estimate, 2),
355
+					" "),
356
+					side=1, adj=1, line=-2.5, cex=0.8)
357
+
358
+		list0 <- rownames(resultMatrix)[c(which(resultMatrix[,"level1Up"]==1), 
359
+																			which(resultMatrix[,"level1Down"]==1),
360
+																			which(resultMatrix[,"level2Up"]==1), 
361
+																			which(resultMatrix[,"level2Down"]==1))]
362
+		variable <- which(rownames(resultMatrix)%in%list0)
363
+
364
+		if (length(list0)>1)
365
+			mtext(paste("Spearman DEGs:",
366
+						round(cor.test(vector1[list0], 
367
+													 vector2[list0],
368
+													 method="spearman")$estimate, 2),
369
+						" "), 
370
+						side=1, adj=1, line=-1.5, cex=0.8)
371
+
372
+		leg <- NULL
373
+		leg.col <- NULL
374
+		list1 <- rownames(resultMatrix)[which(resultMatrix[,"level1Up"]==1 | 
375
+																					resultMatrix[,"level1Down"]==1)]
376
+		variable <- which(rownames(resultMatrix) %in% list1)
377
+		color <- "steelblue3"
378
+		points(originalFCvalX[variable], originalFCvalY[variable],
379
+					 col=color, pch=20, cex=0.8)
380
+		leg <- c(leg, paste(substr(a[1], 10, nchar(a[1])-1), "only"))
381
+		leg.col <- c(leg.col, color)
382
+
383
+		list2 <- rownames(resultMatrix)[which(resultMatrix[,"level2Up"]==1 | 
384
+																					resultMatrix[,"level2Down"]==1)]
385
+		variable <- which(rownames(resultMatrix) %in% list2)
386
+		color <- "gold1"
387
+		points(originalFCvalX[variable], originalFCvalY[variable], 
388
+					 col=color, pch=20, cex=0.8)
389
+		leg <- c(leg, paste(substr(a[2], 10, nchar(a[2])-1), "only"))
390
+		leg.col <- c(leg.col, color)
391
+
392
+		list3 <- rownames(resultMatrix)[which(resultMatrix[,"UpUp"]==1 |
393
+																					resultMatrix[,"DownDown"]==1)]
394
+		variable <- which(rownames(resultMatrix) %in% list3)
395
+		color <- "chartreuse3"
396
+		points(originalFCvalX[variable], originalFCvalY[variable], 
397
+					 col=color, pch=20, cex=0.8)
398
+		leg <- c(leg, "homodirectional")
399
+		leg.col <- c(leg.col, color)
400
+
401
+		list4 <- rownames(resultMatrix)[which(resultMatrix[,"UpDown"]==1 | 
402
+																					resultMatrix[,"DownUp"]==1)]
403
+		variable <- which(rownames(resultMatrix) %in% list4)
404
+		color <- "red2"
405
+		points(originalFCvalX[variable], originalFCvalY[variable],
406
+					 col=color, pch=20, cex=0.8)
407
+		leg <- c(leg, "opposite change")
408
+		leg.col <- c(leg.col, color)
409
+
410
+		abline(h=1, lty=2, col="darkgray")
411
+		abline(v=1, lty=2, col="darkgray")
412
+		legend("topleft", leg, fill=leg.col, bty="n")
413
+
414
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
415
+		if (length(listm) > 0) {
416
+			variable <- which(rownames(resultMatrix) %in% track)
417
+			points(originalFCvalX[variable], originalFCvalY[variable], 
418
+						 col="white", pch=4, cex=0.5)
419
+			text(originalFCvalX[variable], originalFCvalY[variable], 
420
+					 labels=listm, pos=3, cex=0.7)
421
+		}
422
+
423
+		# if a file device is open, close it.
424
+		if (!(outputformat == "on screen")) dev.off()
425
+	}
426
+)
427
+
428
+
429
+## Implementation of the SDplot method
430
+setMethod("SDplot", "DEGs",
431
+	function(object, outputformat="on screen", track="") {
432
+	
433
+		resultMatrix <- DEGs.table(object)
434
+		a <- grep(c("FC"), colnames(resultMatrix), value=TRUE)
435
+		b <- grep(c("sd"), colnames(resultMatrix), value=TRUE)
436
+
437
+		vector1 <- rowMeans(resultMatrix[,b[1:2]])
438
+		vector2 <- resultMatrix[,a[1]] 
439
+
440
+		vector3 <- rowMeans(resultMatrix[,b[3:4]])
441
+		vector4 <- resultMatrix[,a[2]] 
442
+
443
+		xlabel="SD"
444
+		ylabel="log2 FC"
445
+
446
+		if (outputformat == "pdf") pdf(file="SDplot.DEGs.pdf")
447
+		if (outputformat == "ps") postscript(file="SDplot.DEGs.ps")
448
+		if (outputformat == "jpeg") jpeg(filename="SDplot.DEGs.jpeg")
449
+		
450
+		par(fig=c(0,1,0,1), mar=c(4,4,1,2), mgp=c(2, 0.75, 0))	
451
+		layout(matrix(c(1,2), 2, 1, byrow = TRUE))
452
+		plot(vector1, vector2, xlab=xlabel, ylab=ylabel, 
453
+				 col="grey", pch=20, frame.plot=TRUE,cex=0.8,
454
+				 main=substr(a[1], 10, nchar(a[1])-1))
455
+
456
+		list1 <- rownames(resultMatrix)[which(resultMatrix[,"level1Up"]==1 | 
457
+																					resultMatrix[,"level1Down"]==1)]
458
+		variable <- which(rownames(resultMatrix) %in% list1)
459
+		color <- "steelblue3"
460
+		points(vector1[variable], vector2[variable], 
461
+					 col=color, pch=20, cex=0.8)
462
+		abline(h=0, lty=2, col="darkgray")
463
+		legend("topright", c("DEGs"), fill=c(color), bty="n", cex=0.8)
464
+
465
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
466
+		if (length(listm) > 0) {
467
+			variable <- which(rownames(resultMatrix) %in% track)
468
+			points(vector1[variable], vector2[variable], col="white", pch=4, cex=0.5)
469
+			text(vector1[variable], vector2[variable], labels=listm, pos=3, cex=0.7)
470
+		}
471
+
472
+		plot(vector3, vector4, xlab=xlabel, ylab=ylabel, 
473
+				 col="grey", pch=20, frame.plot=TRUE, cex=0.8,
474
+				 main=substr(a[2], 10, nchar(a[2])-1))
475
+		abline(h=0, lty=2, col="darkgray")
476
+
477
+		list2 <- rownames(resultMatrix)[which(resultMatrix[,"level2Up"]==1 | 
478
+																					resultMatrix[,"level2Down"]==1)]
479
+		variable <- which(rownames(resultMatrix) %in% list2)
480
+		color <- "gold1"
481
+		points(vector3[variable], vector4[variable],
482
+					 col=color, pch=20, cex=0.8)
483
+		abline(h=0, lty=2, col="darkgray")
484
+		legend("topright", c("DEGs"), fill=c(color), bty="n", cex=0.8)
485
+
486
+		listm <- rownames(resultMatrix)[rownames(resultMatrix) %in% track]
487
+		if (length(listm) > 0) {
488
+			variable <- which(rownames(resultMatrix) %in% track)
489
+			points(vector3[variable], vector4[variable], 
490
+						 col="white", pch=4, cex=0.5)
491
+			text(vector3[variable], vector4[variable], labels=listm, pos=3, cex=0.7)
492
+		}
493
+
494
+		# if a file device is open, close it.
495
+		if (!(outputformat=="on screen")) dev.off()
496
+	}
497
+)
498
+
499
+
500
+## Implementation of the GOenrichment method
501
+setMethod("GOEnrichment", "DEGs",
502
+	function(object, ontology="all", classOfDEGs="both", test.method="classic", 
503
+					 test.threshold = 0.05, mult.cor=TRUE) {
504
+		
505
+		resultMatrix <- DEGs.table(object)
506
+		a <- grep(c("FC"), colnames(resultMatrix), value=TRUE)
507
+
508
+		myInterestedGenes1stlevel <- rownames(resultMatrix)[
509
+																	which(resultMatrix[,"level1Up"]==1 | 
510
+																				resultMatrix[,"level1Down"]==1)]
511
+		myInterestedGenes2ndlevel <- rownames(resultMatrix)[
512
+																	which(resultMatrix[,"level2Up"]==1 | 
513
+																				resultMatrix[,"level2Down"]==1)]
514
+		
515
+		if (classOfDEGs == "up") {
516
+			myInterestedGenes1stlevel <- rownames(resultMatrix)[
517
+																		which(resultMatrix[,"level1Up"]==1)]
518
+			myInterestedGenes2ndlevel <- rownames(resultMatrix)[
519
+																		which(resultMatrix[,"level2Up"]==1)]
520
+		}
521
+			
522
+		if (classOfDEGs == "down") {
523
+			myInterestedGenes1stlevel <- rownames(resultMatrix)[
524
+																		which(resultMatrix[,"level1Down"]==1)]
525
+			myInterestedGenes2ndlevel <- rownames(resultMatrix)[
526
+																		which(resultMatrix[,"level2Down"]==1)]
527
+		}
528
+
529
+		if (length(myInterestedGenes1stlevel) == 0 && 
530
+				length(myInterestedGenes2ndlevel) == 0) {
531
+			print("There are no genes selected for enrichment analysis.")
532
+			print("Please try to use a less stringent gene filter.")
533
+		}
534
+		else {
535
+			background <- unique((toTable(org.Hs.egSYMBOL)[,"symbol"]))
536
+			label.level.DEGs <- c(substr(a[1], 10, nchar(a[1])-1),
537
+														substr(a[2], 10, nchar(a[2])-1))
538
+
539
+			# Choice of the ontology
540
+			
541
+			if (ontology !=  "all") {
542
+				onto.first <- createspecifictable(background, 
543
+											myInterestedGenes1stlevel, ontology, 
544
+											test.method, test.threshold, mult.cor, 
545
+											label.level.DEGs[1])
546
+											
547
+				onto.second <- createspecifictable(background,
548
+											 myInterestedGenes2ndlevel, ontology, 
549
+											 test.method, test.threshold, mult.cor, 
550
+											 label.level.DEGs[2])
551
+											 
552
+				final.table <- rbind(onto.first, onto.second)
553
+						
554
+				allRes <- new("GOsets", enriched.table = final.table, 
555
+											label.level.enriched = label.level.DEGs)
556
+			} 
557
+			else {		
558
+				BP.first <- createspecifictable(background, myInterestedGenes1stlevel, 
559
+																				"BP", test.method, test.threshold, 
560
+																				mult.cor, label.level.DEGs[1])
561
+				MF.first <- createspecifictable(background, myInterestedGenes1stlevel, 
562
+																				"MF", test.method, test.threshold,
563
+																				mult.cor, label.level.DEGs[1])
564
+				CC.first <- createspecifictable(background, myInterestedGenes1stlevel, 
565
+																				"CC", test.method, test.threshold,
566
+																				mult.cor, label.level.DEGs[1])
567
+				
568
+				BP.second <- createspecifictable(background, myInterestedGenes2ndlevel, 
569
+																				"BP", test.method, test.threshold,
570
+																				mult.cor, label.level.DEGs[2])
571
+				MF.second <- createspecifictable(background, myInterestedGenes2ndlevel, 
572
+																				"MF", test.method, test.threshold,
573
+																				mult.cor, label.level.DEGs[2])
574
+				CC.second <- createspecifictable(background, myInterestedGenes2ndlevel,
575
+																				"CC", test.method, test.threshold,
576
+																				mult.cor, label.level.DEGs[2])
577
+
578
+				final.table <- rbind(BP.first, MF.first, CC.first, 
579
+														 BP.second, MF.second, CC.second)
580
+				allRes <- new("GOsets", enriched.table=final.table, 
581
+											label.level.enriched=label.level.DEGs)
582
+			}
583
+		}
584
+		return(allRes)
585
+	}
586
+)
587
+
588
+
589
+## Implementation of the RegulatoryEnrichment method
590
+setMethod("RegulatoryEnrichment", "DEGs", 
591
+          function(object, classOfDEGs="both", 
592
+									 significance.threshold= 0.05, mult.cor=TRUE) {
593
+  
594
+	resultMatrix <- DEGs.table(object)
595
+	label.fc <- grep(c("FC"), colnames(resultMatrix), value=TRUE)
596
+	label.level.DEGs <- c(substr(label.fc[1], 10, nchar(label.fc[1])-1),
597
+												substr(label.fc[2], 10, nchar(label.fc[2])-1))
598
+
599
+	
600
+	genes.1stlevel <- rownames(resultMatrix)[
601
+															which(resultMatrix[,"level1Up"]==1 | 
602
+																		resultMatrix[,"level1Down"]==1)]
603
+	genes.2ndlevel <- rownames(resultMatrix)[
604
+															which(resultMatrix[,"level2Up"]==1 | 
605
+																		resultMatrix[,"level2Down"]==1)]
606
+		
607
+	if (classOfDEGs == "up") {
608
+		genes.1stlevel <- rownames(resultMatrix)[
609
+															which(resultMatrix[,"level1Up"]==1)]
610
+		genes.2ndlevel <- rownames(resultMatrix)[
611
+															which(resultMatrix[,"level2Up"]==1)]
612
+	}
613
+			
614
+	if (classOfDEGs == "down") {
615
+		genes.1stlevel <- rownames(resultMatrix)[
616
+															which(resultMatrix[,"level1Down"]==1)]
617
+		genes.2ndlevel <- rownames(resultMatrix)[
618
+															which(resultMatrix[,"level2Down"]==1)]
619
+	}
620
+
621
+	if (length(genes.1stlevel) == 0 && length(genes.2ndlevel) == 0) {
622
+		print("There are no genes selected for enrichment analysis.")
623
+		stop("Please try to use a less stringent gene filter.")
624
+	}
625
+	else {
626
+		enriched.1 <- computeGeneListEnrichment(genes.1stlevel,
627
+											label.level.DEGs[1], significance.threshold, mult.cor)
628
+		enriched.2 <- computeGeneListEnrichment(genes.2ndlevel,
629
+											label.level.DEGs[2], significance.threshold, mult.cor)
630
+		
631
+		return(new("EnrichedSets", enriched.table=rbind(enriched.1, enriched.2), 
632
+												 label.level.enriched=label.level.DEGs))	
633
+	}
634
+}
635
+)
636
+
637
+
638
+## Helper functions
639
+###############################################################################
640
+
641
+createspecifictable <- function(background,myInterestedGenes,ontology,
642
+	test.method,test.threshold,mult.cor,level.label) {
643
+	
644
+	xx <- annFUN.org(ontology, feasibleGenes=background, 
645
+									mapping="org.Hs.eg.db", ID="symbol")
646
+	xxxx <- inverseList(xx)
647
+	geneNames <- names(xxxx)
648
+	geneList <- factor(as.integer(geneNames %in% myInterestedGenes))
649
+	names(geneList) <- geneNames
650
+	str(geneList)
651
+	
652
+	GOdata <- new("topGOdata", ontology = ontology, allGenes = geneList,
653
+								annot = annFUN.gene2GO, gene2GO = xxxx, nodeSize= 5)
654
+		
655
+	if (test.method == "classic") 
656
+		resultFisher <- runTest(GOdata, algorithm = "classic", 
657
+														statistic = "fisher")
658
+	if (test.method == "elim") 
659
+		resultFisher <- runTest(GOdata, algorithm = "elim", 
660
+														statistic = "fisher")
661
+	if (test.method == "weight") 
662
+		resultFisher <- runTest(GOdata, algorithm = "weight", 
663
+														statistic = "fisher")
664
+	if (test.method == "weight01") 
665
+		resultFisher <- runTest(GOdata, algorithm = "weight01",
666
+														statistic = "fisher")
667
+	if (test.method == "parentchild") 
668
+		resultFisher <- runTest(GOdata, algorithm = "parentchild", 
669
+														statistic = "fisher")
670
+
671
+	s <- score(resultFisher)
672
+	
673
+	enrich.table <- GenTable(GOdata, Fisher.result = resultFisher, 
674
+													 orderBy = "Fisher.result", 
675
+													 ranksOf = "Fisher.result", topNodes = length(s))
676
+	enrich.table$Fisher.classic.BH <- p.adjust(
677
+																		as.numeric(enrich.table[,"Fisher.result"]),
678
+																		"BH", 
679
+																		n=length(enrich.table[,"Fisher.result"]))
680
+	
681
+	if (!mult.cor) {
682
+		enrich.table <- enrich.table[
683
+											which(
684
+												enrich.table[,"Fisher.result"] <= test.threshold),]
685
+	}
686
+	else {
687
+		enrich.table <- enrich.table[
688
+											which(
689
+												enrich.table[,"Fisher.classic.BH"] <= test.threshold),]
690
+	}
691
+		
692
+	enrich.table <- enrich.table[order(enrich.table[,"Fisher.result"], 
693
+															 na.last = TRUE, decreasing = FALSE),]
694
+	
695
+	OntologyLevel <- matrix(,nrow=dim(enrich.table)[1], ncol=2)
696
+	colnames(OntologyLevel) <- c("ontology", "level")
697
+	OntologyLevel[,1] <- ontology
698
+	OntologyLevel[,2] <- level.label
699
+	final.table <-  cbind (OntologyLevel, enrich.table)
700
+	colnames(final.table) <- c("ontology", "level","GO.ID", "term", "annotated",
701
+														 "significant", "expected",
702
+														 "pv.fisher", "pv.fisher.BH")
703
+	
704
+	return (final.table)
705
+}
706
+
707
+
708
+computeGeneListEnrichment <- 
709
+		function(DEGs.genes, level.label, significance.threshold, mult.cor) {
710
+	
711
+	# explicitly define the two tables (contained in tRanslatomeDataset)
712
+	# to avoid Rcmd check complaining as it does not see these	
713
+	regulatory.elements.regulated <- NULL
714
+	regulatory.elements.counts <- NULL
715
+	
716
+	# we need data contained in our tRanslatome dataset, so load it
717
+	data(tRanslatomeSampleData, envir=environment())
718
+	
719
+		enrichments <- c()
720
+		
721
+		for (i in c(1:nrow(regulatory.elements.regulated))) {
722
+			regulated.list <- strsplit(
723
+								as.character(regulatory.elements.regulated[i,2]), 
724
+								split=",")[[1]]
725
+			
726
+			regel.name <- as.character(regulatory.elements.regulated[i,1])
727
+			regelIdx <- which(regulatory.elements.counts[,1] == regel.name)
728
+			regel.regulated <- as.integer(regulatory.elements.counts[regelIdx,2])
729
+			regel.nonregulated <- as.integer(regulatory.elements.counts[regelIdx,3])
730
+			
731
+			# get genes from input list regulated by this regulator
732
+			regel.input.regulated <- DEGs.genes[
733
+																		which(DEGs.genes %in% regulated.list)] 
734
+			
735
+			# build the contingency table for fisher test
736
+			cont.table <- matrix(nrow=2, ncol=2)
737
+			# number of input genes regulated by this element
738
+			cont.table[1,1] <- length(regel.input.regulated)
739
+			# number of input genes NOT regulated by this element
740
+			cont.table[1,2] <- length(DEGs.genes) - cont.table[1,1]
741
+			# number of non-input genes regulated by this element
742
+			cont.table[2,1] <- as.integer(regel.regulated) - cont.table[1,1]
743
+			# number of non-input genes NOT regulated by this element
744
+			cont.table[2,2] <- (regel.regulated + regel.nonregulated) - 
745
+												 cont.table[2,1]
746
+			
747
+			# compute the fisher test p-value and report the regulator only if
748
+			# the enrichment p-value is below the given significance threshold
749
+			pvalue <- fisher.test(cont.table)$p.value
750
+			if (pvalue <= significance.threshold) 
751
+				enrichments <- rbind(enrichments, c(regel.name, level.label, 
752
+														 cont.table[1,1],
753
+														 paste(regel.input.regulated, collapse=","),
754
+														 pvalue))
755
+		} 
756
+	
757
+	colnames(enrichments) <- c("ID", "level", "number","list","pv.fisher")															 	
758
+	# if requested, adjust the p-value for multiple testing and add
759
+	# the corrected p-value to the enrichments table
760
+	if (mult.cor)
761
+		enrichments <- cbind(enrichments, 
762
+												 "pv.fisher.BH"=p.adjust(enrichments[,5], 
763
+												 method="BH",n=nrow(regulatory.elements.counts)))
764
+												
765
+	return(as.data.frame(enrichments[order(as.numeric(enrichments[,5])),],stringsAsFactors=F))
766
+}
0 767
\ No newline at end of file
1 768
new file mode 100644
... ...
@@ -0,0 +1,235 @@
1
+## Class declaration
2
+###############################################################################
3
+
4
+setClass("EnrichedSets", representation(enriched.table="data.frame", 
5
+																				label.level.enriched="character"))
6
+
7
+
8
+## Generics declaration (getters, setters and methods)
9
+###############################################################################
10
+
11
+setGeneric("enriched.table", signature="object", 
12
+          function(object) standardGeneric("enriched.table"))
13
+          
14
+setGeneric("label.level.enriched", signature="object", 
15
+          function(object) standardGeneric("label.level.enriched"))
16
+
17
+setGeneric("Radar", signature="object", 
18
+          function(object, outputformat="on screen",
19
+                   n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE, ...) 
20
+          standardGeneric("Radar"))
21
+          
22
+setGeneric("Heatmap", signature="object", 
23
+          function(object, outputformat="on screen", 
24
+                   n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE, ...)
25
+          standardGeneric("Heatmap"))
26
+ 
27
+
28
+## Implementation of getters
29
+###############################################################################
30
+
31
+setMethod("enriched.table", "EnrichedSets",
32
+	function(object) {
33
+		object@enriched.table
34
+	}
35
+)
36
+
37
+setMethod("label.level.enriched", "EnrichedSets",
38
+	function(object) {
39
+		object@label.level.enriched
40
+	}
41
+)
42
+
43
+
44
+## Methods implementation
45
+###############################################################################
46
+
47
+setMethod("show", "EnrichedSets",
48
+	function(object) {
49
+		print(apply(object@enriched.table, c(1,2), function(x){
50
+							if(nchar(x) > 30) 
51
+								return(paste(substr(x,1,30),"...",sep=""))
52
+							else
53
+								return(x)
54
+							})
55
+				 )
56
+		print(object@label.level.enriched)
57
+	}
58
+)
59
+
60
+
61
+setMethod("Radar", "EnrichedSets", 
62
+	function(object, outputformat="on screen", 
63
+           n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE) {
64
+	
65
+		pvalList1stlevel <- NULL
66
+		pvalList2ndlevel <- NULL
67
+
68
+		if (outputformat == "pdf") pdf(file="Enrichment.Radar.pdf")
69
+		if (outputformat == "postscript") postscript(file="Enrichment.Radar.ps")
70
+		if (outputformat == "jpeg") 
71
+      jpeg(filename="Enrichment.Radar.jpeg", width = 650, height = 500, 
72
+					 units = "px", pointsize = 12, bg = "transparent")
73
+		if (outputformat == "png") 
74
+			png(filename="Enrichment.Radar.png", width = 650, height = 500, 
75
+					units = "px", pointsize = 12, bg = "transparent")
76
+
77
+		Atable <- enriched.table(object)[,c("level","ID","list","pv.fisher") ]
78
+		if (mult.cor) 
79
+			Atable <- enriched.table(object)[,c("level","ID","list", "pv.fisher.BH") ] 	 
80
+
81
+		indexfirst <- apply(Atable, 1, 
82
+												function(row) all(label.level.enriched(object)[1] %in% row))
83
+		indexsecond <- apply(Atable, 1, 
84
+												function(row) all(label.level.enriched(object)[2] %in% row))
85
+		List1stlevel <- Atable[indexfirst, -c(1)] 
86
+		List2ndlevel <- Atable[indexsecond, -c(1)]
87
+
88
+		uniquelist <- unique(c(List1stlevel[1:n.nodes.1stlevel,"ID"], 
89
+												   List2ndlevel[1:n.nodes.2ndlevel,"ID"]))
90
+		uniquelist <- uniquelist[!is.na(uniquelist)]
91
+
92
+		if (length(List1stlevel[,1]) == 0 | length(List2ndlevel[,1]) == 0) {
93
+			print("One of the lists is empty!")
94
+		}
95
+		else {
96
+			for(i in 1:length(uniquelist)) {
97
+				print(uniquelist[i])
98
+				
99
+				if ((uniquelist[i] %in% List1stlevel[,"ID"]) == TRUE) {
100
+					pvalList1stlevel <- c(pvalList1stlevel, List1stlevel[
101
+																which(List1stlevel[,"ID"] == uniquelist[i]),
102
+																3])
103
+				}
104
+				else {
105
+					pvalList1stlevel <- c(pvalList1stlevel, 1)
106
+				}
107
+				
108
+				if ((uniquelist[i] %in% List2ndlevel[,"ID"]) == TRUE) {
109
+					pvalList2ndlevel  <- c(pvalList2ndlevel, List2ndlevel[
110
+																 which(List2ndlevel[,"ID"]==uniquelist[i]), 
111
+																 3])
112
+				}
113
+				else {
114
+					pvalList2ndlevel <- c(pvalList2ndlevel, 1)
115
+				}
116
+			}
117
+		}
118
+
119
+		matrix <- matrix(c(-log(as.numeric(pvalList1stlevel), base=10), 
120
+											 -log(as.numeric(pvalList2ndlevel), base=10)), 
121
+										 nrow=2, byrow=TRUE)
122
+		print(matrix)
123
+
124
+		if ((sum(matrix == 0)) == (nrow(matrix) * ncol(matrix))) {
125
+			print("There are no significant enriched regulators in both levels.")
126
+			print("Try to set mult.cor to FALSE.")
127
+		}
128
+		else {
129
+			par(fig=c(0, 1, 0, 1), mar=c(8, 8, 8, 8), mgp=c(2, 0.75, 0))
130
+			radial.plot(matrix, labels=uniquelist, rp.type="p", 
131
+									grid.unit="-log10 p-value", line.col=c("steelblue3","gold1"),
132
+									show.grid.labels=3, lwd=3, start=1,
133
+									point.symbols=18, point.col="black")
134
+									
135
+			legend("topleft", c(label.level.enriched(object)[1], label.level.enriched(object)[2]), 
136
+						 fill=c("steelblue3", "gold1"), bty="n")
137
+		}
138
+		
139
+		# if a file device is open, close it.
140
+		if (!(outputformat == "on screen")) dev.off()
141
+	}
142
+)
143
+
144
+
145
+setMethod("Heatmap", "EnrichedSets",
146
+	function(object, outputformat="on screen", 
147
+					 n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE) {
148
+	
149
+		pvalList1stlevel <- NULL
150
+		pvalList2ndlevel <- NULL
151
+
152
+		if (outputformat == "pdf") pdf(file="Enrichment.Heatmap.pdf")
153
+		if (outputformat =="postscript") postscript(file="Enrichment.Heatmap.ps")
154
+		if (outputformat == "jpeg") 
155
+			jpeg(filename="Enrichment.Heatmap.jpeg", width = 650, height = 500, 
156
+					 units = "px", pointsize = 12, bg = "transparent")
157
+		if (outputformat == "png") 
158
+			png(filename="Enrichment.Heatmap.png", width = 650, height = 500, 
159
+					units = "px", pointsize = 12, bg = "transparent")
160
+
161
+		Atable <- enriched.table(object)[,c("level", "ID", "list", "pv.fisher")]
162
+		if (mult.cor) 
163
+			Atable <- enriched.table(object)[,c("level", "ID", "list", "pv.fisher.BH")] 
164
+
165
+		indexfirst <- apply(Atable, 1, 
166
+												function(row) all(label.level.enriched(object)[1] %in% row))
167
+		indexsecond <- apply(Atable, 1, 
168
+												function(row) all(label.level.enriched(object)[2] %in% row))
169
+		List1stlevel <- Atable[indexfirst, -c(1) ] 
170
+		List2ndlevel <- Atable[indexsecond, -c(1)]
171
+
172
+		uniquelist <- unique(c(List1stlevel[1:n.nodes.1stlevel, "ID"], 
173
+													 List2ndlevel[1:n.nodes.2ndlevel, "ID"]))
174
+		uniquelist <- uniquelist[!is.na(uniquelist)]
175
+		
176
+		if (length(List1stlevel[,1]) == 0 | length(List2ndlevel[,1]) == 0) {
177
+			print("One of the lists is empty!")
178
+		}
179
+		else {
180
+			for(i in 1:length(uniquelist)) {
181
+
182
+				print(uniquelist[i])
183
+				
184
+				if ((uniquelist[i] %in% List1stlevel[,"ID"]) == TRUE) {
185
+					pvalList1stlevel  <-  c(pvalList1stlevel, List1stlevel[
186
+																 which(List1stlevel[,"ID"]==uniquelist[i]),
187
+																 3])
188
+				}
189
+				else {
190
+					pvalList1stlevel <- c(pvalList1stlevel, 1)
191
+				}
192
+					
193
+				if ((uniquelist[i] %in% List2ndlevel[,"ID"]) == TRUE) {
194
+					pvalList2ndlevel  <- c(pvalList2ndlevel, List2ndlevel[
195
+																 which(List2ndlevel[,"ID"]==uniquelist[i]),
196
+																 3])
197
+				}
198
+				else {
199
+					pvalList2ndlevel <- c(pvalList2ndlevel, 1)
200
+				}
201
+			}
202
+		}
203
+
204
+		matrix <- matrix(c(-log(as.numeric(pvalList1stlevel), base=10), 
205
+											 -log(as.numeric(pvalList2ndlevel), base=10)), 
206
+										 ncol=2)
207
+		print(matrix)
208
+
209
+		print(dim(matrix))
210
+		print(length(uniquelist))
211
+
212
+		if ((sum(matrix == 0)) == (nrow(matrix) * ncol(matrix))) {
213
+			print("There are no significant enriched regulators in both levels.")
214
+			print("Try to set mult.cor to FALSE.")
215
+		}
216
+		else {
217
+			rownames(matrix) <- uniquelist
218
+			colnames(matrix) <- c(label.level.enriched(object)[1], 
219
+														label.level.enriched(object)[2])
220
+			print(matrix)
221
+			heatmap.2(matrix, col = RGBColVec(50)[c(26:50)], key=TRUE, 
222
+								margins = c(7,20), keysize=1.5, denscol="white", na.rm = TRUE,
223
+							  scale="none", dendrogram="both", trace="none", 
224
+							  Colv=TRUE, Rowv=TRUE, labRow=NULL, labCol=NULL, 
225
+							  cexRow=1, cexCol=1)
226
+		}
227
+		
228
+		# if a file device is open, close it.
229
+		if (!(outputformat == "on screen")) dev.off()
230
+	}
231
+)
232
+
233
+
234
+## Helper functions
235
+###############################################################################
0 236
new file mode 100644
... ...
@@ -0,0 +1,391 @@
1
+## Class declaration
2
+###############################################################################
3
+
4
+setClass("GOsets", contains="EnrichedSets")
5
+
6
+
7
+## Generics declaration (getters, setters and methods)
8
+###############################################################################
9
+
10
+setGeneric("GOComparison", signature="object", 
11
+          function(object) standardGeneric("GOComparison"))
12
+
13
+
14
+## Implementation of getters
15
+###############################################################################
16
+
17
+
18
+## Methods implementation
19
+###############################################################################
20
+
21
+setMethod("Radar", "GOsets", 
22
+	function(object, outputformat="on screen", 
23
+           n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE,
24
+           ontology="MF") {
25
+	
26
+		pvalList1stlevel <- NULL
27
+		pvalList2ndlevel <- NULL
28
+
29
+		if (outputformat == "pdf") pdf(file="Radar.pdf")
30
+		if (outputformat == "postscript") postscript(file="Radar.ps")
31
+		if (outputformat == "jpeg") 
32
+      jpeg(filename="Radar.jpeg", width = 650, height = 500, units = "px", 
33
+					 pointsize = 12, bg = "transparent")
34
+		if (outputformat == "png") 
35
+			png(filename="Radar.png", width = 650, height = 500, units = "px", 
36
+					pointsize = 12, bg = "transparent")
37
+
38
+		if (length(unique(enriched.table(object)[,1])) == 1 ) {
39
+			Atable <- enriched.table(object)[,c("level", "GO.ID", 
40
+																					"term", "pv.fisher.BH")]
41
+			if (!mult.cor) 
42
+				Atable <- enriched.table(object)[,c("level", "GO.ID", 
43
+																					"term", "pv.fisher") ] 	 
44
+		}
45
+		else {
46
+			# take the rows relevant only to the selected ontology
47
+			index <- apply(enriched.table(object), 1, 
48
+										function(row) all(ontology %in% row))
49
+			Atable <- enriched.table(object)[index,
50
+																c("level", "GO.ID", "term", "pv.fisher.BH")]
51
+			if (!mult.cor) 
52
+				Atable <- enriched.table(object)[index,
53
+																c("level", "GO.ID", "term", "pv.fisher")] 
54
+		}
55
+
56
+		indexfirst <- apply(Atable, 1, 
57
+									function(row) all(label.level.enriched(object)[1] %in% row))
58
+		indexsecond <- apply(Atable, 1, 
59
+									function(row) all(label.level.enriched(object)[2] %in% row))
60
+		List1stlevel <- Atable[indexfirst, -c(1) ] 
61
+		List2ndlevel <- Atable[indexsecond, -c(1)]
62
+
63
+		uniquelist <- unique(c(List1stlevel[1:n.nodes.1stlevel,"GO.ID"], 
64
+												   List2ndlevel[1:n.nodes.2ndlevel,"GO.ID"]))
65
+		uniquelistTerms <- unique(c(List1stlevel[1:n.nodes.1stlevel,"term"], 
66
+																List2ndlevel[1:n.nodes.2ndlevel,"term"]))
67
+
68
+		uniquelist <- uniquelist[!is.na(uniquelist)]
69
+		uniquelistTerms <- uniquelistTerms[!is.na(uniquelistTerms)]
70
+
71
+		if (length(List1stlevel[,1]) == 0 | length(List2ndlevel[,1]) == 0) {
72
+			print("One of the lists is empty!")
73
+		}
74
+		else {
75
+			for(i in 1:length(uniquelist)) {
76
+				print(uniquelist[i])
77
+				print(uniquelistTerms[i])
78
+				
79
+				if ((uniquelist[i] %in% List1stlevel[,"GO.ID"]) == TRUE) {
80
+					pvalList1stlevel <- c(pvalList1stlevel, List1stlevel[
81
+																which(List1stlevel[,"GO.ID"] == uniquelist[i]),
82
+																3])
83
+				}
84
+				else {
85
+					pvalList1stlevel <- c(pvalList1stlevel, 1)
86
+				}
87
+				
88
+				if ((uniquelist[i] %in% List2ndlevel[,"GO.ID"]) == TRUE) {
89
+					pvalList2ndlevel  <- c(pvalList2ndlevel, List2ndlevel[
90
+																 which(List2ndlevel[,"GO.ID"]==uniquelist[i]), 
91
+																 3])
92
+				}
93
+				else {
94
+					pvalList2ndlevel <- c(pvalList2ndlevel, 1)
95
+				}
96
+			}
97
+		}
98
+
99
+		matrix <- matrix(c(-log(as.numeric(pvalList1stlevel), base=10), 
100
+											 -log(as.numeric(pvalList2ndlevel), base=10)), 
101
+										 nrow=2, byrow=TRUE)
102
+		print(matrix)
103
+
104
+		if ((sum(matrix == 0)) == (nrow(matrix) * ncol(matrix))) {
105
+			print("There are no significant enriched terms in both levels.")
106
+			print("Try to set mult.cor to FALSE.")
107
+		}
108
+		else {
109
+			par(fig=c(0, 1, 0, 1), mar=c(8, 8, 8, 8), mgp=c(2, 0.75, 0))
110
+			radial.plot(matrix, labels=uniquelistTerms, rp.type="p", 
111
+									grid.unit="-log10 p-value", line.col=c("steelblue3","gold1"),
112
+									show.grid.labels=3, lwd=3, start=1,
113
+									point.symbols=18, point.col="black")
114
+									
115
+			legend("topleft", c(label.level.enriched(object)[1], 
116
+													label.level.enriched(object)[2]), 
117
+						 fill=c("steelblue3", "gold1"), bty="n")
118
+		}
119
+		
120
+		# if a file device is open, close it.
121
+		if (!(outputformat == "on screen")) dev.off()
122
+	}
123
+)
124
+
125
+
126
+setMethod("Heatmap", "GOsets",
127
+	function(object, outputformat="on screen", 
128
+					 n.nodes.1stlevel="5", n.nodes.2ndlevel="5", mult.cor=TRUE, 
129
+					 ontology="MF") {
130
+	
131
+		pvalList1stlevel <- NULL
132
+		pvalList2ndlevel <- NULL
133
+
134
+		if (outputformat == "pdf") pdf(file="GO.Heatmap.pdf")
135
+		if (outputformat =="postscript") postscript(file="GO.Heatmap.ps")
136
+		if (outputformat == "jpeg") 
137
+			jpeg(filename="GO.Heatmap.jpeg", width = 650, height = 500, units = "px",
138
+					 pointsize = 12, bg = "transparent")
139
+		if (outputformat == "png") 
140
+			png(filename="GO.Heatmap.png", width = 650, height = 500, units = "px", 
141
+					pointsize = 12, bg = "transparent")
142
+
143
+		if (length(unique(enriched.table(object)[,1])) == 1 ) {
144
+			Atable <- enriched.table(object)[,c("level", "GO.ID", 
145
+																					"term", "pv.fisher.BH")]
146
+			if (!mult.cor) 
147
+				Atable <- enriched.table(object)[,c("level", "GO.ID", 
148
+																						"term", "pv.fisher") ] 			 
149
+		}
150
+		else {
151
+			# take the rows relevant only to the selected ontology
152
+			index <- apply(enriched.table(object), 1, 
153
+										 function(row) all(ontology %in% row))
154
+			Atable <- enriched.table(object)[index, 
155
+																 c("level","GO.ID","term","pv.fisher.BH")]
156
+			if (!mult.cor) 
157
+				Atable <- enriched.table(object)[index,
158
+																	 c("level","GO.ID","term","pv.fisher")]  
159
+		}
160
+
161
+		indexfirst <- apply(Atable, 1, 
162
+									function(row) all(label.level.enriched(object)[1] %in% row))
163
+		indexsecond <- apply(Atable, 1, 
164
+									function(row) all(label.level.enriched(object)[2] %in% row))
165
+		List1stlevel <- Atable[indexfirst, -c(1) ] 
166
+		List2ndlevel <- Atable[indexsecond, -c(1)]
167
+
168
+		uniquelist <- unique(c(List1stlevel[1:n.nodes.1stlevel, "GO.ID"], 
169
+													 List2ndlevel[1:n.nodes.2ndlevel, "GO.ID"]))
170
+		uniquelistTerms <- unique(c(List1stlevel[1:n.nodes.1stlevel, "term"], 
171
+																List2ndlevel[1:n.nodes.2ndlevel, "term"]))
172
+
173
+		uniquelist <- uniquelist[!is.na(uniquelist)]
174
+		uniquelistTerms <- uniquelistTerms[!is.na(uniquelistTerms)]
175
+
176
+		if (length(List1stlevel[,1]) == 0 | length(List2ndlevel[,1]) == 0) {
177
+			print("One of the lists is empty!")
178
+		}
179
+		else {
180
+			for(i in 1:length(uniquelist)) {
181
+
182
+				print(uniquelist[i])
183
+				print(uniquelistTerms[i])
184
+				
185
+				if ((uniquelist[i] %in% List1stlevel[,"GO.ID"]) == TRUE) {
186
+					pvalList1stlevel  <-  c(pvalList1stlevel, List1stlevel[
187
+																 which(List1stlevel[,"GO.ID"]==uniquelist[i]),
188
+																 3])
189
+				}
190
+				else {
191
+					pvalList1stlevel <- c(pvalList1stlevel, 1)
192
+				}
193
+					
194
+				if ((uniquelist[i] %in% List2ndlevel[,"GO.ID"]) == TRUE) {
195
+					pvalList2ndlevel  <- c(pvalList2ndlevel, List2ndlevel[
196
+																 which(List2ndlevel[,"GO.ID"]==uniquelist[i]),
197
+																 3])
198
+				}
199
+				else {
200
+					pvalList2ndlevel <- c(pvalList2ndlevel, 1)
201
+				}
202
+			}
203
+		}
204
+
205
+		matrix <- matrix(c(-log(as.numeric(pvalList1stlevel), base=10), 
206
+											 -log(as.numeric(pvalList2ndlevel), base=10)), 
207
+										 ncol=2)
208
+		print(matrix)
209
+
210
+		print(dim(matrix))
211
+		print(dim(uniquelistTerms))
212
+
213
+		if ((sum(matrix == 0)) == (nrow(matrix) * ncol(matrix))) {
214
+			print("There are no significant enriched terms in both levels.")
215
+			print("Try to set mult.cor to FALSE.")
216
+		}
217
+		else {
218
+			rownames(matrix) <- uniquelistTerms
219
+			colnames(matrix) <- c(label.level.enriched(object)[1], 
220
+														label.level.enriched(object)[2])
221
+			print(matrix)
222
+			heatmap.2(matrix, col = RGBColVec(50)[c(26:50)], key=TRUE, 
223
+								margins = c(7,20), keysize=1.5, denscol="white", na.rm = TRUE,
224
+							  scale="none", dendrogram="both", trace="none", 
225
+							  Colv=TRUE, Rowv=TRUE, labRow=NULL, labCol=NULL, 
226
+							  cexRow=1, cexCol=1)
227
+		}
228
+		
229
+		# if a file device is open, close it.
230
+		if (!(outputformat == "on screen")) dev.off()
231
+	}
232
+)
233
+
234
+
235
+setMethod("GOComparison", "GOsets",
236
+	function(object) {
237
+		resulttable <- NULL
238
+
239
+		# identity - stays the same for both of the cases
240
+		# first find the duplicates, then when encountering that duplicate 
241
+		# check wether it is already covered, 
242
+		# 	if not, leave that row as it is, put the label "both levels" 
243
+		# and put it in "coveredDuplicates", 
244
+		# 	if it is already covered mark it for deletion
245
+		
246
+		Atable <- enriched.table(object)[,c("ontology", "level", "GO.ID", "term")]
247
+		duplicates <- Atable[duplicated(Atable[,"GO.ID"]), "GO.ID"]
248
+		Atable$level <- as.character(Atable$level)
249
+		coveredDuplicates <- c()
250
+		rowsForDeletion <- c()
251
+
252
+		for (i in 1:nrow(Atable)) {
253
+			if (any(duplicates == Atable[i, "GO.ID"])) {
254
+				if (all(coveredDuplicates != Atable[i, "GO.ID"])) {
255
+					Atable[i,"level"] <- "both levels"
256
+					coveredDuplicates <- c(coveredDuplicates, Atable[i,"GO.ID"])
257
+				}
258
+				else {
259
+					rowsForDeletion <- c(rowsForDeletion, i)
260
+				}
261
+			}
262
+			else {
263
+				Atable[i,"level"] <- paste(Atable[i,"level"],"only")
264
+			}
265
+		}
266
+
267
+		Atable <- Atable[-c(rowsForDeletion),]
268
+
269
+		if (length(unique(enriched.table(object)[,1])) == 1 ) {
270
+			Btable <- enriched.table(object)[,c("level","GO.ID","term")]
271
+			ontology <- as.character(unique(enriched.table(object)[,1]))
272
+
273
+			indexfirst <- apply(Btable, 1, 
274
+									function(row) all(label.level.enriched(object)[1] %in% row))
275
+			indexsecond <- apply(Btable, 1, 
276
+									function(row) all(label.level.enriched(object)[2] %in% row))
277
+			List1stlevel <- Btable[indexfirst,-c(1)] 
278
+			List2ndlevel <- Btable[indexsecond,-c(1)]
279
+
280
+			if (length(List1stlevel[,1]) == 0 | length(List2ndlevel[,1]) == 0) {
281
+				print("One of the lists is empty!")
282
+			}
283
+			else {
284
+				resulttable1to2 <- comparisonBetweenTwoLists(List1stlevel, 
285
+													List2ndlevel, ontology, 
286
+													paste(label.level.enriched(object)[1], " to ", 
287
+																label.level.enriched(object)[2]))
288
+																
289
+				resulttable2to1 <- comparisonBetweenTwoLists(List2ndlevel, 
290
+													List1stlevel, ontology, 
291
+													paste(label.level.enriched(object)[2], " to ",
292
+																label.level.enriched(object)[1]))
293
+				
294
+				#bind the two result tables
295
+				Similarity.Matrix <- rbind(resulttable1to2, resulttable2to1)
296
+
297
+				totalsimilarity <- c(mean(c(as.numeric(resulttable1to2[,7]), 
298
+																		as.numeric(resulttable2to1[,7]))))
299
+				names(totalsimilarity) <- c(
300
+									as.character(unique(enriched.table(object)[,1])))
301
+			}
302
+		}
303
+		else {	
304
+			ontologyCC <- specificOntologyResult(object, "CC")
305
+			ontologyMF <- specificOntologyResult(object, "MF")
306
+			ontologyBP <- specificOntologyResult(object, "BP")
307
+			
308
+			totalsimilarity <- c(mean(c(as.numeric(ontologyCC$result1[,7]), 
309
+																	as.numeric(ontologyCC$result2[,7]))),
310
+                           mean(c(as.numeric(ontologyMF$result1[,7]), 
311
+																	as.numeric(ontologyMF$result2[,7]))),
312
+                           mean(c(as.numeric(ontologyBP$result1[,7]), 
313
+																	as.numeric(ontologyBP$result2[,7]))))
314
+			names(totalsimilarity) <- c("CC","MF","BP")
315
+
316
+			Similarity.Matrix <- rbind(ontologyCC$result1, ontologyCC$result2, 
317
+																 ontologyMF$result1, ontologyMF$result2, 
318
+																 ontologyBP$result1, ontologyBP$result2)
319
+		}
320
+		
321
+		# return the obtained results table
322
+		return(new("GOsims", 
323
+							similarity.matrix = Similarity.Matrix, 
324
+							identity.matrix = as.matrix(Atable), 
325
+							average.similarity.scores = totalsimilarity))
326
+	}
327
+)
328
+
329
+
330
+## Helper functions
331
+###############################################################################
332
+
333
+specificOntologyResult <- function(object, ontology) {
334
+
335
+	index <- apply(enriched.table(object), 1, 
336
+								 function(row) all(ontology %in% row))
337
+	ontologyTable <- enriched.table(object)[index,c("level", "GO.ID", "term") ] 
338
+	
339
+	indexfirst <- apply(ontologyTable, 1,