library(limma)
data.directory<-"http://eh3.uc.edu/teaching/cfg/2006/data/"
targets<-readTargets("http://eh3.uc.edu/teaching/cfg/2006/data/NTargets.txt")
spottypes<-readSpotTypes("http://eh3.uc.edu/teaching/cfg/2006/data/SpotTypes.txt")
LimmaDataNickel<-read.maimages(files=targets$experiment,source="genepix", path = data.directory, 
        names=paste("Cy5_",targets$cy5,"_VS_Cy3_",targets$cy3,sep=""),
	columns=list(Gf = "F532 Median",Gb ="B532 Median", Rf = "F635 Median", Rb = "B635 Median"),
	annotation=c("Name","ID","Block","Row","Column"),wt.fun=wtflags(0))

LimmaDataNickel$weights[LimmaDataNickel$genes$ID=="Blank"]<-0
LimmaDataNickel$weights[LimmaDataNickel$genes$ID=="empty"]<-0
LimmaDataNickel$printer<-getLayout(LimmaDataNickel$genes)
LimmaDataNickel$genes$Status<-controlStatus(spottypes,LimmaDataNickel)
NBLimmaDataNickel<-backgroundCorrect(LimmaDataNickel,method="none")
NNBLimmaDataNickel<-normalizeWithinArrays(NBLimmaDataNickel,method="none")
design<-c(-1,1,1)
LimmaFit<-lmFit(NNBLimmaDataNickel,design)
complete.data<-apply(NNBLimmaDataNickel$weights,1,sum)==3
LimmaFitTstat<-abs(LimmaFit$coefficients/(LimmaFit$sigma*LimmaFit$stdev.unscaled))
LimmaFitPvalue<-2*pt(LimmaFitTstat,LimmaFit$df.residual,lower.tail=FALSE)

N<-sum(!is.na(LimmaFitPvalue))
FDRpvalue<-p.adjust(LimmaFitPvalue,method="fdr")
BONFpvalue<-p.adjust(LimmaFitPvalue,method="bonferroni")
Qvalue<-(N/rank(LimmaFitPvalue))*LimmaFitPvalue

All.Sigs<-data.frame(list(p.value=LimmaFitPvalue,Qvalue=Qvalue,FDRpvalue=FDRpvalue,BONFpvalue=BONFpvalue))
All.Sigs[1:10,]
OrderingPvalues<-order(LimmaFitPvalue)
All.Sigs<-All.Sigs[OrderingPvalues,]
All.Sigs[1:10,]

lowess(x=NNBLimmaDataNickel$A,y=NNBLimmaDataNickel$M)
NickelLoess<-loess(NNBLimmaDataNickel$M[,1]~NNBLimmaDataNickel$A[,1])
attributes(NickelLoess)
length(NickelLoess$fitted)
length(NNBLimmaDataNickel$A
NickelLoess$fitted[1:10]
NickelLoess$x[1:10]

par(mfrow=c(1,1))
plotMA(NNBLimmaDataNickel,array=1,zero.weights=F,legend=F)
lines(c(-5,100),c(0,0),lw=5,col="blue")
points(NNBLimmaDataNickel$A[,1],NickelLoess$fitted,col="red")

NLNBLimmaDataNickel<-normalizeWithinArrays(NBLimmaDataNickel,method="loess")
par(mfrow=c(1,1))
plotMA(NLNBLimmaDataNickel,array=1,zero.weights=F,legend=F)
lines(c(-5,100),c(0,0),lw=5,col="blue")

LoessLimmaFit<-lmFit(NLNBLimmaDataNickel,design)
LoessLimmaFitTstat<-abs(LoessLimmaFit$coefficients/(LoessLimmaFit$sigma*LoessLimmaFit$stdev.unscaled))
LoessLimmaFitPvalue<-2*pt(LoessLimmaFitTstat,LoessLimmaFit$df.residual,lower.tail=FALSE)
LoessFDRpvalue<-p.adjust(LoessLimmaFitPvalue,method="fdr")

plot(LimmaFit$coefficients[complete.data],-log(FDRpvalue[complete.data],base=10),main="Not Normalized",
xlab="Mean Difference",ylab="-log10(FDRpvalue)",ylim=c(0,0.6))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = NULL, equilogs = TRUE)

plot(LoessLimmaFit$coefficients[complete.data],-log(LoessFDRpvalue[complete.data],base=10),main="Loess Normalized",
xlab="Mean Difference",ylab="-log10(FDRpvalue)",ylim=c(0,0.6))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = NULL, equilogs = TRUE)

EBayesLoessLimmaFit<-eBayes(LoessLimmaFit)
attributes(EBayesLoessLimmaFit)
EBayesFDR<-p.adjust(EBayesLoessLimmaFit$p.value,method="fdr")

plot(LoessLimmaFit$coefficients[complete.data],-log(LoessFDRpvalue[complete.data],base=10),main="Simple Paired t-test",
xlab="Mean Difference",ylab="-log10(FDRpvalue)",ylim=c(0,3))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = NULL, equilogs = TRUE)

plot(EBayesLoessLimmaFit$coefficients[complete.data],-log(EBayesFDR[complete.data],base=10),main="Empirical Bayes t-test",
xlab="Mean Difference",ylab="-log10(FDRpvalue)",ylim=c(0,3))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = NULL, equilogs = TRUE)

LinearModelData<-cbind(NLNBLimmaDataNickel,NLNBLimmaDataNickel)
dim(LinearModelData)
LinearModelData$M[1,]
LinearModelData$M<-cbind((NLNBLimmaDataNickel$A+NLNBLimmaDataNickel$M/2),NLNBLimmaDataNickel$A-(NLNBLimmaDataNickel$M/2))
LinearModelData$M[1,]

designMatrix<-cbind(Ctl=c(1,0,0,0,1,1),Nic=c(0,1,1,1,0,0))
designMatrix
contrast<-c(-1,1)

FitLM<-lmFit(LinearModelData,designMatrix)
FitContrast<-contrasts.fit(FitLM,contrast)
attributes(FitContrast)
ContrastTstat<-abs(FitContrast$coefficients)/(FitContrast$sigma*FitContrast$stdev.unscaled)
ContrastPvalue<-2*pt(ContrastTstat,FitContrast$df.residual,lower.tail=FALSE)

MNic<-apply(LinearModelData$M[,Nic],1,mean,na.rm=TRUE)
VNic<-apply(LinearModelData$M[,Nic],1,var,na.rm=TRUE)
MCtl<-apply(LinearModelData$M[,Ctl],1,mean,na.rm=TRUE)
VCtl<-apply(LinearModelData$M[,Ctl],1,var,na.rm=TRUE)
NNic<-apply(!is.na(LinearModelData$M[,Nic]),1,sum,na.rm=TRUE)
NCtl<-apply(!is.na(LinearModelData$M[,Ctl]),1,sum,na.rm=TRUE)

VNicCtl<-(((NNic-1)*VNic)+((NCtl-1)*VCtl))/(NCtl+NNic-2)

DF<-NNic+NCtl-2

TStat<-abs(MNic-MCtl)/((VNicCtl*((1/NNic)+(1/NCtl)))^0.5)
TPvalue<-2*pt(TStat,DF,lower.tail=FALSE)
TStat[1]
TPvalue[1]

par(mfrow=c(1,2))
plot(ContrastTstat[complete.data],TStat[complete.data])
plot(FitContrast$coefficients[complete.data],(MNic-MCtl)[complete.data])
plot(-log10(ContrastPvalue[complete.data]),-log10(TPvalue[complete.data]))