## Methods implementation ############################################################################### setMethod("computeDEGs", "TranslatomeDataset", function(object, method="limma", significance.threshold= 0.05, FC.threshold= 0, log.transformed = FALSE, mult.cor=TRUE) { # check input parameters for correctness if (method != "none" & (length(object@cond.a)<2 | length(object@cond.b)<2 | length(object@cond.c)<2 | length(object@cond.d)<2)) stop('You need at least two replicates for each condition to calculate DEGs with statistical methods!') if (!(method %in% c("limma", "SAM", "t-test", "TE", "RP", "ANOTA", "DESeq", "edgeR", "none"))) stop('This method is not recognized!') # conditions for the two levels (first is 1,2 and second is 3,4) cond1 <- object@expr.matrix[,object@cond.a] cond2 <- object@expr.matrix[,object@cond.b] cond3 <- object@expr.matrix[,object@cond.c] cond4 <- object@expr.matrix[,object@cond.d] if (!(object@data.type == "ngs") && !log.transformed) { cond1 <- log(cond1, base=2) cond2 <- log(cond2, base=2) cond3 <- log(cond3, base=2) cond4 <- log(cond4, base=2) } cond <- cbind(cond1, cond2) cond.vector <- c(rep(0, ncol(cond1)), rep(1, ncol(cond2))) cond.2 <- cbind(cond3,cond4) cond.2.vector <- c(rep(0, ncol(cond3)), rep(1, ncol(cond4))) # if the chosen method is translational efficiency, put together samples # as first & second level case and first & second level control # meaning: Pol Case vs Sub/Tot Case and Pol Ctrl vs Sub/Tot Ctrl if (method == "TE") { cond <- cbind(cond1, cond3) cond.vector <- c(rep(0, ncol(cond1)), rep(1, ncol(cond3))) cond.2 <- cbind(cond2,cond4) cond.2.vector <- c(rep(0, ncol(cond2)), rep(1, ncol(cond4))) } #Calculation of FC, avg, and sd for the first level FC <- apply(cond, 1, function(x) mean(x[which(cond.vector == 1)], na.rm=TRUE) - mean(x[which(cond.vector == 0)], na.rm=TRUE)) if (object@data.type == "ngs") FC <- apply(cond, 1, function(x) log(mean(x[which(cond.vector == 1)], na.rm=TRUE) / mean(x[which(cond.vector == 0)], na.rm=TRUE), base=2)) avg.trt <- apply(cond, 1, function(x) mean(x[which(cond.vector == 1)], na.rm=TRUE)) sd.trt <- apply(cond, 1, function(x) sd(x[which(cond.vector == 1)], na.rm=TRUE)) avg.ctr <- apply(cond, 1, function(x) mean(x[which(cond.vector == 0)], na.rm=TRUE)) sd.ctr <- apply(cond, 1, function(x) sd(x[which(cond.vector == 0)], na.rm=TRUE)) #Calculation of FC, avg, and sd for the second level FC2 <- apply(cond.2, 1, function(x) mean(x[which(cond.2.vector == 1)],na.rm=TRUE) - mean(x[which(cond.2.vector == 0)],na.rm=TRUE)) if (object@data.type == "ngs") FC <- apply(cond.2, 1, function(x) log(mean(x[which(cond.2.vector == 1)],na.rm=TRUE) / mean(x[which(cond.2.vector == 0)],na.rm=TRUE),base=2)) avg.trt2 <- apply(cond.2, 1, function(x) mean(x[which(cond.2.vector == 1)], na.rm=TRUE)) sd.trt2 <- apply(cond.2, 1, function(x) sd(x[which(cond.2.vector == 1)], na.rm=TRUE)) avg.ctr2 <- apply(cond.2, 1, function(x) mean(x[which(cond.2.vector == 0)], na.rm=TRUE)) sd.ctr2 <- apply(cond.2, 1, function(x) sd(x[which(cond.2.vector == 0)], na.rm=TRUE)) #execution of the selected significance computation method #(sig.matrix is set with the default values for "none" method) sig.matrix <- matrix(0,nrow=nrow(cond),ncol=4) if (method == "RP") sig.matrix <- methodRP(cond, cond.2, cond.vector, cond.2.vector, mult.cor) if (method == "t-test") sig.matrix <- methodTTest(cond, cond.2, cond.vector, cond.2.vector) if (method == "TE") # compute translational efficiency p-values with Limma as if it was # the normal condition, but cond and cond.2 have been built in a # different way (tot/sub + pol ctrl) and (tot/sub + pol case) sig.matrix <- methodLimma(cond, cond.2, cond.vector, cond.2.vector) if (method == "SAM") sig.matrix <- methodSAM(cond, cond.2, cond.vector, cond.2.vector) if (method == "limma") sig.matrix <- methodLimma(cond, cond.2, cond.vector, cond.2.vector) if (method == "ANOTA") sig.matrix <- methodANOTA(cond, cond.2, cond.vector, cond.2.vector) if (method == "DESeq") sig.matrix <- methodDESeq(cond, cond.2, cond.vector, cond.2.vector) if (method == "edgeR") sig.matrix <- methodEdgeR(cond, cond.2, cond.vector, cond.2.vector) #generation of the final matrix col1 <- ifelse(mult.cor, 2, 1) col2 <- ifelse(mult.cor, 4, 3) level1Up <- sig.matrix[,col1] < significance.threshold & FC > FC.threshold level1Down <- sig.matrix[,col1] < significance.threshold & FC < -FC.threshold level2Up <- sig.matrix[,col2] < significance.threshold & FC2 > FC.threshold level2Down <- sig.matrix[,col2] < significance.threshold & FC2 < -FC.threshold UpUp <- level1Up == 1 & level2Up == 1 DownUp <- level1Down == 1 & level2Up == 1 UpDown <- level1Up == 1 & level2Down == 1 DownDown <- level1Down == 1 & level2Down == 1 final.matrix <- cbind(FC, avg.ctr, sd.ctr, avg.trt, sd.trt, sig.matrix[,1:2], FC2, avg.ctr2, sd.ctr2, avg.trt2, sd.trt2, sig.matrix[,3:4], level1Up, level1Down, level2Up, level2Down, UpUp, DownUp, UpDown, DownDown) colnames(final.matrix) <- c( LOG2.FC1=paste("log2 FC (", object@label.level[1], ")", sep=""), AVG11=paste("avg(", object@label.level[1], ", ", object@label.condition[1], ")", sep=""), SD11=paste("sd(", object@label.level[1], ", ", object@label.condition[1], ")", sep=""), AVG12=paste("avg(", object@label.level[1], ", ", object@label.condition[2], ")", sep=""), SD12=paste("sd(", object@label.level[1], ", ", object@label.condition[2], ")", sep=""), PVAL.1=paste(method, ".pval(", object@label.level[1], ")", sep=""), PVAL.MTC.1=paste(method, ".pval.mtc(", object@label.level[1], ")", sep=""), LOG2.FC2=paste("log2 FC (", object@label.level[2], ")", sep=""), AVG21=paste("avg(", object@label.level[2], ", ", object@label.condition[1], ")", sep=""), SD21=paste("sd(", object@label.level[2], ", ", object@label.condition[1], ")", sep=""), AVG22=paste("avg(", object@label.level[2], ", ", object@label.condition[2], ")", sep=""), SD22=paste("sd(", object@label.level[2], ", ", object@label.condition[2], ")", sep=""), PVAL.2=paste(method, ".pval(", object@label.level[2], ")", sep=""), PVAL.MTC.2=paste(method, ".pval.mtc(", object@label.level[2], ")", sep=""), "level1Up", "level1Down", "level2Up", "level2Down", "UpUp", "DownUp", "UpDown", "DownDown") # store and return the obtained DEGs object@DEGs <- new("DEGs", DEGs.table=final.matrix, method=method, significance.threshold=significance.threshold, FC.threshold=FC.threshold, label.level.DEGs=object@label.level, label.condition=object@label.condition)