R/LogNMulti.R

In wiscstatman/EBSeq: An R package for gene and isoform differential expression analysis of RNA-seq data

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+600)
		rownames(FMat)=rownames(FPartiLog)=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)


		if(length(zGood)==0){
		#Min=min(min(F0Log[which(F0Log!=-Inf)]), 
		#	min(F1Log[which(F1Log!=-Inf)]))
		tmpMat=FPartiLog
		tmpMean=apply(tmpMat,1,mean)
		FLogMdf=FPartiLog-tmpMean
		FMdf=exp(FLogMdf)

		FmultiPMdf=FMdf*PInMat
		DenomMdf=rowSums(FmultiPMdf)
		ZEach=apply(FmultiPMdf,2,function(i)i/DenomMdf)
		zNaNName1Mdf=names(DenomMdf)[is.na(DenomMdf)]
		# other NAs in LikeFun
		LFMdf=ZEach*(log(FmultiPMdf))
		zNaNMoreMdf=rownames(LFMdf)[which(is.na(rowSums(LFMdf)))]
		zNaNNameMdf=unique(c(zNaNName1Mdf,zNaNMoreMdf))
		zGood=which(!rownames(LFMdf)%in%zNaNNameMdf)



	
		}

		ZEachGood=ZEach[zGood,]
		###Update P
        PFromZ=colSums(ZEach[zGood,])/length(zGood)
        FGood=FPartiLog[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
    }