GeneMeta.R
In GeneMeta: MetaAnalysis for High Throughput Experiments

### R code from vignette source 'GeneMeta.Rnw'

###################################################
### code chunk number 1: loaddata
###################################################
library(GeneMeta)
library(RColorBrewer)
#load("~/Bioconductor/Projects/GraphCombine/MetaBreast/data/Nevins.RData")
data(Nevins)


###################################################
### code chunk number 2: Splitting
###################################################
set.seed(1609)
thestatus  <- pData(Nevins)[,"ER.status"]
group1     <- which(thestatus=="pos")
group2     <- which(thestatus=="neg")
rrr        <- c(sample(group1, floor(length(group1)/2)),
		sample(group2,ceiling(length(group2)/2)))
Split1     <- Nevins[,rrr]
Split2     <- Nevins[,-rrr]


###################################################
### code chunk number 3: phenoData
###################################################
#obtain classes
Split1.ER<-pData(Split1)[,"ER.status"]
levels(Split1.ER) <- c(0,1)
Split1.ER<- as.numeric(as.character(Split1.ER))
Split2.ER<-pData(Split2)[,"ER.status"]
levels(Split2.ER) <- c(0,1)
Split2.ER<- as.numeric(as.character(Split2.ER))


###################################################
### code chunk number 4: calcd
###################################################
#calculate d for Split1
d.Split1     <- getdF(Split1, Split1.ER)

#adjust d value
d.adj.Split1    <- dstar(d.Split1, length(Split1.ER))
var.d.adj.Spli1 <- sigmad(d.adj.Split1, sum(Split1.ER==0), sum(Split1.ER==1))


#calculate d for Split2
d.Split2 <- getdF(Split2, Split2.ER)

#adjust d value
d.adj.Split2     <- dstar(d.Split2, length(Split2.ER))
var.d.adj.Split2 <- sigmad(d.adj.Split2, sum(Split2.ER==0), sum(Split2.ER==1))



###################################################
### code chunk number 5: calcQ
###################################################
#calculate Q
mymns  <- cbind(d.adj.Split1, d.adj.Split2)
myvars <- cbind(var.d.adj.Spli1,var.d.adj.Split2)
my.Q   <- f.Q(mymns, myvars)
mean(my.Q)


###################################################
### code chunk number 6: Qplo
###################################################
hist(my.Q,breaks=50,col="red")


###################################################
### code chunk number 7: grapics
###################################################
########### graphics ################

num.studies<-2

#quantiles of the chisq distribution
chisqq <- qchisq(seq(0, .9999, .001), df=num.studies-1)
tmp<-quantile(my.Q, seq(0, .9999, .001))
qqplot(chisqq, tmp, ylab="Quantiles of Sample",pch="*",
    xlab="Quantiles of Chi square", main="QQ Plot")
lines(chisqq, chisqq, lty="dotted",col="red")


###################################################
### code chunk number 8: FEMmodel
###################################################
 muFEM = mu.tau2(mymns, myvars)
 sdFEM = var.tau2(myvars)
 ZFEM = muFEM/sqrt(sdFEM)


###################################################
### code chunk number 9: qnormPlotFEM
###################################################
 qqnorm(ZFEM,pch="*")
 qqline(ZFEM,col="red")


###################################################
### code chunk number 10: estimatedEffects
###################################################
my.tau2.DL<-tau2.DL(my.Q, num.studies, my.weights=1/myvars)

#obtain new variances s^2+tau^2
myvarsDL <- myvars + my.tau2.DL


#compute
muREM <- mu.tau2(mymns, myvarsDL)


#cumpute mu(tau)
varREM <- var.tau2(myvarsDL)

ZREM <- muREM/sqrt(varREM)


###################################################
### code chunk number 11: compMU
###################################################
plot(muFEM, muREM, pch=".")
abline(0,1,col="red")


###################################################
### code chunk number 12: theTau
###################################################
hist(my.tau2.DL,col="red",breaks=50,main="Histogram of tau")


###################################################
### code chunk number 13: claculationAllinOne
###################################################
esets     <- list(Split1,Split2,Nevins)
data.ER   <-pData(Nevins)[,"ER.status"]
levels(data.ER) <- c(0,1)
data.ER<- as.numeric(as.character(data.ER))
classes   <- list(Split1.ER,Split2.ER,data.ER)
theScores <- zScores(esets,classes,useREM=FALSE,CombineExp=1:2)


###################################################
### code chunk number 14: show results
###################################################
theScores[1:2,]


###################################################
### code chunk number 15: originalvsCombination
###################################################
plot(theScores[,"zSco_Ex_3"],theScores[,"zSco"],pch="*",xlab="original score",ylab="combined score")
abline(0,1,col="red")


###################################################
### code chunk number 16: IDRplot
###################################################
IDRplot(theScores,Combine=1:2,colPos="blue", colNeg="red")


###################################################
### code chunk number 17: calculationAllinOne
###################################################
ScoresFDR <- zScoreFDR(esets, classes, useREM=FALSE, nperm=50,CombineExp=1:2)


###################################################
### code chunk number 18: calculation whole set
###################################################
names(ScoresFDR)


###################################################
### code chunk number 19: showSet
###################################################
ScoresFDR$pos[1:2,]


###################################################
### code chunk number 20: gettingFDR
###################################################
FDRwholeSettwo <- sort(ScoresFDR$"two.sided"[,"FDR"])
experimentstwo <- list()
for(j in 1:3){
experimentstwo[[j]] <- sort(ScoresFDR$"two.sided"[,paste("FDR_Ex_",j,sep="")])
}


###################################################
### code chunk number 21: bb2
###################################################
theNewC <- brewer.pal(3,"Set2")


###################################################
### code chunk number 22: FDRtwo
###################################################
#####################
#                   #
#two sided z values #
#                   #
#####################

plot(FDRwholeSettwo,pch="*",col="red",ylab="FDR",xlab="Number of genes")
for(j in 1:3)
points(experimentstwo[[j]],pch="*", col=theNewC[j])
legend(4000,0.4,c("Combined set","Split 1" , "Split 2" ,"original set"), c("red",theNewC[1:3]))


###################################################
### code chunk number 23: theFDRplots
###################################################
#par(mfrow=c(2,2))
#CountPlot(ScoresFDR,Score="FDR",kindof="neg",cols=c("red",theNewC),
#	main="Negative FDR", xlab="FDR threshold", ylab="Number of genes",CombineExp=1:2)
#CountPlot(ScoresFDR,Score="FDR",cols=c("red",theNewC),kindof="pos",
#main="Positive FDR", xlab="FDR threshold", ylab="Number of genes",Combine=1:2)
CountPlot(ScoresFDR,Score="FDR",kindof="two.sided",cols=c("red",theNewC),main="two sided FDR", xlab="FDR threshold",ylab="Number of genes",CombineExp=3)