Browse code

change to DESeq2

erik.dassi@gmail.com authored on 25/11/2020 15:48:58
Showing3 changed files

... ...
@@ -1,11 +1,11 @@
1 1
 Package: tRanslatome
2 2
 Type: Package
3 3
 Title: Comparison between multiple levels of gene expression
4
-Version: 1.25.1
5
-Date: 2020-03-03
4
+Version: 1.26.1
5
+Date: 2020-11-26
6 6
 Author: Toma Tebaldi, Erik Dassi, Galena Kostoska
7 7
 Maintainer: Toma Tebaldi <tebaldi@science.unitn.it>, Erik Dassi <erik.dassi@unitn.it>
8
-Depends: R (>= 2.15.0), methods, limma, sigPathway, anota, DESeq,
8
+Depends: R (>= 2.15.0), methods, limma, sigPathway, anota, DESeq2,
9 9
         edgeR, RankProd, topGO, org.Hs.eg.db, GOSemSim, Heatplus,
10 10
         gplots, plotrix, Biobase
11 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, Translational Efficiency, t-test, 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. 
... ...
@@ -1,5 +1,5 @@
1 1
 import(Biobase)
2
-import(DESeq)
2
+import(DESeq2)
3 3
 import(GOSemSim)
4 4
 import(Heatplus)
5 5
 import(RankProd)
... ...
@@ -1,11 +1,11 @@
1 1
 ## Class declaration
2 2
 ###############################################################################
3 3
 
4
-setClass("TranslatomeDataset", 
5
-	representation(expr.matrix="matrix", 
6
-								 cond.a="character", cond.b="character", 
7
-								 cond.c="character", cond.d="character", 
8
-								 label.level="character", label.condition="character", 
4
+setClass("TranslatomeDataset",
5
+	representation(expr.matrix="matrix",
6
+								 cond.a="character", cond.b="character",
7
+								 cond.c="character", cond.d="character",
8
+								 label.level="character", label.condition="character",
9 9
 								 data.type="character", DEGs="DEGs")
10 10
 )
11 11
 
... ...
@@ -15,34 +15,34 @@ setClass("TranslatomeDataset",
15 15
 
16 16
 setGeneric("getExprMatrix",
17 17
 	function(object) standardGeneric("getExprMatrix"))
18
-	
18
+
19 19
 setGeneric("getConditionA",
20 20
 	function(object) standardGeneric("getConditionA"))
21
-	
21
+
22 22
 setGeneric("getConditionB",
23 23
 	function(object) standardGeneric("getConditionB"))
24
-	
24
+
25 25
 setGeneric("getConditionC",
26 26
 	function(object) standardGeneric("getConditionC"))
27
-	
27
+
28 28
 setGeneric("getConditionD",
29 29
 	function(object) standardGeneric("getConditionD"))
30
-	
30
+
31 31
 setGeneric("getDataType",
32 32
 	function(object) standardGeneric("getDataType"))
33
-	
34
-setGeneric("getConditionLabels", 
33
+
34
+setGeneric("getConditionLabels",
35 35
 	function(object) standardGeneric("getConditionLabels"))
36
-	
36
+
37 37
 setGeneric("getLevelLabels",
38 38
 	function(object) standardGeneric("getLevelLabels"))
39
-	
39
+
40 40
 setGeneric("getDEGs",
41 41
 	function(object) standardGeneric("getDEGs"))
42 42
 
43 43
 setGeneric("computeDEGs", signature="object",
44
-	function(object, method="limma", significance.threshold= 0.05, 
45
-					 FC.threshold= 0,	log.transformed = FALSE, mult.cor=TRUE) 
44
+	function(object, method="limma", significance.threshold= 0.05,
45
+					 FC.threshold= 0,	log.transformed = FALSE, mult.cor=TRUE)
46 46
 	standardGeneric("computeDEGs"))
47 47
 
48 48
 
... ...
@@ -89,21 +89,21 @@ setMethod("getDEGs", "TranslatomeDataset",
89 89
 ## Methods implementation
90 90
 ###############################################################################
91 91
 
92
-setMethod("computeDEGs", "TranslatomeDataset", 
93
-	function(object, method="limma", significance.threshold= 0.05, 
92
+setMethod("computeDEGs", "TranslatomeDataset",
93
+	function(object, method="limma", significance.threshold= 0.05,
94 94
 					 FC.threshold= 0,	log.transformed = FALSE, mult.cor=TRUE) {
95
-	
95
+
96 96
 		# check input parameters for correctness
97
-		if (method != "none" & 
98
-				(length(object@cond.a)<2 | length(object@cond.b)<2 | 
99
-				 length(object@cond.c)<2 | length(object@cond.d)<2)) 
100
-			stop('You need at least two replicates for each condition 
97
+		if (method != "none" &
98
+				(length(object@cond.a)<2 | length(object@cond.b)<2 |
99
+				 length(object@cond.c)<2 | length(object@cond.d)<2))
100
+			stop('You need at least two replicates for each condition
101 101
 			to calculate DEGs with statistical methods!')
102
-		
103
-		if (!(method %in% 
104
-					c("limma", "t-test", "TE", "RP", "ANOTA", "DESeq", "edgeR", "none"))) 
102
+
103
+		if (!(method %in%
104
+					c("limma", "t-test", "TE", "RP", "ANOTA", "DESeq", "edgeR", "none")))
105 105
 			stop('This method is not recognized!')
106
-		
106
+
107 107
 		# conditions for the two levels (first is 1,2 and second is 3,4)
108 108
 		cond1 <- object@expr.matrix[,object@cond.a]
109 109
 		cond2 <- object@expr.matrix[,object@cond.b]
... ...
@@ -112,17 +112,17 @@ setMethod("computeDEGs", "TranslatomeDataset",
112 112
 
113 113
 		if (!(object@data.type == "ngs") && !log.transformed) {
114 114
 			cond1 <- log(cond1, base=2)
115
-			cond2 <- log(cond2, base=2)  
115
+			cond2 <- log(cond2, base=2)
116 116
 			cond3 <- log(cond3, base=2)
117 117
 			cond4 <- log(cond4, base=2)
118
-		}	
119
-	
118
+		}
119
+
120 120
 		cond <- cbind(cond1, cond2)
121 121
 		cond.vector <- c(rep(0, ncol(cond1)), rep(1, ncol(cond2)))
122 122
 
123 123
 		cond.2 <- cbind(cond3,cond4)
124 124
 		cond.2.vector <- c(rep(0, ncol(cond3)), rep(1, ncol(cond4)))
125
-		
125
+
126 126
 		# if the chosen method is translational efficiency, put together samples
127 127
 		# as first & second level case and first & second level control
128 128
 		# meaning: Pol Case vs Sub/Tot Case and Pol Ctrl vs Sub/Tot Ctrl
... ...
@@ -133,14 +133,14 @@ setMethod("computeDEGs", "TranslatomeDataset",
133 133
 			cond.2 <- cbind(cond2,cond4)
134 134
 			cond.2.vector <- c(rep(0, ncol(cond2)), rep(1, ncol(cond4)))
135 135
 		}
136
-		
136
+
137 137
 		#Calculation of FC, avg, and sd for the first level
138 138
 		FC <- apply(cond, 1,
139
-					function(x) mean(x[which(cond.vector == 1)], na.rm=TRUE) - 
139
+					function(x) mean(x[which(cond.vector == 1)], na.rm=TRUE) -
140 140
 											mean(x[which(cond.vector == 0)], na.rm=TRUE))
141
-		if (object@data.type  ==  "ngs") 
141
+		if (object@data.type  ==  "ngs")
142 142
         FC <- apply(cond, 1,
143
-					function(x) log(mean(x[which(cond.vector == 1)], na.rm=TRUE) / 
143
+					function(x) log(mean(x[which(cond.vector == 1)], na.rm=TRUE) /
144 144
 													mean(x[which(cond.vector == 0)], na.rm=TRUE), base=2))
145 145
 		avg.trt <- apply(cond, 1,
146 146
 								function(x) mean(x[which(cond.vector == 1)], na.rm=TRUE))
... ...
@@ -153,11 +153,11 @@ setMethod("computeDEGs", "TranslatomeDataset",
153 153
 
154 154
 		#Calculation of FC, avg, and sd for the second level
155 155
 		FC2 <- apply(cond.2, 1,
156
-						function(x) mean(x[which(cond.2.vector == 1)],na.rm=TRUE) - 
156
+						function(x) mean(x[which(cond.2.vector == 1)],na.rm=TRUE) -
157 157
 												mean(x[which(cond.2.vector == 0)],na.rm=TRUE))
158
-		if (object@data.type  ==  "ngs") 
158
+		if (object@data.type  ==  "ngs")
159 159
       FC2 <- apply(cond.2, 1,
160
-						function(x) log(mean(x[which(cond.2.vector == 1)],na.rm=TRUE) / 
160
+						function(x) log(mean(x[which(cond.2.vector == 1)],na.rm=TRUE) /
161 161
 														mean(x[which(cond.2.vector == 0)],na.rm=TRUE),base=2))
162 162
     avg.trt2 <- apply(cond.2, 1,
163 163
 									function(x) mean(x[which(cond.2.vector == 1)], na.rm=TRUE))
... ...
@@ -168,30 +168,30 @@ setMethod("computeDEGs", "TranslatomeDataset",
168 168
 		sd.ctr2 <- apply(cond.2, 1,
169 169
 									function(x) sd(x[which(cond.2.vector == 0)], na.rm=TRUE))
170 170
 
171
-		#execution of the selected significance computation method 
171
+		#execution of the selected significance computation method
172 172
 		#(sig.matrix is set with the default values for "none" method)
173 173
 		sig.matrix <- matrix(0,nrow=nrow(cond),ncol=4)
174
-		
175
-		if (method == "RP") 
174
+
175
+		if (method == "RP")
176 176
 			sig.matrix <- methodRP(cond, cond.2, cond.vector, cond.2.vector, mult.cor)
177
-		if (method == "t-test") 
177
+		if (method == "t-test")
178 178
 			sig.matrix <- methodTTest(cond, cond.2, cond.vector, cond.2.vector)
179 179
 		if (method == "TE")
180 180
 			# compute translational efficiency p-values with Limma as if it was
181 181
 			# the normal condition, but cond and cond.2 have been built in a
182 182
 			# different way (tot/sub + pol ctrl) and (tot/sub + pol case)
183 183
 			sig.matrix <- methodLimma(cond, cond.2, cond.vector, cond.2.vector)
184
-		if (method == "limma") 
184
+		if (method == "limma")
185 185
 			sig.matrix <- methodLimma(cond, cond.2, cond.vector, cond.2.vector)
186
-		if (method == "ANOTA") 
186
+		if (method == "ANOTA")
187 187
 			sig.matrix <- methodANOTA(cond, cond.2, cond.vector, cond.2.vector)
188
-		if (method == "DESeq") 
188
+		if (method == "DESeq")
189 189
 			sig.matrix <- methodDESeq(cond, cond.2, cond.vector, cond.2.vector)
190
-		if (method == "edgeR") 
190
+		if (method == "edgeR")
191 191
 			sig.matrix <- methodEdgeR(cond, cond.2, cond.vector, cond.2.vector)
192 192
 
193 193
 		#generation of the final matrix
194
-		
194
+
195 195
 		col1 <- ifelse(mult.cor, 2, 1)
196 196
 		col2 <- ifelse(mult.cor, 4, 3)
197 197
 
... ...
@@ -200,55 +200,55 @@ setMethod("computeDEGs", "TranslatomeDataset",
200 200
 		level2Up <- sig.matrix[,col2] < significance.threshold & FC2 > FC.threshold
201 201
 		level2Down <- sig.matrix[,col2] < significance.threshold & FC2 < -FC.threshold
202 202
 
203
-		UpUp <- level1Up == 1 & level2Up == 1			 
203
+		UpUp <- level1Up == 1 & level2Up == 1
204 204
 		DownUp <- level1Down == 1 & level2Up == 1
205 205
 		UpDown <- level1Up == 1 & level2Down == 1
206 206
 		DownDown <- level1Down == 1 & level2Down == 1
207 207
 
208
-		final.matrix <- cbind(FC, avg.ctr, sd.ctr, 
208
+		final.matrix <- cbind(FC, avg.ctr, sd.ctr,
209 209
 													avg.trt, sd.trt, sig.matrix[,1:2],
210
-													FC2, avg.ctr2, sd.ctr2, 
210
+													FC2, avg.ctr2, sd.ctr2,
211 211
 													avg.trt2, sd.trt2, sig.matrix[,3:4],
212 212
 													level1Up, level1Down, level2Up, level2Down,
213 213
 													UpUp, DownUp, UpDown, DownDown)
214
-													
214
+
215 215
 		colnames(final.matrix) <- c(
216
-			LOG2.FC1=paste("log2 FC (", 
217
-										 object@label.level[1], ")", sep=""), 
218
-			AVG11=paste("avg(", object@label.level[1], ", ", 
219
-									object@label.condition[1], ")", sep=""),  
216
+			LOG2.FC1=paste("log2 FC (",
217
+										 object@label.level[1], ")", sep=""),
218
+			AVG11=paste("avg(", object@label.level[1], ", ",
219
+									object@label.condition[1], ")", sep=""),
220 220
 			SD11=paste("sd(", object@label.level[1], ", ",
221
-								 object@label.condition[1], ")", sep=""), 
221
+								 object@label.condition[1], ")", sep=""),
222 222
 			AVG12=paste("avg(", object@label.level[1], ", ",
223
-									object@label.condition[2], ")", sep=""),  
224
-			SD12=paste("sd(", object@label.level[1], ", ", 
225
-								 object@label.condition[2], ")", sep=""), 
226
-			PVAL.1=paste(method, ".pval(", 
227
-									 object@label.level[1], ")", sep=""),  
228
-			PVAL.MTC.1=paste(method, ".pval.mtc(", 
229
-											 object@label.level[1], ")", sep=""), 
230
-			LOG2.FC2=paste("log2 FC (", 
231
-										 object@label.level[2], ")", sep=""), 
232
-			AVG21=paste("avg(", object@label.level[2], ", ", 
233
-									object@label.condition[1], ")", sep=""),  
234
-			SD21=paste("sd(", object@label.level[2], ", ", 
235
-								 object@label.condition[1], ")", sep=""), 
236
-			AVG22=paste("avg(", object@label.level[2], ", ", 
237
-									object@label.condition[2], ")", sep=""),  
238
-			SD22=paste("sd(", object@label.level[2], ", ", 
239
-								 object@label.condition[2], ")", sep=""), 
240
-			PVAL.2=paste(method, ".pval(", 
241
-									 object@label.level[2], ")", sep=""),  
242
-			PVAL.MTC.2=paste(method, ".pval.mtc(", 
243
-											 object@label.level[2], ")", sep=""), 
244
-			"level1Up",  "level1Down",  "level2Up",  "level2Down",  
223
+									object@label.condition[2], ")", sep=""),
224
+			SD12=paste("sd(", object@label.level[1], ", ",
225
+								 object@label.condition[2], ")", sep=""),
226
+			PVAL.1=paste(method, ".pval(",
227
+									 object@label.level[1], ")", sep=""),
228
+			PVAL.MTC.1=paste(method, ".pval.mtc(",
229
+											 object@label.level[1], ")", sep=""),
230
+			LOG2.FC2=paste("log2 FC (",
231
+										 object@label.level[2], ")", sep=""),
232
+			AVG21=paste("avg(", object@label.level[2], ", ",
233
+									object@label.condition[1], ")", sep=""),
234
+			SD21=paste("sd(", object@label.level[2], ", ",
235
+								 object@label.condition[1], ")", sep=""),
236
+			AVG22=paste("avg(", object@label.level[2], ", ",
237
+									object@label.condition[2], ")", sep=""),
238
+			SD22=paste("sd(", object@label.level[2], ", ",
239
+								 object@label.condition[2], ")", sep=""),
240
+			PVAL.2=paste(method, ".pval(",
241
+									 object@label.level[2], ")", sep=""),
242
+			PVAL.MTC.2=paste(method, ".pval.mtc(",
243
+											 object@label.level[2], ")", sep=""),
244
+			"level1Up",  "level1Down",  "level2Up",  "level2Down",
245 245
 			"UpUp",  "DownUp",  "UpDown",  "DownDown")
246 246
 
247 247
 		# store and return the obtained DEGs
248
-		object@DEGs <- new("DEGs", DEGs.table=final.matrix, method=method, 
249
-											 significance.threshold=significance.threshold, 
250
-											 FC.threshold=FC.threshold, 
251
-											 label.level.DEGs=object@label.level, 
248
+		object@DEGs <- new("DEGs", DEGs.table=final.matrix, method=method,
249
+											 significance.threshold=significance.threshold,
250
+											 FC.threshold=FC.threshold,
251
+											 label.level.DEGs=object@label.level,
252 252
 											 label.condition=object@label.condition)
253 253
 		return(object@DEGs)
254 254
 	}
... ...
@@ -275,78 +275,78 @@ setMethod("show", "TranslatomeDataset",
275 275
 
276 276
 # Constructor method for users
277 277
 newTranslatomeDataset <- function(expr.matrix, cond.a, cond.b, cond.c, cond.d,
278
-																	 data.type="array", 
279
-																	 label.level=c("1st level","2nd level"), 
278
+																	 data.type="array",
279
+																	 label.level=c("1st level","2nd level"),
280 280
 																	 label.condition=c("control","treated")) {
281 281
 
282 282
 	# check input parameters for completeness and correctness
283
-  if (missing(expr.matrix) | missing(cond.a) | missing(cond.b) | 
284
-			missing(cond.c) | missing(cond.d)) 
283
+  if (missing(expr.matrix) | missing(cond.a) | missing(cond.b) |
284
+			missing(cond.c) | missing(cond.d))
285 285
 		stop('Some of the mandatory arguments are missing!')
286
-  
287
-  # if the input dataset is a Biobase ExpressionSet, 
286
+
287
+  # if the input dataset is a Biobase ExpressionSet,
288 288
   # extract the expression matrix from it
289 289
   finalMatrix = expr.matrix
290 290
   if (is(expr.matrix, "ExpressionSet")) finalMatrix = exprs(expr.matrix)
291
-  
292
-	return(new(Class="TranslatomeDataset",expr.matrix = finalMatrix, 
293
-						 cond.a=cond.a, cond.b=cond.b, cond.c=cond.c, cond.d=cond.d, 
294
-             data.type=data.type, label.level=label.level, 
291
+
292
+	return(new(Class="TranslatomeDataset",expr.matrix = finalMatrix,
293
+						 cond.a=cond.a, cond.b=cond.b, cond.c=cond.c, cond.d=cond.d,
294
+             data.type=data.type, label.level=label.level,
295 295
              label.condition=label.condition))
296 296
 }
297 297
 
298 298
 # Implementation of the RP helper function
299 299
 
300 300
 methodRP <- function(cond, cond.2, cond.vector, cond.2.vector, mult.cor) {
301
-	
301
+
302 302
 	num.perm = ifelse(mult.cor, 1000, 10)
303 303
 
304
-	RP <- RP(cond, cond.vector, num.perm=num.perm, 
304
+	RP <- RP(cond, cond.vector, num.perm=num.perm,
305 305
 					 logged=TRUE, gene.names=rownames(cond))
306 306
 	rp.pval.1 <- apply(RP$pval, 1, min)
307 307
 	rp.pfp.1 <- apply(RP$pfp, 1, min)
308
-	
309
-	RP2 <- RP(cond.2, cond.2.vector, num.perm=num.perm, 
308
+
309
+	RP2 <- RP(cond.2, cond.2.vector, num.perm=num.perm,
310 310
 						logged=TRUE, gene.names=rownames(cond.2))
311 311
 	rp.pval.2 <- apply(RP2$pval, 1, min)
312 312
 	rp.pfp.2 <-  apply(RP2$pfp, 1, min)
313
-	
314
-	# build the significance p-values matrix and return it	
313
+
314
+	# build the significance p-values matrix and return it
315 315
 	return(cbind(rp.pval.1, rp.pfp.1, rp.pval.2, rp.pfp.2))
316 316
 }
317 317
 
318 318
 # Implementation of the t-test helper function
319 319
 methodTTest <- function(cond, cond.2, cond.vector, cond.2.vector) {
320
-	
320
+
321 321
 	# t.test.pval for first level
322 322
 	t.test.pval <- calcTStatFast(cond, cond.vector)$pval
323
-	t.test.pval.adj <- p.adjust(t.test.pval, 
323
+	t.test.pval.adj <- p.adjust(t.test.pval,
324 324
 															method="BH", n=length(t.test.pval))
325 325
 
326 326
 	# t.test.pval for second level
327 327
 	t.test.pval2 <- calcTStatFast(cond.2, cond.2.vector)$pval
328
-	t.test.pval.adj2 <- p.adjust(t.test.pval2, 
328
+	t.test.pval.adj2 <- p.adjust(t.test.pval2,
329 329
 															 method="BH", n=length(t.test.pval2))
330 330
 
331
-	# build the significance p-values matrix and return it	
332
-	return(cbind(t.test.pval, t.test.pval.adj, 
333
-							 t.test.pval2, t.test.pval.adj2))	
331
+	# build the significance p-values matrix and return it
332
+	return(cbind(t.test.pval, t.test.pval.adj,
333
+							 t.test.pval2, t.test.pval.adj2))
334 334
 }
335 335
 
336 336
 
337 337
 # Implementation of the limma helper function
338 338
 methodLimma <- function(cond, cond.2, cond.vector, cond.2.vector) {
339
-	
340
-	eset <- data.frame(cond.a=cond[,cond.vector==0], 
341
-										 cond.b=cond[,cond.vector==1], 
342
-										 cond.c=cond.2[,cond.2.vector==0], 
339
+
340
+	eset <- data.frame(cond.a=cond[,cond.vector==0],
341
+										 cond.b=cond[,cond.vector==1],
342
+										 cond.c=cond.2[,cond.2.vector==0],
343 343
 										 cond.d=cond.2[,cond.2.vector==1])
344
-										  	
345
-	Target <- c(rep("cond.a", sum(cond.vector==0)), 
346
-							rep("cond.b", sum(cond.vector==1)), 
347
-							rep("cond.c", sum(cond.2.vector==0)), 
344
+
345
+	Target <- c(rep("cond.a", sum(cond.vector==0)),
346
+							rep("cond.b", sum(cond.vector==1)),
347
+							rep("cond.c", sum(cond.2.vector==0)),
348 348
 							rep("cond.d", sum(cond.2.vector==1)))
349
-							
349
+
350 350
 	colnames(eset) <- Target
351 351
 	FileName <- colnames(eset)
352 352
 	targets <- data.frame(FileName=FileName, Target=Target)
... ...
@@ -360,101 +360,108 @@ methodLimma <- function(cond, cond.2, cond.vector, cond.2.vector) {
360 360
 	cont.firstlevel <- makeContrasts("cond.b - cond.a", levels = design)
361 361
 	fitfirstlevel <- contrasts.fit(fit, cont.firstlevel)
362 362
 	fitfirstlevel <- eBayes(fitfirstlevel)
363
-	BH.firstlevel <- p.adjust(fitfirstlevel$F.p.value, method = "BH") 
363
+	BH.firstlevel <- p.adjust(fitfirstlevel$F.p.value, method = "BH")
364 364
 
365 365
 	cont.secondlevel <- makeContrasts("cond.d - cond.c", levels = design)
366 366
 	fitsecondlevel <- contrasts.fit(fit, cont.secondlevel)
367 367
 	fitsecondlevel <- eBayes(fitsecondlevel)
368
-	BH.secondlevel <- p.adjust(fitsecondlevel$F.p.value, method = "BH") 
369
-		
370
-	# build the significance p-values matrix and return it	
371
-	return(cbind(fitfirstlevel$F.p.value, BH.firstlevel, 
368
+	BH.secondlevel <- p.adjust(fitsecondlevel$F.p.value, method = "BH")
369
+
370
+	# build the significance p-values matrix and return it
371
+	return(cbind(fitfirstlevel$F.p.value, BH.firstlevel,
372 372
 							 fitsecondlevel$F.p.value, BH.secondlevel))
373 373
 }
374 374
 
375 375
 # Implementation of the ANOTA helper function
376 376
 methodANOTA <- function(cond, cond.2, cond.vector, cond.2.vector) {
377 377
 	#library(anota)
378
-	
379
-	if (length(cond.vector) != length(cond.2.vector) | 
378
+
379
+	if (length(cond.vector) != length(cond.2.vector) |
380 380
 			!all(cond.vector == cond.2.vector))
381
-		stop('The ANOTA method cannot be applied if the two levels 
381
+		stop('The ANOTA method cannot be applied if the two levels
382 382
 		have a different experimental design!')
383
-	
383
+
384 384
 	#first level
385
-	anotaQcOut <- anotaPerformQc(dataT = cond, dataP = cond.2, 
385
+	anotaQcOut <- anotaPerformQc(dataT = cond, dataP = cond.2,
386 386
 															 phenoVec = cond.vector)
387
-	anotaSigGeneOut <- anotaGetSigGenes(dataT = cond, dataP = cond.2, 
387
+	anotaSigGeneOut <- anotaGetSigGenes(dataT = cond, dataP = cond.2,
388 388
 															 phenoVec = cond.vector, anotaQcObj = anotaQcOut)
389 389
 
390 390
 	pvalues.1st <- anotaSigGeneOut$apvStats[[1]][,"apvP"]
391
-	pvalues.1st.BH <- p.adjust(pvalues.1st, method = "BH") 
391
+	pvalues.1st.BH <- p.adjust(pvalues.1st, method = "BH")
392 392
 
393 393
 	#second level
394
-	anotaQcOut.2nd <- anotaPerformQc(dataT = cond.2, dataP = cond, 
394
+	anotaQcOut.2nd <- anotaPerformQc(dataT = cond.2, dataP = cond,
395 395
 																	 phenoVec = cond.vector)
396
-	anotaSigGeneOut.2nd <- anotaGetSigGenes(dataT = cond.2, dataP = cond, 
397
-																	 phenoVec = cond.vector, 
396
+	anotaSigGeneOut.2nd <- anotaGetSigGenes(dataT = cond.2, dataP = cond,
397
+																	 phenoVec = cond.vector,
398 398
 																	 anotaQcObj = anotaQcOut.2nd)
399 399
 
400 400
 	pvalues.2nd <- anotaSigGeneOut.2nd$apvStats[[1]][,"apvP"]
401
-	pvalues.2nd.BH <- p.adjust(pvalues.2nd, method = "BH") 
401
+	pvalues.2nd.BH <- p.adjust(pvalues.2nd, method = "BH")
402 402
 
403
-	# build the significance p-values matrix and return it	
403
+	# build the significance p-values matrix and return it
404 404
 	return(cbind(pvalues.2nd, pvalues.2nd.BH, pvalues.1st, pvalues.1st.BH))
405 405
 }
406 406
 
407 407
 # Implementation of the DESeq helper function (for NGS data)
408 408
 methodDESeq <- function(cond, cond.2, cond.vector, cond.2.vector) {
409
-	#require(DESeq)
410
-	
411
-	cond.1.deseq <- newCountDataSet(cond, cond.vector)
409
+	#require(DESeq2)
410
+	metadata.1 <- as.data.frame(cbind(cond=cond.vector))
411
+	metadata.1$cond <- factor(metadata.1$cond)
412
+
413
+	cond.1.deseq <- DESeqDataSetFromMatrix(cond, metadata.1, design=~cond)
412 414
 	cond.1.deseq <- estimateSizeFactors(cond.1.deseq)
413 415
   # if GLM fit fails, switch to local fitting
414 416
 	tmp <- tryCatch(cond.1.deseq <- estimateDispersions(cond.1.deseq),
415 417
 								 error=function(ex){ return(-1) })
416 418
 	if (is.numeric(tmp))
417
-		cond.1.deseq <- estimateDispersions(cond.1.deseq, 
419
+		cond.1.deseq <- estimateDispersions(cond.1.deseq,
418 420
 																			fitType="local",
419
-																			sharingMode="fit-only")	
420
-	res.1 <- nbinomTest(cond.1.deseq, "0", "1")
421
-	
422
-	cond.2.deseq <- newCountDataSet(cond.2, cond.2.vector)
421
+																			sharingMode="fit-only")
422
+	res.1 <- nbinomWaldTest(cond.1.deseq)
423
+	res.1 <- results(res.1)
424
+
425
+	metadata.2 <- as.data.frame(cbind(cond=cond.2.vector))
426
+	metadata.2$cond <- factor(metadata.2$cond)
427
+
428
+	cond.2.deseq <- DESeqDataSetFromMatrix(cond.2, metadata.2, design=~cond)
423 429
 	cond.2.deseq <- estimateSizeFactors(cond.2.deseq)
424 430
   # if GLM fit fails, switch to local fitting
425
-  tmp <- tryCatch(cond.2.deseq <-estimateDispersions(cond.2.deseq), 
431
+  tmp <- tryCatch(cond.2.deseq <-estimateDispersions(cond.2.deseq),
426 432
 								 error=function(ex){ return(-1) })
427 433
 	if (is.numeric(tmp))
428
-		cond.2.deseq <- estimateDispersions(cond.2.deseq, 
434
+		cond.2.deseq <- estimateDispersions(cond.2.deseq,
429 435
 																			fitType="local",
430
-																			sharingMode="fit-only")	
431
-	res.2 <- nbinomTest(cond.2.deseq, "0", "1")
432
-	
433
-	# build the significance p-values matrix and return it	
434
-	return(cbind(res.1$pval, res.1$padj, res.2$pval, res.2$padj))
436
+																			sharingMode="fit-only")
437
+	res.2 <- nbinomWaldTest(cond.2.deseq)
438
+	res.2 <- results(res.2)
439
+
440
+	# build the significance p-values matrix and return it
441
+	return(cbind(res.1$pvalue, res.1$padj, res.2$pvalue, res.2$padj))
435 442
 }
436 443
 
437 444
 # Implementation of the edgeR helper function (for NGS data)
438 445
 methodEdgeR <- function(cond, cond.2, cond.vector, cond.2.vector) {
439 446
 	#require(edgeR)
440
-	
447
+
441 448
 	cond.1.edger <- DGEList(counts=cond, group=cond.vector)
442 449
 	cond.1.edger <- calcNormFactors(cond.1.edger)
443 450
 	cond.1.edger <- estimateCommonDisp(cond.1.edger)
444 451
 	cond.1.edger <- estimateTagwiseDisp(cond.1.edger)
445 452
 	res.1 <- exactTest(cond.1.edger)
446
-	
453
+
447 454
 	cond.2.edger <- DGEList(counts=cond.2, group=cond.2.vector)
448 455
 	cond.2.edger <- calcNormFactors(cond.2.edger)
449 456
 	cond.2.edger <- estimateCommonDisp(cond.2.edger)
450 457
 	cond.2.edger <- estimateTagwiseDisp(cond.2.edger)
451 458
 	res.2 <- exactTest(cond.2.edger)
452
-	
453
-	# build the significance p-values matrix and return it	
454
-	return(cbind(res.1$table$PValue, 
459
+
460
+	# build the significance p-values matrix and return it
461
+	return(cbind(res.1$table$PValue,
455 462
 							 p.adjust(res.1$table$PValue,
456
-												method="BH",n=length(res.1$table$PValue)), 
457
-               res.2$table$PValue, 
463
+												method="BH",n=length(res.1$table$PValue)),
464
+               res.2$table$PValue,
458 465
                p.adjust(res.2$table$PValue,
459 466
 												method="BH",n=length(res.2$table$PValue))))
460 467
 }
461 468
\ No newline at end of file