EBSeq/R/MergeGene.R

   1 MergeGene <-
   2 function(GeneSIMout, Num, Path="./"){
   3 NumSample=ncol(GeneSIMout[[i]]$generateData)
   4
   5 NumGene=rep(0,Num)
   6 for (i in 1:Num)NumGene[i]=nrow(GeneSIMout[[i]]$generateData)
   7
   8 MinNumGene=min(NumGene)
   9 AproxNumDE=length(GeneSIMout[[1]]$TrueDE)
  10
  11 GeneMergeTable=matrix(rep(0,12),nrow=6)
  12         for(i in 1:Num)GeneMergeTable=GeneMergeTable+GeneSIMout[[i]][[1]]
  13         GeneMergeTable=GeneMergeTable/Num
  14         GeneMergeTable=round(GeneMergeTable,2)
  15
  16         GeneMergeDVD=rep(0,2)
  17           for(i in 1:Num)GeneMergeDVD=GeneMergeDVD+GeneSIMout[[i]][[3]]
  18                   GeneMergeDVD=round(GeneMergeDVD/Num,2)
  19
  20           GeneMergePhi=matrix(rep(0,2),nrow=2)
  21                   for(i in 1:Num)GeneMergePhi=GeneMergePhi+GeneSIMout[[i]][[4]]
  22                           GeneMergePhi=round(GeneMergePhi/Num,2)
  23 ## Write
  24 TXTname=paste(paste(Path,paste("Gene","DVD",GeneMergeDVD[1], GeneMergeDVD[2],"Phi",GeneMergePhi[1], GeneMergePhi[2],"Sample",NumSample,sep="_"),sep=""),".txt",sep="")
  25 write.table(GeneMergeTable, file=TXTname)
  26
  27
  28 ####### Note everytime # DE genes and # total genes may different. (since NA issue)
  29   GeneMergeFD=matrix(rep(0,5*MinNumGene),ncol=5)
  30   GeneMergeFD.p=matrix(rep(0,5*MinNumGene),ncol=5)
  31   GeneMergeTP.p=matrix(rep(0,5*MinNumGene),ncol=5)
  32   GeneMergeFN.p=matrix(rep(0,5*MinNumGene),ncol=5)
  33   GeneMergeTN.p=matrix(rep(0,5*MinNumGene),ncol=5)
  34
  35   GeneMergeFDR=matrix(rep(0,5*MinNumGene),ncol=5)
  36   GeneMergeTPR=matrix(rep(0,5*MinNumGene),ncol=5)
  37   GeneMergeFPR=matrix(rep(0,5*MinNumGene),ncol=5)
  38
  39
  40   for(i in 1:Num){
  41         # Make sure names in the same order
  42         # Get FD number for each number of genes found
  43     TotalNum=nrow(GeneSIMout[[i]]$generateData)
  44         NumDE=length(GeneSIMout[[i]]$TrueDE)
  45         EBSeqNames=names(GeneSIMout[[i]]$EBSeqPP)
  46     tmpMatrix=cbind(GeneSIMout[[i]]$DESeqP[EBSeqNames],GeneSIMout[[i]]$edgeRP[EBSeqNames], exp(GeneSIMout[[i]]$BaySeqPP[EBSeqNames,2]),GeneSIMout[[i]]$BBSeqP[EBSeqNames],GeneSIMout[[i]]$EBSeqPP)
  47         # Bayseq and EBseq are PP. Others are p value
  48     tmpFD=TopCts(tmpMatrix, c(0,0,1,0,1), GeneSIMout[[i]]$TrueDE[GeneSIMout[[i]]$TrueDE%in%EBSeqNames], MinNumGene)
  49     # Get percentage for FP, TP, TN, FN!
  50         tmpFD.p=tmpFD/TotalNum
  51         # TP = Find - FD
  52         tmpTP.p=(c(1:MinNumGene)-tmpFD)/TotalNum
  53         # FN = TrueDE - TP
  54         tmpFN.p=NumDE/TotalNum - tmpTP.p
  55         # TN = TrueEE - FD
  56         tmpTN.p=(TotalNum-NumDE)/TotalNum - tmpFD.p
  57
  58
  59         tmpFDR=tmpFD.p/(tmpFD.p+tmpTP.p)
  60         tmpFPR=tmpFD.p/(tmpFD.p+tmpTN.p)
  61         tmpTPR=tmpTP.p/(tmpFN.p+tmpTP.p)
  62         GeneMergeFDR=GeneMergeFDR+tmpFDR
  63         GeneMergeTPR=GeneMergeTPR+tmpTPR
  64         GeneMergeFPR=GeneMergeFPR+tmpFPR
  65
  66     GeneMergeFD.p=GeneMergeFD.p+tmpFD.p
  67         GeneMergeTP.p=GeneMergeTP.p+tmpTP.p
  68         GeneMergeFN.p=GeneMergeFN.p+tmpFN.p
  69         GeneMergeTN.p=GeneMergeTN.p+tmpTN.p
  70
  71         GeneMergeFD=GeneMergeFD+tmpFD
  72  }
  73   GeneMergeFD=GeneMergeFD/Num
  74   GeneMergeFD.p=GeneMergeFD.p/Num
  75   GeneMergeTP.p=GeneMergeTP.p/Num
  76   GeneMergeFN.p=GeneMergeFN.p/Num
  77   GeneMergeTN.p=GeneMergeTN.p/Num
  78
  79   GeneMergeFDR=GeneMergeFDR/Num
  80   GeneMergeTPR=GeneMergeTPR/Num
  81   GeneMergeFPR=GeneMergeFPR/Num
  82
  83
  84 PlotTopName=paste(paste(Path,paste("Top","Gene","DVD",GeneMergeDVD[1], GeneMergeDVD[2],"Phi",GeneMergePhi[1], GeneMergePhi[2],"Sample",NumSample, sep="_"),sep=""),".pdf",sep="")
  85
  86 TrueDELength=length(GeneSIMout[[i]]$TrueDE[GeneSIMout[[i]]$TrueDE%in%EBSeqNames])
  87 pdf(PlotTopName)
  88   PlotTopCts(TrueDELength,GeneMergeFD[1:TrueDELength,],c("DESeq","edgeR","BaySeq","BBSeq","EBSeq"))
  89 dev.off()
  90
  91
  92 PlotFDName=paste(paste(Path,paste("FDTP","Gene","DVD",GeneMergeDVD[1], GeneMergeDVD[2],"Phi",GeneMergePhi[1], GeneMergePhi[2],"Sample",NumSample,sep="_"),sep=""),".pdf",sep="")
  93 pdf(PlotFDName)
  94   PlotFDTP(MinNumGene,GeneMergeFDR, GeneMergeTPR, c("DESeq","edgeR","BaySeq","BBSeq","EBSeq"))
  95 dev.off()
  96
  97 PlotFPName=paste(paste(Path,paste("FPRTP","Gene","DVD",GeneMergeDVD[1], GeneMergeDVD[2],"Phi",GeneMergePhi[1], GeneMergePhi[2],"Sample",NumSample,sep="_"),sep=""),".pdf",sep="")
  98 pdf(PlotFPName)
  99   PlotFPTP(MinNumGene,GeneMergeFPR, GeneMergeTPR, c("DESeq","edgeR","BaySeq","BBSeq","EBSeq"))
 100   dev.off()
 101
 102
 103 out=list(GeneMergeTable=GeneMergeTable, GeneMergeDVD=GeneMergeDVD, GeneMergePhi=GeneMergePhi, GeneMergeFD=GeneMergeFD)
 104
 105
 106 }
 107