######################### ## Exercise Simpleaffy ## ######################### ## 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.xls") 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.xls", sep="\t", quote=F, col.names = NA) significant <- pairwise.filter(results, fc=log2(2), min.present.no=4, tt= 0.01, present.by.group=FALSE) sig.rma <- data.frame(means(significant), fc(significant), tt(significant), calls(significant)) write.table(sig.rma, file="sig.rma.xls", sep="\t", quote=F, col.names = NA) pdf(file="my_comp_rma.pdf"); 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.xls") 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.xls", sep="\t", quote=F, col.names = NA) significant <- pairwise.filter(results, fc=log2(2), min.present.no=4, tt= 0.01, present.by.group=FALSE) sig.mas5 <- data.frame(means(significant), fc(significant), tt(significant), calls(significant)) write.table(sig.mas5, file="sig.mas5.xls", sep="\t",quote=F, col.names = NA) pdf(file="my_comp_mas.pdf"); plot(significant,type="scatter"); dev.off() ###### ## Identifiy overlap between two methods ## overlap <-(merge(sig.rma, sig.mas5, by.x = "row.names", by.y = "row.names", all = FALSE)) write.table(overlap, file="overlap.xls", sep="\t",quote=F, col.names = NA) ###### ## Quality control steps ## qc <- qc(raw.data, x.mas) # Calls "qc" function which generates object containing scale factors, GAPDH-Actin_3'/5' ratios, target intensity, % present, average, min, max, mean background int and more; for more details type "?qc". x.mas@description@preprocessing # Creates QC summary. pdf(file="qc.pdf"); plot(qc); dev.off() # Creates summary plot of QC data. See description on page 5 of vignette "simpleaffy". ######