setwd("/Users/cliu/Desktop/comparison/Seq-weights/final/3vs3") d=c(1,1.2,1.5,2,3,4,8) filenames=dir(pattern="final")[1:30] ROCname=c("voomonly","puremodel","remove","aw2","newaw2","weightsonly") titles1=rep("NA",30) titles1[12]=c("(A)") titles1[18]=c("(B)") titles1[30]=c("(C)") pchs=c(1,2,20,8,5,6) files1=apply(as.matrix(filenames),1,function(x){strsplit(x,".rda")[[1]][1]}) listnames=matrix(NA,5,length(filenames)) for(i in 1:nrow(listnames)){ listnames[i,]=paste(files1,ROCname[i],sep="_") } names=rbind(listnames,listnames,listnames,listnames,listnames,listnames,listnames) dname=as.vector(rbind(d,d,d,d,d)) finalname=apply(names,2,function(x){paste(x,dname,sep="_")}) dim(finalname) ROCresult=matrix(NA,42,ncol(finalname)) #fill information by column# special=c(7:18,25:30) for(i in 1:ncol(ROCresult)){ cat(i) tmp1=filenames[i] load(tmp1) i=which(tmp1==filenames) #tmp1=filenamesall[i+30] # load(tmp1) fd =fd/nsim fdpure=fdpure/nsim fdonly=fdonly/nsim fdvfremove=fdvfremove/nsim vfda2=vfda2/nsim vfdanew2=vfdanew2/nsim for( j in 1:length(d)){ ROCresult[6*(j-1)+1,i]=cumsum(fd[1:200,j])[200] ROCresult[6*(j-1)+2,i]=cumsum(fdpure[1:200,j])[200] ROCresult[6*(j-1)+3,i]=cumsum(vfda2[1:200,j])[200] ROCresult[6*(j-1)+4,i]=cumsum(vfdanew2[1:200,j])[200] ROCresult[6*(j-1)+5,i]=cumsum(fdvfremove[1:200,j])[200] ROCresult[6*(j-1)+6,i]=cumsum(fdonly[1:200,j])[200] } } pdf("FC-3sel.pdf",height=3,width=8) par(mfrow=c(1,3),mai=c(0.55,0.55,0.25,0.05)) group=c(rep(1:6,7)) names=c(rep(1.3,6),rep(1.5,6),rep(2,6),rep(3,6),rep(4,6)) sizes=rep(c(0.5,0.6,0.7,0.8,0.9,1),5) n=c(12,18,30) for(i in n){ plot(c(1,7),c(ceiling(min(ROCresult[,i]))-200,ceiling(max(ROCresult[,i]))+20),type="n",ylab="cumulative False discoveries",xlab="Sample variability",axes=F,main=titles1[i]) axis(1, at=c(1,2,3,4,5,6,7), labels=d) axis(2) colssel=c("red","blue","green","orange","black","purple") pchs=c(8,2,20,6,5,1) cexs=c(1.3,1.3,0.8,1.3,1.3,1.3) for(j in 1:6){ points(ROCresult[group==j,i],col=colssel[j],pch=pchs[j],cex=cexs[j]) lines(ROCresult[group==j,i],col=colssel[j],lwd=1) } } legend("topleft",legend=c("No Weights","Voom","Sample","Voom + Sample","Voom + Block","Sample Removal"),text.col=c("blue","red","purple","green","orange","black"),pch=c(2,8,1,20,6,5),lty=1,col=c("blue","red","purple","green","orange","black"),cex=0.8) dev.off() ###################################### setwd("/Users/cliu/Desktop/comparison/Seq-weights/final/3vs3") d=c(1,1.2,1.5,2,3,4,8) filenames=dir(pattern="final")[1:30] pdf("fig6box.pdf",height=8,width=8) par(mfrow=c(2,2)) group=c(rep(1:5,7)) names=c(rep(1.3,6),rep(1.5,6),rep(2,6),rep(3,6),rep(4,6)) sizes=rep(c(0.5,0.6,0.7,0.8,0.9,1),5) titles1=rep("NA",30) titles1[18]=c("(A)") titles1[30]=c("(B)") nq=c(18,30) for(l in nq){ if(l==18){ tmp1=filenames[l] load(tmp1) x1=c(seq(0.55,1.25,by=0.12),seq(1.55,2.25,by=0.12),seq(2.55,3.25,by=0.12),seq(3.55,4.25,by=0.12),seq(4.55,5.25,by=0.12),seq(5.55,6.25,by=0.12),seq(6.55,7.25,by=0.12)) bw=0.08 n=seq(1,42,by=6) boxplot(nd,boxwex=bw,at=x1[n],axes=F,xlim=c(0,8),ylim=c(0,200),col="red",ylab="Number of genes detected (FDR<0.1)",xlab="Sample variability") boxplot(ndpure,boxwex=bw,at=x1[n+1],col="blue",axes=F,add=T) boxplot(vnda2,boxwex=bw,at=x1[n+2],col="green",axes=F,add=T) boxplot(vndanew2,boxwex=bw,at=x1[n+3],col="orange",axes=F,add=T) axis(1, at=c(0,1,2,3,4,5,6,7,8), labels=c(NA,d,NA)) axis(2,at=c(0,50,100,150,200)) boxplot(ndvfremove,boxwex=bw,at=x1[n+4],col="gray",axes=F,add=T) boxplot(ndonly,boxwex=bw,at=x1[n+5],col="purple",axes=F,add=T) title(main=titles1[l]) legend("topright",legend=c("No Weights","Voom","Sample","Voom + Sample","Voom + Block","Sample Removal"),text.col=c("blue","red","purple","green","orange","grey"),cex=0.8) }else{ tmp1=filenames[l] load(tmp1) x1=c(seq(0.55,1.25,by=0.12),seq(1.55,2.25,by=0.12),seq(2.55,3.25,by=0.12),seq(3.55,4.25,by=0.12),seq(4.55,5.25,by=0.12),seq(5.55,6.25,by=0.12),seq(6.55,7.25,by=0.12)) bw=0.08 n=seq(1,42,by=6) boxplot(nd,boxwex=bw,at=x1[n],axes=F,xlim=c(0,8),ylim=c(100,250),col="red",ylab="Number of genes detected (FDR<0.1)",xlab="Sample variability") boxplot(ndpure,boxwex=bw,at=x1[n+1],col="blue",axes=F,add=T) boxplot(vnda2,boxwex=bw,at=x1[n+2],col="green",axes=F,add=T) boxplot(vndanew2,boxwex=bw,at=x1[n+3],col="orange",axes=F,add=T) axis(1, at=c(0,1,2,3,4,5,6,7,8), labels=c(NA,d,NA)) axis(2,at=c(100,150,200,250)) boxplot(ndvfremove,boxwex=bw,at=x1[n+4],col="gray",axes=F,add=T) boxplot(ndonly,boxwex=bw,at=x1[n+5],col="purple",axes=F,add=T) title(main=titles1[l]) } } titles1[18]=c("(C)") titles1[30]=c("(D)") for(l in nq){ if(l==18){ tmp=filenames[l] load(tmp) fd =fd/nsim fdpure=fdpure/nsim vfdanew2=vfdanew2/nsim vfda2=vfda2/nsim fdvfremove=fdvfremove/nsim fdonly=fdonly/nsim newfd=newfdpure=newvfdanew2=newvfda2=newfdvfremove=newfdonly=matrix(NA,100,7) for(j in 1:ncol(newfd)){ for(i in 1:100){ if(nd[i,j]==0){ newfd[i,j]=0 }else{ newfd[i,j]=fd[nd[i,j],j] } } } for(j in 1:ncol(newfd)){ for(i in 1:100){ if(ndpure[i,j]==0){ newfdpure[i,j]=0 }else{ newfdpure[i,j]=fdpure[ndpure[i,j],j] } } } for(j in 1:ncol(newfd)){ for(i in 1:100){ if(vndanew2[i,j]==0){ newvfdanew2[i,j]=0 }else{ newvfdanew2[i,j]=vfdanew2[vndanew2[i,j],j] } } } for(j in 1:ncol(newfd)){ for(i in 1:100){ if(vnda2[i,j]==0){ newvfda2[i,j]=0 }else{ newvfda2[i,j]=vfda2[vnda2[i,j],j] } } } for(j in 1:ncol(newfd)){ for(i in 1:100){ if(ndvfremove[i,j]==0){ newfdvfremove[i,j]=0 }else{ newfdvfremove[i,j]=fdvfremove[ndvfremove[i,j],j] } } } for(j in 1:ncol(newfd)){ for(i in 1:100){ if(ndonly[i,j]==0){ newfdonly[i,j]=0 }else{ newfdonly[i,j]=fdonly[ndonly[i,j],j] } } } newfd= newfd/nd newfd[is.infinite(newfd)]=NA newfdpure= newfdpure/ndpure newfdpure[is.infinite(newfdpure)]=NA newvfdanew2=newvfdanew2/vndanew2 newvfdanew2[is.infinite(newvfdanew2)]=NA newvfda2=newvfda2/vnda2 newvfda2[is.infinite(newvfda2)]=NA newfdvfremove= newfdvfremove/ndvfremove newfdvfremove[is.infinite(newfdvfremove)]=NA newfdonly=newfdonly/ndonly newfdonly[is.infinite(newfdonly)]=NA boxplot(newfd[,1:6],boxwex=bw,at=x1[n[1:6]],axes=F,xlim=c(0,8),ylim=c(0,0.2),col="red",ylab="Empirical FDR",xlab="Sample variability") boxplot(newfdpure[,1:4],boxwex=bw,at=x1[(n+1)[1:4]],col="blue",axes=F,add=T) boxplot(newvfda2,boxwex=bw,at=x1[n+2],col="green",axes=F,add=T) boxplot(newvfdanew2,boxwex=bw,at=x1[n+3],col="orange",axes=F,add=T) axis(1, at=c(0,1,2,3,4,5,6,7,8), labels=c(NA,d,NA)) axis(2,at=c(0,0.05,0.1,0.15,0.2)) boxplot(newfdvfremove,boxwex=bw,at=x1[n+4],col="gray",axes=F,add=T) boxplot(newfdonly[,1:5],boxwex=bw,at=x1[(n+5)[1:5]],col="purple",axes=F,add=T) title(main=titles1[l]) }else{ tmp=filenames[l] load(tmp) fd =fd/nsim fdpure=fdpure/nsim vfdanew2=vfdanew2/nsim vfda2=vfda2/nsim fdvfremove=fdvfremove/nsim fdonly=fdonly/nsim nd.fd=allnd=matrix(NA,6,7) newfd=newfdpure=newvfdanew2=newvfda2=newfdvfremove=newfdonly=matrix(NA,100,7) for(j in 1:ncol(newfd)){ for(i in 1:100){ newfd[i,j]=fd[ceiling(nd[i,j]+0.0000001),j] newfdpure[i,j]=fdpure[ceiling(ndpure[i,j]+0.000001),j] newvfdanew2[i,j]=vfdanew2[ceiling(vndanew2[i,j]+0.000001),j] newvfda2[i,j]=vfda2[ceiling(vnda2[i,j]+0.000001),j] newfdvfremove[i,j]=fdvfremove[ceiling(ndvfremove[i,j]+0.000001),j] newfdonly[i,j]=fdonly[ceiling(ndonly[i,j]+0.000001),j] } } newfd= newfd/nd newfdpure= newfdpure/ndpure newvfdanew2=newvfdanew2/vndanew2 newvfda2=newvfda2/vnda2 newfdvfremove= newfdvfremove/ndvfremove newfdonly=newfdonly/ndonly boxplot(newfd,boxwex=bw,at=x1[n],axes=F,xlim=c(0,8),ylim=c(0,0.2),col="red",ylab="Empirical FDR",xlab="Sample variability") boxplot(newfdpure,boxwex=bw,at=x1[(n+1)],col="blue",axes=F,add=T) boxplot(newvfda2,boxwex=bw,at=x1[n+2],col="green",axes=F,add=T) boxplot(newvfdanew2,boxwex=bw,at=x1[n+3],col="orange",axes=F,add=T) axis(1, at=c(0,1,2,3,4,5,6,7,8), labels=c(NA,d,NA)) axis(2,at=c(0,0.05,0.1,0.15,0.2)) boxplot(newfdvfremove,boxwex=bw,at=x1[n+4],col="gray",axes=F,add=T) boxplot(newfdonly,boxwex=bw,at=x1[(n+5)],col="purple",axes=F,add=T) title(main=titles1[l]) } } dev.off() ################################################################################################################ ################################################################################################################ ################################################################################################################ #null simulation plot # load("/Users/cliu/Desktop/comparison/Seq-weights/nullsimulation.rda") nd5 =nd5/ngenes/nsim #voom only# ndpure5=ndpure5/ngenes/nsim# pure model# vnda25=vnda25/ngenes/nsim #weight all (6 outlier)# vndanew25=vndanew25/ngenes/nsim #voom + weight by group# ndvfremove5=ndvfremove5/nsim/ngenes #remove# allnd=cbind(nd5,ndpure5,vnda25,vndanew25,ndvfremove5) pdf("NullSimulation.pdf") plot(c(1,7),c(0.002,0.012),type="n", ylab="Proportion of genes with p-value < 0.01",xlab="variability",main="null simulation",axes=F) axis(1, at=c(1,2,3,4,5,6,7), labels=d) axis(2) colssel=c("red","blue","green","orange","black") for(j in 1:5){ points(allnd[,j],col=colssel[j]) lines(allnd[,j],col=colssel[j],lwd=2) } abline(h=0.01,lty=2,col="grey") legend("bottomleft",legend=c("Voom-only","No-weights","Sample-weigts","Adjusted-Weights","Remove-outlier"),text.col=colssel) dev.off()