set.seed(2014) library(limma) library(edgeR) d=c(1,1.2,1.5,2,3,4,8) libsize=c(0.5,0.6,0.7,0.8,0.9,1) testfc=c(1.3,1.5,2,3,4) savelibsize=c(5,6,7,8,9,10) savetestfc=c(13,15,2,3,4) for(mm in 1:6){ x <- read.delim("forCharity_RawReadCounts.txt") x <- as.matrix(x[,-1]) ngenes <- 10000 for(pp in 1:5){ cat(pp,"") u <- runif(100) nsim <- 100 fd =fdpure=fdvfremove=vfda2=vfdanew2=fdonly=matrix(0,ngenes,length(d)) nd =ndpure=ndvfremove=vnda2=vndanew2=ndonly=matrix(0,nsim,length(d)) nd5 =ndpure5=ndvfremove5=vnda25=vndanew25=nd5only=rep(0,length(d)) aw=aw1=aw2=awold2=awnew2=awonly=matrix(NA,length(d)*nsim,10) u <- runif(100) for(j in 1:length(d)){ cat(j,"") for (k in 1:nsim) { prop <- goodTuringProportions(x) n1 <- 5 n2 <- 5 nlibs <- n1+n2 expected.lib.size <-c(rep(20e6,9),20e6*libsize[mm]) is.expr <- rowSums(prop>1e-6)>=6 prop <- prop[is.expr,] i <- seq(from=0,to=1,length=length(prop)) f <- approxfun(i,sort(prop),rule=2) baselineprop <- f( (1:ngenes)/(ngenes+1) ) baselineprop <- baselineprop/sum(baselineprop) i <- sample(1:ngenes,200) i1 <- i[1:100] #seprate two groups for simulate groups factor in further analysis # i2 <- i[101:200] fc <- testfc[pp] baselineprop1 <- baselineprop2 <- baselineprop #calculate baseline propration for each gene baselineprop1[i1] <- baselineprop1[i1]*fc #add fc as fold changes# baselineprop2[i2] <- baselineprop2[i2]*fc #same as above# mu0.1 <- matrix(baselineprop1,ngenes,1) %*% matrix(expected.lib.size[1:n1],1,n1) mu0.2 <- matrix(baselineprop2,ngenes,1) %*% matrix(expected.lib.size[(n1+1):(n1+n2)],1,n2) #mu0.1 <- matrix(baselineprop,ngenes,1) %*% matrix(expected.lib.size[1:n1],1,n1) #mu0.2 <- matrix(baselineprop,ngenes,1) %*% matrix(expected.lib.size[(n1+1):(n1+n2)],1,n2) mu0 <- cbind(mu0.1,mu0.2) status <- rep(0,ngenes) status[i1] <- -1 status[i2] <- 1 BCV0 <- 0.12+1/sqrt(mu0) df.BCV <- 5 BCV <- BCV0*sqrt(df.BCV/rchisq(ngenes,df=df.BCV)) if(NCOL(BCV)==1) BCV <- matrix(BCV,ngenes,nlibs) shape <- 1/BCV^2 squareBCV=BCV^2 squareBCV1=squareBCV squareBCV1[,10]=squareBCV[,10]*d[j] shape1 <- 1/squareBCV1 scale <- mu0/shape1 mu <- matrix(rgamma(ngenes*nlibs,shape=shape1,scale=scale),ngenes,nlibs) #random generation gamma distribution: counts <- matrix(rpois(ngenes*nlibs,lambda=mu),ngenes,nlibs) #random poisson distribution, colSums(counts) keep <- rowSums(counts)>=10 nkeep <- sum(keep) counts2 <- counts[keep,] # Design group <- factor(c(1,1,1,1,1,2,2,2,2,2)) design <- model.matrix(~group) y <- voom(counts2,design,plot=FALSE) ypure=y ypure$weights=NULL y$genes$Status <- status[keep] ######################### voom + no weights####################### ypure$genes$Status <- status[keep] vfpure = lmFit(ypure, design) vfpure = eBayes(vfpure) opure <- order(vfpure$lods[,2], decreasing=TRUE) rankingpure <- rep(0, ngenes) rankingpure[1:nkeep] <- status[keep][opure] rankingpure[(nkeep+1):ngenes] <- status[!keep] fdpure[,j] <- fdpure[,j] + cumsum(1-abs(rankingpure)) ndpure[k,j]<- sum(p.adjust(vfpure$p.value[,2],method="BH")<0.1) # ndpure5[j]<- ndpure5[j] + sum(vfpure$p.value[,2]<0.01) ######################### voom####################### y$genes$Status <- status[keep] vf = lmFit(y, design) vf = eBayes(vf) o <- order(vf$lods[,2], decreasing=TRUE) ranking <- rep(0, ngenes) ranking[1:nkeep] <- status[keep][o] ranking[(nkeep+1):ngenes] <- status[!keep] fd[,j] <- fd[,j] + cumsum(1-abs(ranking)) nd[k,j]<- sum(p.adjust(vf$p.value[,2],method="BH")<0.1) #nd5[j]<- nd5[j] + sum(vf$p.value[,2]<0.01) ######################### v weights matrix only####################### y <- voom(counts2,design,plot=FALSE) yonly=y yonly$weights=matrix(1,nrow(counts2),ncol(counts2)) yonly$genes$Status <- status[keep] awonly[k+(j-1)*nsim,]=arrayWeights(yonly,design) wtsonly=(asMatrixWeights(awonly[k+(j-1)*nsim,], dim(yonly)))*(yonly$w) attr(wtsonly, "arrayweights") = NULL vfonly = lmFit(yonly, design,weights=wtsonly) vfonly = eBayes(vfonly) oonly <- order(vfonly$lods[,2], decreasing=TRUE) rankingonly <- rep(0, ngenes) rankingonly[1:nkeep] <- status[keep][oonly] rankingonly[(nkeep+1):ngenes] <- status[!keep] fdonly[,j] <- fdonly[,j] + cumsum(1-abs(rankingonly)) ndonly[k,j]<- sum(p.adjust(vfonly$p.value[,2],method="BH")<0.1) #nd5only[j]<- nd5only[j] + sum(vfonly$p.value[,2]<0.01) ######################### voom + sample weights ####################### ya<-voomWithQualityWeights(counts2, design, normalization="none",plot=F) vfita2= lmFit(ya, design) vfita2=eBayes(vfita2) o <- order(vfita2$lods[,2], decreasing=TRUE) vrankinga2 <- rep(0, ngenes) vrankinga2[1:nkeep] <- status[keep][o] vrankinga2[(nkeep+1):ngenes] <- status[!keep] vfda2[,j] <- vfda2[,j]+ cumsum(1-abs(vrankinga2)) vnda2[k,j] <- sum(p.adjust(vfita2$p.value[,2],method="BH")<0.1) # vnda25[j] <- vnda25[j] + sum(vfita2$p.value[,2]<0.01) ######################### voom + block weights ####################### Z <- cbind(c(1,1,1,1,1,0,0,0,0,-1),c(0,0,0,0,0,1,1,1,1,-1),c(0,0,0,0,0,0,0,0,0,1)) ynew2<-voomWithQualityWeights(counts2, design,normalization="none", var.design=Z, plot=F) vfitanew2= lmFit(ynew2, design) vfitanew2=eBayes(vfitanew2) o <- order(vfitanew2$lods[,2], decreasing=TRUE) vrankinganew2 <- rep(0, ngenes) vrankinganew2[1:nkeep] <- status[keep][o] vrankinganew2[(nkeep+1):ngenes] <- status[!keep] vfdanew2[,j] <- vfdanew2[,j]+ cumsum(1-abs(vrankinganew2)) vndanew2[k,j] <- sum(p.adjust(vfitanew2$p.value[,2],method="BH")<0.1) # vndanew25[j] <- vndanew25[j] + sum(vfitanew2$p.value[,2]<0.01) ######remove the outlier###### keep <- rowSums(counts)>=10 nkeep <- sum(keep) counts2 <- counts[keep,][,1:9] # Design group <- factor(c(1,1,1,1,1,2,2,2,2)) design <- model.matrix(~group) y <- voom(counts2,design,plot=FALSE) y$genes$Status <- status[keep] vfremove = lmFit(y, design) vfremove = eBayes(vfremove) o <- order(vfremove$lods[,2], decreasing=TRUE) ranking <- rep(0, ngenes) ranking[1:nkeep] <- status[keep][o] ranking[(nkeep+1):ngenes] <- status[!keep] fdvfremove[,j] <- fdvfremove[,j] + cumsum(1-abs(ranking)) ndvfremove[k,j] <- sum(p.adjust(vfremove$p.value[,2],method="BH")<0.1) # ndvfremove5[j] <- ndvfremove5[j]+ sum(vfremove$p.value[,2]<0.01) } } save(list=ls(),file=paste("finalweight5vs5",savetestfc[pp],savelibsize[mm],".rda",sep="")) } }