EBSeq/demo/EBSeq.R

   1 library(EBSeq)
   2 set.seed(13)
   3
   4 # Section 3.1
   5
   6 GeneGenerate=GeneSimu(DVDconstant=4, DVDqt1=NULL, DVDqt2=NULL,
   7   Conditions=rep(c(1,2),each=5), NumofSample=10, NumofGene=10000,
   8   DEGeneProp=.1, Phiconstant=NULL, Phi.qt1=.1, Phi.qt2=.9,
   9   Meanconstant=NULL, OnlyData=T)
  10 GeneData=GeneGenerate$data
  11 GeneTrueDENames=GeneGenerate$TrueDE
  12 str(GeneData)
  13 str(GeneTrueDENames)
  14
  15 Sizes=MedianNorm(GeneData)
  16
  17 EBres=EBTest(Data=GeneData,
  18   Conditions=as.factor(rep(c(1,2),each=5)),sizeFactors=Sizes, maxround=5)
  19
  20 PP=GetPP(EBres)
  21 str(PP)
  22 DEfound=names(PP)[which(PP>=.95)]
  23 str(DEfound)
  24 sum(DEfound%in%GeneTrueDENames)
  25
  26 QQP(QList=EBres$QList1, AlphaResult=EBres[[1]][5,1],
  27   BetaResult=EBres[[2]][5,1], name="Gene Simulation", AList="F", GroupName=NULL)
  28 DenNHist(QList=EBres$QList1, Alpha=EBres[[1]][5,1], Beta=EBres[[2]][5,1],
  29   name="Gene Simulation", AList="F", GroupName=NULL)
  30
  31 # Section 3.2
  32
  33 IsoGenerate=IsoSimu(DVDconstant=NULL, DVDqt1=.97, DVDqt2=.98,
  34   Conditions=as.factor(rep(c(1,2),each=5)), NumofSample=10,
  35   NumofIso=c(1000,2000,3000), DEIsoProp=.1, Phiconstant=NULL,
  36   Phi.qt1=.25, Phi.qt2=.75, OnlyData=T )
  37 str(IsoGenerate)
  38
  39 IsoMat=do.call(rbind,IsoGenerate$data)
  40 str(IsoMat)
  41
  42 IsoSizes=MedianNorm(IsoMat)
  43
  44 IsoNames=rownames(IsoMat)
  45 str(IsoNames)
  46 GeneNames=paste("Gene",c(1:3000),sep="_")
  47 IsosGeneNames=c(GeneNames[1:1000],rep(GeneNames[1001:2000],each=2),
  48   rep(GeneNames[2001:3000],each=3))
  49 NgList=GetNg(IsoNames, IsosGeneNames)
  50 IsoNgTrun=NgList$IsoformNgTrun
  51 IsoNgTrun[c(1:3,1001:1003,3001:3003)]
  52
  53 IsoEBres=EBTest(Data=IsoMat, NgVector=IsoNgTrun,
  54   Conditions=as.factor(rep(c(1,2),each=5)),sizeFactors=IsoSizes, maxround=5)
  55 IsoPP=GetPP(IsoEBres)
  56 str(IsoPP)
  57 IsoDE=IsoPP[which(IsoPP>=.95)]
  58 str(IsoDE)
  59 sum(names(IsoDE)%in%IsoGenerate$TrueDE)
  60
  61 par(mfrow=c(2,2))
  62 PolyFitValue=vector("list",3)
  63 for(i in 1:3)
  64   PolyFitValue[[i]]=PolyFitPlot(IsoEBres$C1Mean[[i]],
  65     IsoEBres$C1EstVar[[i]],5)
  66
  67 PolyAll=PolyFitPlot(unlist(IsoEBres$C1Mean), unlist(IsoEBres$C1EstVar),5)
  68 lines(log10(IsoEBres$C1Mean[[1]][PolyFitValue[[1]]$sort]),
  69   PolyFitValue[[1]]$fit[PolyFitValue[[1]]$sort],col="yellow")
  70 lines(log10(IsoEBres$C1Mean[[2]][PolyFitValue[[2]]$sort]),
  71   PolyFitValue[[2]]$fit[PolyFitValue[[2]]$sort],col="pink")
  72 lines(log10(IsoEBres$C1Mean[[3]][PolyFitValue[[3]]$sort]),
  73   PolyFitValue[[3]]$fit[PolyFitValue[[3]]$sort],col="green")
  74 legend("topleft",c("All Isoforms","Ng = 1","Ng = 2","Ng = 3"),
  75   col=c("red","yellow","pink","green"),lty=1,lwd=3,box.lwd=2)
  76
  77 par(mfrow=c(2,2))
  78 QQP(QList=IsoEBres$QList1, AlphaResult=IsoEBres[[1]][5,],
  79  BetaResult=IsoEBres[[2]][5,],
  80  name="Isoforms", AList="F", GroupName=paste("Ng = ",c(1:3),sep=""))
  81
  82 DenNHist(QList=IsoEBres$QList1, Alpha=IsoEBres[[1]][5,],
  83   Beta=IsoEBres[[2]][5,],
  84   name="Isoforms", AList="F", GroupName=paste("Ng = ",c(1:3),sep=""))
  85
  86 # Section 3.3
  87
  88 Conditions=c("C1","C1","C2","C2","C3","C3")
  89 PosParti=GetPatterns(Conditions)
  90 PosParti
  91
  92 Parti=PosParti[-3,]
  93 Parti
  94
  95 MultiData=GeneMultiSimu(Conditions=Conditions,AllParti=Parti,
  96           NumofSample=6,NumofGene=1000,DEGeneProp=c(.7,.1,.1,.1),
  97           DVDqt1=.98,DVDqt2=.99,Phi.qt1=.25,Phi.qt2=.75)
  98 str(MultiData)
  99
 100 MultiSize=MedianNorm(MultiData$data)
 101 MultiRes=EBMultiTest(MultiData$data,NgVector=NULL,Conditions=Conditions,
 102            AllParti=Parti, sizeFactors=MultiSize, maxround=5)
 103 MultiPP=GetMultiPP(MultiRes)
 104 names(MultiPP)
 105 MultiPP$PP[1:10,]
 106 MultiPP$MAP[1:10]
 107 MultiPP$Patterns
 108 sum(MultiPP$MAP==MultiData$Patterns)
 109
 110 # EOF