EBSeq/R/TPFDRplot.R

   1 TPFDRplot <-
   2 function(DESeqP, EBZ, TrueDE, main, FDR=NULL){
   3         Seq=seq(0.001,0.5,by=0.001)
   4         DETPR=rep(0,length(Seq))
   5         EBTPR=rep(0,length(Seq))
   6         DEFDR=rep(0,length(Seq))
   7         EBFDR=rep(0,length(Seq))
   8         DETPNum=rep(0,length(Seq))
   9     EBTPNum=rep(0,length(Seq))
  10     DEFDNum=rep(0,length(Seq))
  11     EBFDNum=rep(0,length(Seq))
  12         for (i in 1:length(Seq)){
  13                 DESeqOnes=names(DESeqP)[DESeqP<=Seq[i]]
  14                 if (length(FDR)==0) EBOnes=names(EBZ)[EBZ>=crit.fun(1-EBZ, Seq[i])]
  15                 else if (FDR=="H") EBOnes=names(EBZ)[EBZ>=(1-Seq[i])]
  16                         else EBOnes=names(EBZ)[EBZ>=FDR[i]]
  17
  18                 DETPNum[i]=sum(DESeqOnes%in%TrueDE)
  19                 EBTPNum[i]=sum(EBOnes%in%TrueDE)
  20                 DEFDNum[i]=sum(!DESeqOnes%in%TrueDE)
  21                 EBFDNum[i]=sum(!EBOnes%in%TrueDE)
  22
  23                 DETPR[i]=DETPNum[i]/length(TrueDE)
  24                 EBTPR[i]=EBTPNum[i]/length(TrueDE)
  25                 DEFDR[i]=DEFDNum[i]/length(TrueDE)
  26                 EBFDR[i]=EBFDNum[i]/length(TrueDE)
  27         }
  28         plot(Seq,DETPR,ylim=c(0,1),xlim=c(0,.5),type="l",col="red", main=paste(main, "TPR"),xlab="controled FDR level", ylab="TPR",lwd=2)
  29         lines(Seq,EBTPR,col="blue",lwd=2)
  30         legend("bottomright",lwd=2, col=c("red","blue"), c("DESeq","EBSeq"))
  31
  32         plot(Seq,DEFDR,ylim=c(0,1),xlim=c(0,.5),type="l",col="red", main=paste(main, "FDR"),xlab="controled FDR level", ylab="TPR",lwd=2)
  33         lines(Seq,EBFDR,col="blue",lwd=2)
  34         legend("topleft", lwd=2, col=c("red","blue"), c("DESeq","EBSeq"))
  35
  36
  37         output=cbind( DETPR,EBTPR, DEFDR,EBFDR,DETPNum,EBTPNum,DEFDNum,EBFDNum)
  38 }
  39