library(affy)
vignette("affy")

setwd("C:\\z\\Courses\\S416\\GEOdata\\Batch1CEL")

D=ReadAffy()
D
attributes(D)
sampleNames(D)

hist(D)
hist(D[,1:3])

boxplot(D)

boxplot(D,las=3)
help(par)

boxplot(D,axes=F,xlab="GeneChip",ylab="log2 PM")
axis(1,labels=substr(sampleNames(D),8,9),at=1:12)
axis(2)
box()

image(D[,1])

gn=geneNames(D)
gn[1:10]

pm(D, gn[10])

cbind(pm(D, gn[10])[,1],mm(D, gn[10])[,1])

plot(rep(1:11,2),c(pm(D, gn[10])[,1],mm(D, gn[10])[,1]),
     pch=rep(c("p","m"),each=11),col=(rep(c(4,2),each=11)),
     ylim=c(0,1000),xlab="Probe",ylab="Intensity")

mean(pm(D)<mm(D))

library(geneplotter)
MAplot(D[,c(1,2,10)], pairs = TRUE, plot.method = "smoothScatter")

r=rma(D)
r

d=exprs(r)
dim(d)
head(d)

#write.exprs(r, file="rmabatch1.txt")

diet=gl(4,3)
diet
#1=chimp, 2=McDonald's, 3=cafe, 4=mouse

boxplot(d)
boxplot(as.data.frame(d))

boxplot(as.data.frame(d),col=diet+1)

get.p=function(y)
{
  out=lm(y~diet)
  a=anova(out)
  df=a[2,1]
  b=coef(out)
  v=vcov(out)
  m=c(0,1,0,0)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p12=2*(1-pt(abs(tstat),df))
  m=c(0,0,1,0)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p13=2*(1-pt(abs(tstat),df))
  m=c(0,0,0,1)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p14=2*(1-pt(abs(tstat),df))
  m=c(0,-1,1,0)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p23=2*(1-pt(abs(tstat),df))
  m=c(0,-1,0,1)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p24=2*(1-pt(abs(tstat),df))
  m=c(0,0,-1,1)
  tstat=t(m)%*%b/sqrt(t(m)%*%v%*%m)
  p34=2*(1-pt(abs(tstat),df))
  p=c(p12,p13,p14,p23,p24,p34)
  p
}

get.p(d[1,])

#p=t(apply(d,1,get.p))


hist(p[,1],col=4,xlab="p-value",ylab="Number of Genes",
     main="Chimp vs. McDonald's")
x11()
hist(p[,4],col=3,xlab="p-value",ylab="Number of Genes",
     main="McDonald's vs. Cafe")

ls()
source("c:/z/Courses/S416/R/multtest.txt")

estimate.m0(p[,1])

nrow(d)-estimate.m0(p[,1])


nrow(d)-estimate.m0(p[,4])

q1=jabes.q(p[,1])

sum(q1<=0.05)
sum(q1<=0.15)
sum(q1<=0.20)

plot(p[,1],q1,xlab="p-value",ylab="q-value",
     ylim=c(0,max(q1)),main="Chimp vs. McDonald's")

q4=jabes.q(p[,4])
sum(q4<=0.20)
min(q4)

qvals=apply(p,2,jabes.q)

apply(qvals<=0.05,2,sum)

plot(p[,3],qvals[,3],xlab="p-value",ylab="q-value",
     ylim=c(0,max(qvals[,3])),main="Chimp vs. Mouse")

plot(p[p[,3]<=0.05,3],qvals[p[,3]<=0.05,3],xlab="p-value",ylab="q-value",
     main="Chimp vs. Mouse")

fdrcuts=seq(.01,0.1,by=0.01)
cbind(fdrcuts,NumberOfGenes=apply(outer(qvals[,3],fdrcuts,"<="),2,sum))

nrow(d)-estimate.m0(p[,3])


out=ub.mix(p[,3])
out

plot.ub.mix(p[,3],out)


out=ub.mix(p[,1])
out

plot.ub.mix(p[,1],out)

out=ub.mix(p[,1],c(.75,.5,3))
out

plot.ub.mix(p[,1],out)

ppde1=ppde(p[,1],out)

plot(p[,1],ppde1,xlab="p-value",ylab="Posterior Probability of DE",
     main="Chimp vs. McDonald's")

ppdecuts=seq(.6,.95,by=0.05)
cbind(ppdecuts,NumberOfGenes=apply(outer(ppde1,ppdecuts,">="),2,sum))

sum(q1<=.20)

row.names(d)[q1<=0.2]

library(limma)

design=model.matrix(~diet)
colnames(design)=c("mu","a2","a3","a4")

fit=lmFit(d,design)

contr.mat=makeContrasts(a2,a3,a4,a3-a2,a4-a2,a4-a3,levels=design)

fit2=contrasts.fit(fit,contr.mat)

fit3=eBayes(fit2)

attributes(fit3)

fit3$df.prior
fit3$s2.prior

plot(log(fit3$sigma),log(sqrt(fit3$s2.post)),
     xlab="Log Sqrt Original Variance Estimate",
     ylab="Log Sqrt Empirical Bayes Variance Estimate")
lines(c(-99,99),c(-99,99),col=2)

lsrs20=log(sqrt(fit3$s2.prior))
abline(h=lsrs20,col=4)
abline(v=lsrs20,col=4)

dim(fit3$p.value)

plot(p[,1],fit3$p.value[,1],
     xlab="Original P-Value",
     ylab="Empirical Bayes P-Value",
     main="Chimp vs. McDonald's")
lines(c(-99,99),c(-99,99),col=2)

limmaq=apply(fit3$p.value,2,jabes.q)

fdrcuts=seq(.01,0.1,by=0.01)
cbind(fdrcuts,NumberOfGenes=apply(outer(qvals[,1],fdrcuts,"<="),2,sum))

cbind(fdrcuts,NumberOfGenes=apply(outer(limmaq[,1],fdrcuts,"<="),2,sum))