################################################### ### chunk number 1: loaddata ################################################### #line 101 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" library(GeneMeta) library(RColorBrewer) #load("~/Bioconductor/Projects/GraphCombine/MetaBreast/data/Nevins.RData") data(Nevins) ################################################### ### chunk number 2: Splitting ################################################### #line 119 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" 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] ################################################### ### chunk number 3: phenoData ################################################### #line 132 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" #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)) ################################################### ### chunk number 4: calcd ################################################### #line 151 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" #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)) ################################################### ### chunk number 5: calcQ ################################################### #line 173 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" #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) ################################################### ### chunk number 6: Qplo ################################################### #line 182 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" hist(my.Q,breaks=50,col="red") ################################################### ### chunk number 7: grapics ################################################### #line 190 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" ########### 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") ################################################### ### chunk number 8: FEMmodel ################################################### #line 209 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" muFEM = mu.tau2(mymns, myvars) sdFEM = var.tau2(myvars) ZFEM = muFEM/sqrt(sdFEM) ################################################### ### chunk number 9: qnormPlotFEM ################################################### #line 218 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" qqnorm(ZFEM,pch="*") qqline(ZFEM,col="red") ################################################### ### chunk number 10: estimatedEffects ################################################### #line 229 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" 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) ################################################### ### chunk number 11: compMU ################################################### #line 254 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" plot(muFEM, muREM, pch=".") abline(0,1,col="red") ################################################### ### chunk number 12: theTau ################################################### #line 262 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" hist(my.tau2.DL,col="red",breaks=50,main="Histogram of tau") ################################################### ### chunk number 13: claculationAllinOne ################################################### #line 276 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" 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) ################################################### ### chunk number 14: show results ################################################### #line 287 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" theScores[1:2,] ################################################### ### chunk number 15: originalvsCombination ################################################### #line 314 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" plot(theScores[,"zSco_Ex_3"],theScores[,"zSco"],pch="*",xlab="original score",ylab="combined score") abline(0,1,col="red") ################################################### ### chunk number 16: IDRplot ################################################### #line 333 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" IDRplot(theScores,Combine=1:2,colPos="blue", colNeg="red") ################################################### ### chunk number 17: calculationAllinOne ################################################### #line 345 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" ScoresFDR <- zScoreFDR(esets, classes, useREM=FALSE, nperm=50,CombineExp=1:2) ################################################### ### chunk number 18: calculation whole set ################################################### #line 351 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" names(ScoresFDR) ################################################### ### chunk number 19: showSet ################################################### #line 361 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" ScoresFDR$pos[1:2,] ################################################### ### chunk number 20: gettingFDR ################################################### #line 371 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" FDRwholeSettwo <- sort(ScoresFDR$"two.sided"[,"FDR"]) experimentstwo <- list() for(j in 1:3){ experimentstwo[[j]] <- sort(ScoresFDR$"two.sided"[,paste("FDR_Ex_",j,sep="")]) } ################################################### ### chunk number 21: bb2 ################################################### #line 379 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" theNewC <- brewer.pal(3,"Set2") ################################################### ### chunk number 22: FDRtwo ################################################### #line 384 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" ##################### # # #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])) ################################################### ### chunk number 23: theFDRplots ################################################### #line 408 "vignettes/GeneMeta/inst/doc/GeneMeta.Rnw" #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)