###need to run with wd having the files for hwk2
mypar <- function(a=1,b=1){
  par(mar=c(2.5,2.5,1.6,1.1),mgp=c(1.5,.5,0))
  par(mfrow=c(a,b))
}

library(Biobase)
library(geneplotter)
library(genefilter)
eset = read.exprSet(exprs="exprs.txt",phenoData="phenoData.txt",description="miame.txt",annotation="annotation.txt")

##later we will be using log exclusively
exprs(eset) <- log2(exprs(eset))

bitmap("figure-01.png",res=300,pointsize=20)
mypar(1,1)
plot(2^exprs(eset)[,1],2^exprs(eset)[,2],log="xy",xlab="Expression on array1",ylab="Expression on array 2",cex=.25,pch=16,yaxt="n")
axis(2,at=10^c(1,3,5),label=c("1","100","10000"))
##hack to make lines at FC=2 or 1/2
lines(2^seq(-2,16),2^seq(-1,17),col="red")
lines(2^seq(-1,17),2^seq(-2,16),col="red")
dev.off()

bitmap("figure-02.png",res=300,pointsize=20)
mypar(1,1)
plot(2^exprs(eset)[,1],2^exprs(eset)[,2],xlab="Expression on array1",ylab="Expression on array 2",cex=.25,pch=16)
abline(1024,-1)
dev.off()

bitmap("figure-03.png",res=300,pointsize=20)
mypar(1,1)
x=rnorm(10000,0,1/4)
hist(x,nclass=35,col="grey",main="Gene 1",xlim=c(-2.5,2.5))
dev.off()

bitmap("figure-04.png",res=300,pointsize=20,width=12)
y=rnorm(10000,0,1)
mypar(1,2)
hist(x,nclass=35,col="grey",main="Gene 1",xlim=c(-2.5,2.5))
hist(y,nclass=35,col="blue",main="Gene 2",xlim=c(-2.5,2.5))
dev.off()

x=exprs(eset)[,1]
y=exprs(eset)[,2]

bitmap("figure-05.png",res=300,pointsize=20)
mypar(1,1)
plot(x,y,xlab="Log expression on array 1",ylab="Log expression on array 2",cex=.25,pch=16)
abline(-1,1,col="red")
abline(1,1,col="red")
dev.off()

M=y-x
A=(x+y)/2

bitmap("figure-06.png",res=300,pointsize=20)
mypar(1,1)
plot(A,M,cex=.25,pch=16)
YLIM=range(M) ##saving for later
abline(h=-1,col="red")
abline(h=1,col="red")
dev.off()

tt = rowttests(eset,eset$pop)

M=tt$dm
A=rowMeans(exprs(eset))

bitmap("figure-07.png",res=300,pointsize=20)
mypar(1,1)
plot(A,M,ylim=YLIM,cex=.25,pch=16)
abline(h=-1,col="red")
abline(h=1,col="red")
dev.off()

bitmap("figure-08.png",res=300,pointsize=20)
mypar(1,1)
smoothScatter(A,M,ylim=YLIM)
abline(h=-1,col="red")
abline(h=1,col="red")
dev.off()


##three genes. i used interactive identify on the volcano thuse code not shown
Index=c( 4977,10771,3238)
x=rep(c(0.9,1.1,1.9,2.1,2.9,3.1),rep(3,6))
cols=rep(c(2,4,2,4,2,4),rep(3,6))
y=as.vector(t(exprs(eset)[Index,]))

bitmap("figure-09.png",res=300,pointsize=20,width=9)
mypar(1,1)
plot(jitter(x),y,col=cols,pch=16,xaxt="n",xlim=c(.7,3.3),ylab="Log expression",xlab="")
axis(1,1:3,c("Gene 1","Gene 2","Gene 3"))
abline(v=1:2+.5)
dev.off()
YLIM=c(-4,4)

bitmap("figure-10.png",res=300,pointsize=20,width=12)
mypar(1,2)
plot(A[-Index],M[-Index],cex=.25,pch=16,ylim=YLIM,xlab="A",ylab="M")
text(A[Index],M[Index],1:3,col="red",cex=1.1)
abline(h=c(-1,1),col="red")
plot(M[-Index],-log10(tt$p.value)[-Index],xlim=YLIM,ylab="-Log (base 10) p-value",cex=.25,pch=16,xlab="M")
text(M[Index],-log10(tt$p.value)[Index],1:3,col="red",cex=1.1)
abline(v=c(-1,1),col="red")
abline(h=c(-2,2),col="red")
dev.off()