EBSeq/R/IsoSimuAt.R

   1 IsoSimuAt<-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 #Ng paired 2012 feb 1
   4 if(is.null(NormFactor)) NormFactor=rep(1,NumofSample)
   5 data(IsoEBresultGouldBart2)
   6
   7 MeansC1=IsoEBresultGouldBart2$C1Mean
   8 MeansC2=IsoEBresultGouldBart2$C2Mean
   9 MeanDVD=sapply(1:9,function(i) MeansC1[[i]]/MeansC2[[i]])
  10 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)]
  11
  12
  13
  14
  15 VarInput=IsoEBresultGouldBart2$VarList
  16 VarInputNg=list(VarInput[[1]],unlist(VarInput[c(2,4,6,8)]),unlist(VarInput[c(3,5,7,9)]))
  17
  18 if(length(NumofIso)==0) NumofIso=sapply(1:3,function(i)length(VarInputNg[[i]]))
  19 PhiInput.raw=IsoEBresultGouldBart2$RList
  20 PhiInput.raw.Ng=list(PhiInput.raw[[1]],unlist(PhiInput.raw[c(2,4,6,8)]),unlist(PhiInput.raw[c(3,5,7,9)]))
  21
  22
  23 if (length(Phiconstant)==0){
  24         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)
  25         PhiIndex=sapply(1:3, function(i)sample(names(PhiLibrary[[i]]),NumofIso[[i]],replace=T),simplify=F)
  26         PhiInputNg=sapply(1:3, function(i)PhiLibrary[[i]][PhiIndex[[i]]])
  27 }
  28 if (length(Phiconstant)!=0)PhiInputNg=sapply(1:3,function(i)rep(Phiconstant,NumofIso[[i]]),simplify=F)
  29
  30 # Wanna DENumbers be proportion to 2
  31 DEIsoNumbers=round(NumofIso*DEIsoProp/2)*2
  32 IsoNames=sapply(1:3,function(i)paste("I",i,c(1:NumofIso[i]),sep="_"),simplify=F)
  33 MeanNg=list(IsoEBresultGouldBart2$MeanList[[1]],unlist(IsoEBresultGouldBart2$MeanList[c(2,4,6,8)]),
  34 unlist(IsoEBresultGouldBart2$MeanList[c(3,5,7,9)]))
  35 MeanInputNg=sapply(1:3, function(i)MeanNg[[i]][PhiIndex[[i]]])
  36
  37 for(i in 1:3){
  38         names(MeanInputNg[[i]])=IsoNames[[i]]
  39         names(PhiInputNg[[i]])=IsoNames[[i]]
  40         }
  41
  42 #########
  43 # data
  44 #########
  45 EEList=sapply(1:3,function(i) sapply(1:NumofIso[[i]], function(j)sapply(1:NumofSample,function(h) rnbinom(1,mu=MeanInputNg[[i]][j]*NormFactor[h], size=PhiInputNg[[i]][j]))),simplify=F)
  46
  47
  48 generateDataraw=vector("list",3)
  49 MeanVector=vector("list",3)
  50 VarVector=vector("list",3)
  51 MOV.post=vector("list",3)
  52
  53
  54 for(g in 1:3){
  55     generateDataraw[[g]]=t(EEList[[g]][,1:NumofIso[g]])
  56         if(length(DVDconstant)==0){
  57                 for(j in 1:NumofIso[g]){
  58                  if (j<=(DEIsoNumbers[g]/2)) generateDataraw[[g]][j,((NumofSample/2)+1):NumofSample]=sapply((NumofSample/2+1):NumofSample, function(h)rnbinom(1, size=PhiInputNg[[g]][j], mu=sample(DVDLibrary,1)*MeanInputNg[[g]][j]*NormFactor[h]), simplify=T)
  59                 if (j>=((DEIsoNumbers[g]/2)+1) & j <=DEIsoNumbers[g]) generateDataraw[[g]][j,1:(NumofSample/2)]=sapply(1:(NumofSample/2),function(h) rnbinom(1, size=MeanInputNg[[g]][j], mu= sample(DVDLibrary,1)*MeanInputNg[[g]][j]*NormFactor[h]),simplify=T)
  60 }
  61          }
  62         if(length(DVDconstant)!=0){
  63         for(j in 1:NumofIso[g]){
  64              if (j<=(DEIsoNumbers[g]/2)) generateDataraw[[g]][j,((NumofSample/2)+1):NumofSample]=sapply((NumofSample/2+1):NumofSample, function(h)rnbinom(1, DVDconstant*MeanInputNg[[g]][j]*NormFactor[h]),simplify=T)
  65              if (j>=((DEIsoNumbers[g]/2)+1) & j <=DEIsoNumbers[g]) generateDataraw[[g]][j,1:(NumofSample/2)]=sapply(1:(NumofSample/2),function(h) rnbinom(1, DVDconstant*MeanInputNg[[g]][j]*NormFactor[h]),simplify=T)
  66                 }
  67         }
  68 rownames(generateDataraw[[g]])=IsoNames[[g]][1:NumofIso[g]]
  69 MeanVector[[g]]=rowMeans(generateDataraw[[g]])
  70 VarVector[[g]]=apply(generateDataraw[[g]],1,var)
  71 MOV.post[[g]]=MeanVector[[g]]/VarVector[[g]]
  72 }
  73
  74
  75 ### Remove MOV=NA
  76 generateData=generateDataraw
  77 for (i in 1:3) generateData[[i]]=generateData[[i]][!is.na(MOV.post[[i]]),]
  78 print(paste("NA MOV's",sum(is.na(unlist(MOV.post)))))
  79 #tmpmean=sapply(1:9,function(i)rowMeans(generateData[[i]]))
  80 #tmpvar=sapply(1:9,function(i)apply(generateData[[i]],1,var))
  81 #source("plot_functions.R")
  82 #CheckSimuNg(tmpmean,tmpvar,c(-1,5),c(-1,7))
  83
  84
  85
  86
  87 ## DE
  88 UseName=sapply(1:3, function(i)rownames(generateData[[i]]),simplify=F)
  89 TrueDE=sapply(1:3, function(i)UseName[[i]][UseName[[i]] %in% rownames(generateData[[i]])[1:DEIsoNumbers[i]]],simplify=F)
  90 TrueDE.unlist=do.call(c,TrueDE)
  91
  92 TrueDELength=sapply(TrueDE,length)
  93
  94 AtNames_Level=vector("list",4)
  95 AtLoc=vector("list",3)
  96 AtFold=vector("list",3)
  97 names(AtNames_Level)=c(4,6,8,10)
  98
  99
 100 for(j in 1:3){
 101 AtLoc[[j]]=sample(c(1:length(Conditions)), TrueDELength[j], replace=T)
 102 AtFold[[j]]=sample(c(4,6,8,10),TrueDELength[j], replace=T)
 103
 104 for(i in 1:TrueDELength[j]){
 105
 106 generateData[[j]][(TrueDELength[j]+i),AtLoc[[j]][i]]=generateData[[j]][(TrueDELength[j]+i),AtLoc[[j]][i]]*AtFold[[j]][i]
 107 AtNames_Level[[as.character(AtFold[[j]][i])]]=c(AtNames_Level[[as.character(AtFold[[j]][i])]],rownames(generateData[[j]])[TrueDELength[j]+i])
 108 }
 109 }
 110 phiuse=sapply(1:3,function(i)PhiInputNg[[i]][UseName[[i]]])
 111 meanuse=sapply(1:3,function(i)MeanInputNg[[i]][UseName[[i]]])
 112
 113 #generateDataNg=list(generateData[[1]], do.call(rbind,generateData[c(2,4,6,8)]), do.call(rbind,generateData[c(3,5,7,9)]))
 114 generateDataNg=generateData
 115
 116 #if(OnlyData==T){
 117
 118 OutName=sapply(1:3,function(i)paste("Iso",i,c(1:nrow(generateDataNg[[i]])),sep="_"))
 119 for(i in 1:3)names(OutName[[i]])=rownames(generateDataNg[[i]])
 120 OutData=generateDataNg
 121 for(i in 1:3)rownames(OutData[[i]])=as.vector(OutName[[i]])
 122 OutTrueDE=as.vector(unlist(OutName)[TrueDE.unlist])
 123 OutAt=as.vector(unlist(OutName)[AtNames <- Level])
 124
 125 output=list(data=OutData, TrueDE=OutTrueDE, Outliers=OutAt)
 126 #       return(output)
 127 #    }
 128         }