EBSeq/R/MergeIso.R

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