## ------------------------------------------------------------------------ rat = read.table("rat_KD.txt", sep = "\t", header = T) dimnames(rat)[[1]] <- rat[,1] rat = rat[,-1] ## ------------------------------------------------------------------------ colnames(rat) t.test(rat[1,1:6], rat[1,7:11]) y <- t.test(rat[2,1:6], rat[2, 7:11]) dim(rat) ttestRat <- function(df, grp1, grp2) { x = df[grp1] y = df[grp2] x = as.numeric(x) y = as.numeric(y) results = t.test(x, y) results$p.value } rawpvalue = apply(rat, 1, ttestRat, grp1 = c(1:6), grp2 = c(7:11)) ## ------------------------------------------------------------------------ hist(rawpvalue) ## ------------------------------------------------------------------------ ##transform our data into log2 base. rat = log2(rat) #calculate the mean of each gene per control group control = apply(rat[,1:6], 1, mean) #calcuate the mean of each gene per test group test = apply(rat[, 7:11], 1, mean) #confirming that we have a vector of numbers class(control) #confirming we have a vector of numbers class(test) #because our data is already log2 transformed, we can take the difference between the means. And this is our log2 Fold Change or log2 Ratio == log2(control / test) foldchange <- control - test ## ------------------------------------------------------------------------ hist(foldchange, xlab = "log2 Fold Change (Control vs Test)") ## ------------------------------------------------------------------------ results = cbind(foldchange, rawpvalue) results = as.data.frame(results) results$probename <- rownames(results) library(ggplot2) volcano = ggplot(data = results, aes(x = foldchange, y = -1*log10(rawpvalue))) volcano + geom_point()