Included EBSeq for downstream differential expression analysis

[rsem.git] / EBSeq / R / LogNMultiEMP.R

diff --git a/EBSeq/R/LogNMultiEMP.R b/EBSeq/R/LogNMultiEMP.R
new file mode 100644 (file)
index 0000000..f852540
--- /dev/null
+++ b/EBSeq/R/LogNMultiEMP.R
@@ -0,0 +1,54 @@
+LogNMulti <-
+function(Input, InputSP, EmpiricalR, EmpiricalRSP, NumOfEachGroup, AlphaIn, BetaIn,  PIn, NoneZeroLength, AllParti, Conditions)
+{
+
+        #For each gene (m rows of Input---m genes)
+        #Save each gene's F0, F1 for further likelihood calculation. 
+               FList=sapply(1:nrow(AllParti),function(i)sapply(1:nlevels(as.factor(AllParti[i,])),
+                                  function(j)f0(do.call(cbind,InputSP[AllParti[i,]==j]),AlphaIn, BetaIn, 
+                                                                        do.call(cbind,EmpiricalRSP[AllParti[i,]==j]), NumOfEachGroup, log=T)),
+                                         simplify=F) 
+               FPartiLog=sapply(FList,rowSums)
+               FMat=exp(FPartiLog)
+               rownames(FMat)=rownames(Input)
+        #Get z
+               #Use data.list in logfunction
+        PInMat=matrix(rep(1,nrow(Input)),ncol=1)%*%matrix(PIn,nrow=1)
+               FmultiP=FMat*PInMat
+               Denom=rowSums(FmultiP)
+               ZEach=apply(FmultiP,2,function(i)i/Denom)
+               zNaNName1=names(Denom)[is.na(Denom)]
+               # other NAs in LikeFun
+               LF=ZEach*(log(FmultiP))
+               zNaNMore=rownames(LF)[which(is.na(rowSums(LF)))]
+               zNaNName=unique(c(zNaNName1,zNaNMore))
+               zGood=which(!rownames(LF)%in%zNaNName)
+               ZEachGood=ZEach[zGood,]
+               ###Update P
+        PFromZ=colSums(ZEach[zGood,])/length(zGood)
+        FGood=FMat[zGood,]
+               ### MLE Part ####
+        # Since we dont wanna update p and Z in this step
+        # Each Ng for one row
+               
+               NumGroupVector=rep(c(1:NoneZeroLength),NumOfEachGroup)
+               
+               NumGroupVector.zGood=NumGroupVector[zGood]
+               NumOfEachGroup.zGood=tapply(NumGroupVector.zGood,NumGroupVector.zGood,length)
+
+        StartValue=c(AlphaIn, BetaIn,PIn[-1])
+               InputSPGood=sapply(1:length(InputSP),function(i)InputSP[[i]][zGood,],simplify=F)
+        EmpiricalRSPGood=sapply(1:length(EmpiricalRSP),function(i)EmpiricalRSP[[i]][zGood,],simplify=F)
+
+               Result<-optim(StartValue,LikefunMulti,InputPool=list(InputSPGood,Input[zGood,],ZEach[zGood,], 
+                                        NoneZeroLength,EmpiricalR[zGood, ],EmpiricalRSPGood, NumOfEachGroup.zGood, AllParti))
+               AlphaNew= Result$par[1]
+               BetaNew=Result$par[2:(1+NoneZeroLength)]
+        PNewNo1=Result$par[(2+NoneZeroLength):length(Result$par)]
+               PNew=c(1-sum(PNewNo1),PNewNo1)
+               ##
+        Output=list(AlphaNew=AlphaNew,BetaNew=BetaNew,PNew=PNew,ZEachNew=ZEach, ZEachGood=ZEachGood, 
+                                       PFromZ=PFromZ, zGood=zGood, zNaNName=zNaNName,FGood=FGood)
+        Output
+    }
+