EBSeq/R/IsoSimu.R

   1 IsoSimu=function(DVDconstant=NULL, DVDqt1=NULL, DVDqt2=NULL, Conditions, NumofSample, NumofIso=NULL, DEIsoProp, Phiconstant=NULL, Phi.qt1=NULL, Phi.qt2=NULL,NormFactor=NULL, OnlyData=T)
   2 {
   3 # 2012 feb 1
   4 # paired simulation
   5 data(IsoEBresultGouldBart2)
   6 if(is.null(NormFactor)) NormFactor=rep(1,NumofSample)
   7
   8 MeansC1=IsoEBresultGouldBart2$C1Mean
   9 MeansC2=IsoEBresultGouldBart2$C2Mean
  10 MeanDVD=sapply(1:9,function(i) MeansC1[[i]]/MeansC2[[i]])
  11 # DVD library with each group here
  12 if (length(DVDconstant)==0) DVDLibrary= unlist(MeanDVD)[unlist(MeanDVD)<quantile(unlist(MeanDVD)[unlist(MeanDVD)!=Inf],DVDqt2) & unlist(MeanDVD)>quantile(unlist(MeanDVD)[unlist(MeanDVD)!=Inf],DVDqt1)]
  13
  14
  15
  16 # If DVD constant, use constant when generate
  17 # If not, use DVDLibrary
  18
  19 VarInput=IsoEBresultGouldBart2$VarList
  20 VarInputNg=list(VarInput[[1]],unlist(VarInput[c(2,4,6,8)]),unlist(VarInput[c(3,5,7,9)]))
  21 #If NumofIso=NULL, empirical # of Iso
  22 #If !=NULL , Input a 9-vector
  23 if(length(NumofIso)==0) NumofIso.raw=sapply(1:3,function(i)length(VarInputNg[[i]]))
  24 if(length(NumofIso)!=0) NumofIso.raw=NumofIso*2
  25
  26 PhiInput.raw=IsoEBresultGouldBart2$RList
  27 PhiInput.raw.Ng=list(PhiInput.raw[[1]],unlist(PhiInput.raw[c(2,4,6,8)]),unlist(PhiInput.raw[c(3,5,7,9)]))
  28
  29
  30 if (length(Phiconstant)==0){
  31         PhiLibrary=sapply(1:3,function(i)PhiInput.raw.Ng[[i]][1/PhiInput.raw.Ng[[i]]<quantile(1/PhiInput.raw.Ng[[i]],Phi.qt2) & 1/PhiInput.raw.Ng[[i]]>quantile(1/PhiInput.raw.Ng[[i]],Phi.qt1)],simplify=F)
  32         PhiIndex=sapply(1:3, function(i)sample(names(PhiLibrary[[i]]),NumofIso.raw[[i]],replace=T),simplify=F)
  33         PhiInputNg=sapply(1:3, function(i)PhiLibrary[[i]][PhiIndex[[i]]])
  34 }
  35 if (length(Phiconstant)!=0)PhiInputNg=sapply(1:3,function(i)rep(Phiconstant,NumofIso.raw[[i]]),simplify=F)
  36
  37 # Wanna DENumbers be proportion to 2
  38 DEIsoNumbers=round(NumofIso.raw*DEIsoProp/2)*2
  39 IsoNames=sapply(1:3,function(i)paste("I",i,c(1:NumofIso.raw[i]),sep="_"),simplify=F)
  40 MeanNg=list(IsoEBresultGouldBart2$MeanList[[1]],unlist(IsoEBresultGouldBart2$MeanList[c(2,4,6,8)]),
  41 unlist(IsoEBresultGouldBart2$MeanList[c(3,5,7,9)]))
  42 MeanInputNg=sapply(1:3, function(i)MeanNg[[i]][PhiIndex[[i]]])
  43
  44 for(i in 1:3){
  45         names(MeanInputNg[[i]])=IsoNames[[i]]
  46         names(PhiInputNg[[i]])=IsoNames[[i]]
  47         }
  48
  49 ##############################
  50 # Get Ng version to every one
  51 ##############################
  52
  53
  54 #########
  55 # data
  56 #########
  57 EEList=sapply(1:3,function(i) sapply(1:NumofIso.raw[[i]], function(j)sapply(1:NumofSample,function(h) rnbinom(1,mu=MeanInputNg[[i]][j]*NormFactor[h], size=PhiInputNg[[i]][j]))),simplify=F)
  58
  59
  60 generateDataraw=vector("list",3)
  61 MeanVector=vector("list",3)
  62 VarVector=vector("list",3)
  63 MOV.post=vector("list",3)
  64
  65
  66 for(g in 1:3){
  67     generateDataraw[[g]]=t(EEList[[g]][,1:NumofIso.raw[g]])
  68         if(length(DVDconstant)==0){
  69                 for(j in 1:NumofIso.raw[g]){
  70                  if (j<=(DEIsoNumbers[g]/2)) generateDataraw[[g]][j,((NumofSample/2)+1):NumofSample]=sapply((NumofSample/2+1):NumofSample, function(h)suppressWarnings(rnbinom(1, size=PhiInputNg[[g]][j], mu=sample(DVDLibrary,1)*MeanInputNg[[g]][j]*NormFactor[h])), simplify=T)
  71                 if (j>=((DEIsoNumbers[g]/2)+1) & j <=DEIsoNumbers[g]) generateDataraw[[g]][j,1:(NumofSample/2)]=sapply(1:(NumofSample/2),function(h) suppressWarnings(rnbinom(1, size=MeanInputNg[[g]][j], mu= sample(DVDLibrary,1)*MeanInputNg[[g]][j]*NormFactor[h])),simplify=T)
  72 }
  73          }
  74         if(length(DVDconstant)!=0){
  75         for(j in 1:NumofIso.raw[g]){
  76              if (j<=(DEIsoNumbers[g]/2)) generateDataraw[[g]][j,((NumofSample/2)+1):NumofSample]=sapply((NumofSample/2+1):NumofSample, function(h)suppressWarnings(rnbinom(1, DVDconstant*MeanInputNg[[g]][j]*NormFactor[h])),simplify=T)
  77              if (j>=((DEIsoNumbers[g]/2)+1) & j <=DEIsoNumbers[g]) generateDataraw[[g]][j,1:(NumofSample/2)]=sapply(1:(NumofSample/2),function(h) wuppressWarnings(rnbinom(1, DVDconstant*MeanInputNg[[g]][j]*NormFactor[h])),simplify=T)
  78                 }
  79         }
  80 rownames(generateDataraw[[g]])=IsoNames[[g]][1:NumofIso.raw[g]]
  81 MeanVector[[g]]=rowMeans(generateDataraw[[g]])
  82 VarVector[[g]]=apply(generateDataraw[[g]],1,var)
  83 MOV.post[[g]]=MeanVector[[g]]/VarVector[[g]]
  84 }
  85
  86
  87 ### Remove MOV=NA
  88 generateData=generateDataraw
  89 for (i in 1:3) generateData[[i]]=generateData[[i]][!is.na(MOV.post[[i]]),]
  90 #print(paste("NA MOV's",sum(is.na(unlist(MOV.post)))))
  91 NumDENow=sapply(1:3, function(i)sum(rownames(generateData[[i]])%in%rownames(generateDataraw[[i]])[1:DEIsoNumbers[i]]))
  92
  93 if(length(NumofIso)!=0){
  94             for(i in 1:3)
  95                 generateData[[i]]=generateData[[i]][c(sample(1:NumDENow[i],round(NumofIso[i]*DEIsoProp),replace=F),round( (dim(generateData[[i]])[1]+1-NumofIso[i]*(1-DEIsoProp)):dim(generateData[[i]])[1])),]
  96 }
  97 generateDataNg=generateData
  98
  99 ## DE
 100 UseName=sapply(1:3, function(i)rownames(generateData[[i]]),simplify=F)
 101 TrueDE=sapply(1:3, function(i)UseName[[i]][UseName[[i]] %in% rownames(generateDataraw[[i]])[1:DEIsoNumbers[i]]],simplify=F)
 102 TrueDE.unlist=do.call(c,TrueDE)
 103
 104 phiuse=sapply(1:3,function(i)PhiInputNg[[i]][UseName[[i]]])
 105 meanuse=sapply(1:3,function(i)MeanInputNg[[i]][UseName[[i]]])
 106
 107 #if(OnlyData==T){
 108
 109 OutName=sapply(1:3,function(i)paste("Iso",i,c(1:nrow(generateDataNg[[i]])),sep="_"))
 110 for(i in 1:3)names(OutName[[i]])=rownames(generateDataNg[[i]])
 111 OutData=generateDataNg
 112 for(i in 1:3)rownames(OutData[[i]])=as.vector(OutName[[i]])
 113 OutTrueDE=as.vector(unlist(OutName)[TrueDE.unlist])
 114 output=list(data=OutData, TrueDE=OutTrueDE)
 115
 116
 117 #output=list(data=generateDataNg, TrueDE=TrueDE.unlist)
 118 return(output)
 119 #    }
 120 # Now only OnlyData=T version
 121 }
 122