# --------------------------------
# Protocol for simpleaffy exercise 
# --------------------------------

# --Required 'covdesc.txt' file for this example--
#
#	treatment
# COLD_CONTROL_12H_SHOOT_REP1.cel	1
# COLD_CONTROL_12H_SHOOT_REP2.cel	1
# COLD_12H_SHOOT_REP1.cel	2
# COLD_12H_SHOOT_REP2.cel	2
# -------

# --RMA Commands--
library(simpleaffy);
raw.data <- read.affy("covdesc.txt");
x.rma <- call.exprs(raw.data, "rma");
write.exprs(x.rma, file="rma_affy_all.txt");
results <- pairwise.comparison(x.rma, "treatment", c("1", "2"), raw.data);
write.table(data.frame(means(results), fc(results), tt(results), calls(results)), file="my_comp.txt", sep="\t");
significant <- pairwise.filter(results, fc=log2(2), min.present.no=4, tt= 0.001, present.by.group=FALSE)
sig.rma <- data.frame(means(significant), fc(significant), tt(significant), calls(significant))
write.table(sig.rma, file="sig.rma.txt", sep="\t")
png(file="my_comp_rma.png"); plot(significant,type="scatter"); dev.off()
# -------

# --GCRMA Commands--
library(simpleaffy);
library(gcrma);
raw.data <- read.affy("covdesc.txt");
x.gcrma <- call.exprs(raw.data, "gcrma");
write.exprs(x.gcrma, file="gcrma_affy_all.txt");
results <- pairwise.comparison(x.gcrma, "treatment", c("1", "2"), raw.data);
write.table(data.frame(means(results), fc(results), tt(results), calls(results)), file="my_comp.txt", sep="\t");
significant <- pairwise.filter(results, fc=log2(2), min.present.no=4, tt= 0.001, present.by.group=FALSE)
sig.gcrma <- data.frame(means(significant), fc(significant), tt(significant), calls(significant))
write.table(sig.gcrma, file="sig.gcrma.txt", sep="\t")
png(file="my_comp_gcrma.png"); plot(significant,type="scatter"); dev.off()
# -------

# --MAS5 Commands--
library(simpleaffy);
raw.data <- read.affy("covdesc.txt");
x.mas <- call.exprs(raw.data, "mas5");
write.exprs(x.mas, file="mas_affy_all.txt");
results <- pairwise.comparison(x.mas, "treatment", c("1", "2"), raw.data);
write.table(data.frame(means(results), fc(results), tt(results), calls(results)), file="my_comp.txt", sep="\t");
significant <- pairwise.filter(results, fc=log2(2), min.present.no=4, tt= 0.001, present.by.group=FALSE)
sig.mas5 <- data.frame(means(significant), fc(significant), tt(significant), calls(significant))
write.table(sig.mas5, file="sig.mas5.txt", sep="\t")
png(file="my_comp_mas.png"); plot(significant,type="scatter"); dev.off()
# -------


# --Identifiy overlap between three methods--
overlap2 <-(merge(sig.rma, sig.gcrma, by.x = "row.names", by.y = "row.names", all = FALSE))
overlap3 <-(merge(overlap2, sig.mas5, by.x = "Row.names", by.y = "row.names", all = FALSE))	
# -------