Nothing
R/iASeq.R
In iASeq: iASeq: integrating multiple sequencing datasets for detecting allele-specific events
Defines functions
singleEMfit
plotBIC
plotMotif
iASeqmotif
ASErawfit
iASeqmotiffit
fdown.loglike
fup.loglike
f0.loglike
Documented in
ASErawfit
f0.loglike
fdown.loglike
fup.loglike
iASeqmotif
iASeqmotiffit
plotBIC
plotMotif
singleEMfit
f0.loglike<-function(xj,nj,cj,dj)
{
lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
fup.loglike<-function(xj,nj,cj,dj,pj)
{
log(1/(1-pj))+log(1-pbeta(pj,xj+cj,nj-xj+dj))+lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
fdown.loglike<-function(xj,nj,cj,dj,pj)
{
log(1/pj)+log(pbeta(pj,xj+cj,nj-xj+dj))+lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
iASeqmotiffit<-function(exprs,studyid,repid,refid,K,iter.max=100,tol=1e-3)
{
set.seed(1)
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
########################################################################
#parameter preparation
########################################################################
pai=rep(0,K+1)
p0j=rep(0,ncol(xij))
p00j=rep(0,ncol(xij))
#need fitting later
c0j<-rep(1,ncol(xij))
d0j<-rep(1,ncol(xij))
cd<-rep(1,studynum)
dd<-rep(1,studynum)
#initialization
pai<-c(0.8,rep(0.2/K,K))
q0<-runif(K*studynum)
q1<-runif(K*studynum)
q2<-runif(K*studynum)
qsum=(q0+q1+q2)
qup=matrix(q1/qsum,K,studynum)
qdown=matrix(q2/qsum,K,studynum)
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
prob=(tempx)/(tempn)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
p00j[sampleid]=mean(prob[which(tempn!=0)])
}
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
addterm=p00j[sampleid]
prob=(tempx+addterm*2)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0j[sampleid]=mu*(mu*(1-mu)/v-1)
d0j[sampleid]=(1-mu)*(mu*(1-mu)/v-1)
p0j[sampleid]=mean(prob[which(tempn!=0)])
}
}
c0j=abs(c0j)
d0j=abs(d0j)
loglike.old<--1e10
########################################################################
#model fitting
########################################################################
loglike0<-list()
loglikeup<-list()
loglikedown<-list()
for(i.iter in 1:iter.max)
{
if(i.iter%%50==0)
{
print(paste("We have run ",i.iter," iterations for K=",K, sep=""))
}
err<-tol+1
for(d in 1:studynum)
{
temploglike0=rep(0,sitenum)
temploglikeup=rep(0,sitenum)
temploglikedown=rep(0,sitenum)
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
sampleid=repseq[j]
temploglike0<-temploglike0+f0.loglike(xij[,sampleid],nij[,sampleid],c0j[sampleid],d0j[sampleid])
temploglikeup<-temploglikeup+fup.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
temploglikedown<-temploglikedown+fdown.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
}
loglike0[[d]]<-temploglike0
loglikeup[[d]]<-temploglikeup
loglikedown[[d]]<-temploglikedown
}
clustlike<-matrix(0,sitenum,K+1)
condlikeup<-list()
condlikedown<-list()
for(d in 1:studynum)
{
condlikeup[[d]]<-matrix(0,sitenum,K+1)
condlikedown[[d]]<-matrix(0,sitenum,K+1)
}
#for ai=0 case
for(d in 1:studynum)
{
clustlike[,1]=clustlike[,1]+loglike0[[d]]
}
clustlike[,1]=clustlike[,1]+log(pai[1])
#for ai!=0 case
templike<-matrix(0,sitenum,3)
for(k in 2:(K+1))
{
for(d in 1:studynum)
{
templike[,1]<-log(qup[k-1,d])+loglikeup[[d]]
templike[,2]<-log(1-qup[k-1,d]-qdown[k-1,d])+loglike0[[d]]
templike[,3]<-log(qdown[k-1,d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
clustlike[,k]<-clustlike[,k]+tempmax+log(tempsum)
condlikeup[[d]][,k]<-templike[,1]/tempsum
condlikedown[[d]][,k]<-templike[,3]/tempsum
}
clustlike[,k]<-clustlike[,k]+log(pai[k])
}
tempmax<-apply(clustlike,1,max)
for(k in 1:(K+1))
{
clustlike[,k]<-exp(clustlike[,k]-tempmax)
}
tempsum<-apply(clustlike,1,sum)
#update motif occurrence rate
for(k in 1:(K+1))
{
clustlike[,k]<-clustlike[,k]/tempsum
}
pai.new<-(apply(clustlike,2,sum)+1)/(sitenum+K+1)
#update motifs
qup.new<-matrix(0,K,studynum)
qdown.new<-matrix(0,K,studynum)
for(k in 1:K)
{
clustpsum<-sum(clustlike[,k+1])
for(d in 1:studynum)
{
qup.new[k,d]<-(sum(clustlike[,k+1]*condlikeup[[d]][,k+1])+1)/(clustpsum+3)
qdown.new[k,d]<-(sum(clustlike[,k+1]*condlikedown[[d]][,k+1])+1)/(clustpsum+3)
}
}
loglike.new<-(sum(tempmax+log(tempsum))+sum(log(pai.new))
+sum(log(qup.new)+log(qdown.new)+log(1-qup.new-qdown.new)))/sitenum
## evaluate convergence
err.pai<-max(abs(pai.new-pai)/pai)
err.qup<-max(abs(qup.new-qup)/qup)
err.qdown<-max(abs(qdown.new-qdown)/qdown)
err<-max(err.pai, err.qup,err.qdown)
pai<-pai.new
qup<-qup.new
qdown<-qdown.new
loglike.old<-loglike.new
if(err<tol)
{
break
}
}
for(d in 1:studynum)
{
condlikeup[[d]]<-matrix(0,sitenum,K+1)
condlikedown[[d]]<-matrix(0,sitenum,K+1)
}
#for ai=0 case
for(d in 1:studynum)
{
clustlike[,1]=clustlike[,1]+loglike0[[d]]
}
clustlike[,1]=clustlike[,1]+log(pai[1])
#for ai!=0 case
templike<-matrix(0,sitenum,3)
for(k in 2:(K+1))
{
for(d in 1:studynum)
{
templike[,1]<-log(qup[k-1,d])+loglikeup[[d]]
templike[,2]<-log(1-qup[k-1,d]-qdown[k-1,d])+loglike0[[d]]
templike[,3]<-log(qdown[k-1,d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
clustlike[,k]<-clustlike[,k]+tempmax+log(tempsum)
condlikeup[[d]][,k]<-templike[,1]/tempsum
condlikedown[[d]][,k]<-templike[,3]/tempsum
}
clustlike[,k]<-clustlike[,k]+log(pai[k])
}
tempmax<-apply(clustlike,1,max)
for(k in 1:(K+1))
{
clustlike[,k]<-exp(clustlike[,k]-tempmax)
}
tempsum<-apply(clustlike,1,sum)
for(k in 1:(K+1))
{
clustlike[,k]<-clustlike[,k]/tempsum
}
p.post=1-(clustlike[,1])
minus=sum(log(pai))+sum(log(qup)+log(qdown)+log(1-qup-qdown))
result<-list(p.post=p.post, motif.prior=pai,motif.qup=qup,motif.qdown=qdown,clustlike=clustlike,
c0j=c0j,d0j=d0j,loglike=loglike.new*sitenum-minus)
}
ASErawfit<-function(exprs,studyid,repid,refid)
{
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
#result<-list(xij=xij,nij=nij,design=experiment.unique)
########################################################################
#parameter preparation
########################################################################
p00d=rep(0,ncol(xij))
p0d=rep(0,studynum)
p0dz=rep(0,studynum)
z<-matrix(0,sitenum,studynum)
b<-matrix(0,sitenum,studynum)
B<-matrix(0,sitenum,studynum)
c0d<-rep(1,studynum)
d0d<-rep(1,studynum)
xid<-matrix(0,nrow=sitenum,ncol=studynum)
nid<-matrix(0,nrow=sitenum,ncol=studynum)
#pool counts
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]+tempx
tempn=nij[,sampleid]+tempn
}
xid[,d]=tempx
nid[,d]=tempn
prob=(tempx)/(tempn)
p00d[d]=mean(prob[which(tempn!=0)])
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]+tempx
tempn=nij[,sampleid]+tempn
}
#naive z statistic
p0dz[d]=mean(tempx/tempn,na.rm=TRUE)
z[,d]=abs(tempx/tempn-p0dz[d])
#naive Bayes statistic
addterm=p00d[d]
prob=(tempx+2*addterm)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0d[d]=mu*(mu*(1-mu)/v-1)
d0d[d]=(1-mu)*(mu*(1-mu)/v-1)
p0d[d]=mean(prob[which(tempn!=0)])
b[,d]=abs((tempx+2*addterm)/(tempn+2)-p0d[d])
B[,d]=abs((tempx+c0d[d])/(tempx+c0d[d]+d0d[d])-p0d[d])
}
result<-list(z=z,b=b,B=B,c0d=c0d,d0d=d0d,p0d=p0d,p0dz=p0dz)
}
iASeqmotif<-function(exprs,studyid,repid,refid,K,iter.max=100,tol=1e-3)
{
fitresult<-list()
for(i in 1:length(K))
{
fitresult[[i]]<-iASeqmotiffit(exprs,studyid,repid,refid,K[i],iter.max=iter.max,tol=tol)
}
if(K[1]>1)
{
studynum=ncol(fitresult[[1]]$motif.qup)
}
else
{
studynum=length(fitresult[[1]]$motif.qup)
}
sitenum=nrow(exprs)
bic<-rep(0,length(K))
loglike<-rep(0,length(K))
for(i in 1:length(K))
{
loglike[i]<-fitresult[[i]]$loglike
}
for(i in 1:length(K))
{
bic[i]<--2*fitresult[[i]]$loglike+(K[i]+K[i]*studynum*2+2*ncol(exprs)/2)*log(sitenum)
}
bestflag=which(bic==min(bic))
result<-list(bestmotif=fitresult[[bestflag]],bic=cbind(K,bic),loglike=cbind(K,loglike))
}
###########################################################
#Motif plot
###########################################################
plotMotif<-function(bestmotif,title="",cutoff)
{
par(mai=c(0.7,1,0.7,0.2))
layout(matrix(1:4,ncol=4))
u<-1:dim(bestmotif$motif.qup)[2]
v<-1:dim(bestmotif$motif.qup)[1]
image(u,v,t(bestmotif$motif.qup),
col=gray(seq(from=1,to=0,by=-0.1)),xlab="StudyID",yaxt = "n",
ylab="Motifs",main=paste(title,"V",sep=""),
cex.axis=2,cex.lab=2,cex.main=2)
axis(2,at=1:length(v),,cex.axis=2,cex.lab=3)
for(i in 1:(length(u)+1))
{
abline(v=(i-0.5))
}
for(i in 1:(length(v)+1))
{
abline(h=(i-0.5))
}
image(u,v,t(bestmotif$motif.qdown),
col=gray(seq(from=1,to=0,by=-0.1)),xlab="StudyID",yaxt = "n",
ylab="Motifs",main=paste(title,"W",sep=""),
cex.axis=2,cex.lab=2,cex.main=2)
axis(2,at=1:length(v),cex.axis=2,cex.lab=3)
for(i in 1:(length(u)+1))
{
abline(v=(i-0.5))
}
for(i in 1:(length(v)+1))
{
abline(h=(i-0.5))
}
Ng=nrow(bestmotif$clustlike)
genecount=floor(bestmotif$motif.p[-1]*Ng)
NK=nrow(bestmotif$motif.qup)
plot(0,0.7,pch=".",xlim=c(0,1.2),ylim=c(0.75,NK+0.25),frame.plot=FALSE,axes=FALSE,xlab="No. of SNPs",ylab=""
,cex.axis=1.5,cex.lab=2,cex.main=3,main=expression(pi))
segments(0,0.7,bestmotif$motif.p[2],0.7)
rect(0,1:NK-0.3,bestmotif$motif.p[-1],1:NK+0.3,col="dark grey")
mtext(1:NK,at=1:NK,side=2,cex=1.5)
text(bestmotif$motif.p[-1]+0.55,1:NK,labels=floor(bestmotif$motif.p[-1]*Ng),cex=1)
clustnum=sapply(1:NK, function(k)
length(which(sapply(1:nrow(bestmotif$clustlike), function(i) which.max(bestmotif$clustlike[i,]))==(k+1) &
bestmotif$clustlike[,(k+1)]>cutoff)))
plot(0,0.7,pch=".",xlim=c(0,1.2),ylim=c(0.75,NK+0.25),frame.plot=FALSE,axes=FALSE,xlab="No. of SNPs",ylab=""
,cex.axis=1.5,cex.lab=2,cex.main=3,main=expression(a[i]))
segments(0,0.7,clustnum[1]/Ng,0.7)
rect(0,1:NK-0.3,clustnum/Ng,1:NK+0.3,col="dark grey")
mtext(1:NK,at=1:NK,side=2,cex=1.5)
text(clustnum/Ng+0.55,1:NK,labels=clustnum,cex=1)
}
###########################################################
#BIC plot
###########################################################
plotBIC<-function(fitted_cormotif)
{
plot(fitted_cormotif$bic[,1], fitted_cormotif$bic[,2], type="b", xlab="Non-null Motif Number", ylab="BIC", main="BIC")
lines(fitted_cormotif$bic[,1], fitted_cormotif$bic[,2])
}
singleEMfit<-function(exprs,studyid,repid,refid,iter.max=100,tol=1e-3)
{
set.seed(1)
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
#result<-list(xij=xij,nij=nij,design=experiment.unique)
########################################################################
#parameter preparation
########################################################################
p0j=rep(0,ncol(xij))
p00j=rep(0,ncol(xij))
#need fitting later
c0j<-rep(1,ncol(xij))
d0j<-rep(1,ncol(xij))
cd<-rep(1,studynum)
dd<-rep(1,studynum)
#initialization
qup=rep(0.1,studynum)
qdown=rep(0.1,studynum)
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
prob=(tempx)/(tempn)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
p00j[sampleid]=mean(prob[which(tempn!=0)])
}
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
addterm=p00j[sampleid]
prob=(tempx+addterm*2)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0j[sampleid]=mu*(mu*(1-mu)/v-1)
d0j[sampleid]=(1-mu)*(mu*(1-mu)/v-1)
p0j[sampleid]=mean(prob[which(tempn!=0)])
}
}
c0j=abs(c0j)
d0j=abs(d0j)
loglike.old<--1e10
########################################################################
#model fitting
########################################################################
loglike0<-list()
loglikeup<-list()
loglikedown<-list()
for(i.iter in 1:iter.max)
{
err<-tol+1
for(d in 1:studynum)
{
temploglike0=rep(0,sitenum)
temploglikeup=rep(0,sitenum)
temploglikedown=rep(0,sitenum)
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
sampleid=repseq[j]
temploglike0<-temploglike0+f0.loglike(xij[,sampleid],nij[,sampleid],c0j[sampleid],d0j[sampleid])
temploglikeup<-temploglikeup+fup.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
temploglikedown<-temploglikedown+fdown.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
}
loglike0[[d]]<-temploglike0
loglikeup[[d]]<-temploglikeup
loglikedown[[d]]<-temploglikedown
}
#compute posterior cluster membership
condlike<-list()
for(d in 1:studynum)
{
condlike[[d]]<-matrix(0,sitenum,3)
}
templike<-matrix(0,sitenum,3)
for(d in 1:studynum)
{
templike[,1]<-log(qup[d])+loglikeup[[d]]
templike[,2]<-log(1-qup[d]-qdown[d])+loglike0[[d]]
templike[,3]<-log(qdown[d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
condlike[[d]][,1]<-templike[,1]/tempsum
condlike[[d]][,2]<-templike[,2]/tempsum
condlike[[d]][,3]<-templike[,3]/tempsum
}
qup.new<-rep(0,studynum)
qdown.new<-rep(0,studynum)
for(d in 1:studynum)
{
qup.new[d]<-(sum(condlike[[d]][,1])+1)/(sitenum+3)
qdown.new[d]<-(sum(condlike[[d]][,3])+1)/(sitenum+3)
}
loglike.new<-(sum(tempmax+log(tempsum))
+sum(log(qup.new)+log(qdown.new)+log(1-qup.new-qdown.new)))/sitenum
## evaluate convergence
err.qup<-max(abs(qup.new-qup)/qup)
err.qdown<-max(abs(qdown.new-qdown)/qdown)
err<-max(err.qup,err.qdown)
qup<-qup.new
qdown<-qdown.new
loglike.old<-loglike.new
if(err<tol)
{
break
}
}
for(d in 1:studynum)
{
condlike[[d]]<-matrix(0,sitenum,3)
}
templike<-matrix(0,sitenum,3)
for(d in 1:studynum)
{
templike[,1]<-log(qup[d])+loglikeup[[d]]
templike[,2]<-log(1-qup[d]-qdown[d])+loglike0[[d]]
templike[,3]<-log(qdown[d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
condlike[[d]][,1]<-templike[,1]/tempsum
condlike[[d]][,2]<-templike[,2]/tempsum
condlike[[d]][,3]<-templike[,3]/tempsum
}
p.post2=matrix(0,sitenum,studynum)
for(d in 1:studynum)
{
p.post2[,d]=condlike[[d]][,2]
}
result<-list(p.study=p.post2,motif.qup=qup,motif.qdown=qdown,condlike=condlike)
}
Try the iASeq package in your browser
Any scripts or data that you put into this service are public.
iASeq documentation built on Nov. 8, 2020, 5:36 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions singleEMfit plotBIC plotMotif iASeqmotif ASErawfit iASeqmotiffit fdown.loglike fup.loglike f0.loglike
Documented in ASErawfit f0.loglike fdown.loglike fup.loglike iASeqmotif iASeqmotiffit plotBIC plotMotif singleEMfit
f0.loglike<-function(xj,nj,cj,dj)
{
lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
fup.loglike<-function(xj,nj,cj,dj,pj)
{
log(1/(1-pj))+log(1-pbeta(pj,xj+cj,nj-xj+dj))+lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
fdown.loglike<-function(xj,nj,cj,dj,pj)
{
log(1/pj)+log(pbeta(pj,xj+cj,nj-xj+dj))+lbeta(xj+cj,nj-xj+dj)-lbeta(cj,dj)
}
iASeqmotiffit<-function(exprs,studyid,repid,refid,K,iter.max=100,tol=1e-3)
{
set.seed(1)
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
########################################################################
#parameter preparation
########################################################################
pai=rep(0,K+1)
p0j=rep(0,ncol(xij))
p00j=rep(0,ncol(xij))
#need fitting later
c0j<-rep(1,ncol(xij))
d0j<-rep(1,ncol(xij))
cd<-rep(1,studynum)
dd<-rep(1,studynum)
#initialization
pai<-c(0.8,rep(0.2/K,K))
q0<-runif(K*studynum)
q1<-runif(K*studynum)
q2<-runif(K*studynum)
qsum=(q0+q1+q2)
qup=matrix(q1/qsum,K,studynum)
qdown=matrix(q2/qsum,K,studynum)
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
prob=(tempx)/(tempn)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
p00j[sampleid]=mean(prob[which(tempn!=0)])
}
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
addterm=p00j[sampleid]
prob=(tempx+addterm*2)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0j[sampleid]=mu*(mu*(1-mu)/v-1)
d0j[sampleid]=(1-mu)*(mu*(1-mu)/v-1)
p0j[sampleid]=mean(prob[which(tempn!=0)])
}
}
c0j=abs(c0j)
d0j=abs(d0j)
loglike.old<--1e10
########################################################################
#model fitting
########################################################################
loglike0<-list()
loglikeup<-list()
loglikedown<-list()
for(i.iter in 1:iter.max)
{
if(i.iter%%50==0)
{
print(paste("We have run ",i.iter," iterations for K=",K, sep=""))
}
err<-tol+1
for(d in 1:studynum)
{
temploglike0=rep(0,sitenum)
temploglikeup=rep(0,sitenum)
temploglikedown=rep(0,sitenum)
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
sampleid=repseq[j]
temploglike0<-temploglike0+f0.loglike(xij[,sampleid],nij[,sampleid],c0j[sampleid],d0j[sampleid])
temploglikeup<-temploglikeup+fup.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
temploglikedown<-temploglikedown+fdown.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
}
loglike0[[d]]<-temploglike0
loglikeup[[d]]<-temploglikeup
loglikedown[[d]]<-temploglikedown
}
clustlike<-matrix(0,sitenum,K+1)
condlikeup<-list()
condlikedown<-list()
for(d in 1:studynum)
{
condlikeup[[d]]<-matrix(0,sitenum,K+1)
condlikedown[[d]]<-matrix(0,sitenum,K+1)
}
#for ai=0 case
for(d in 1:studynum)
{
clustlike[,1]=clustlike[,1]+loglike0[[d]]
}
clustlike[,1]=clustlike[,1]+log(pai[1])
#for ai!=0 case
templike<-matrix(0,sitenum,3)
for(k in 2:(K+1))
{
for(d in 1:studynum)
{
templike[,1]<-log(qup[k-1,d])+loglikeup[[d]]
templike[,2]<-log(1-qup[k-1,d]-qdown[k-1,d])+loglike0[[d]]
templike[,3]<-log(qdown[k-1,d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
clustlike[,k]<-clustlike[,k]+tempmax+log(tempsum)
condlikeup[[d]][,k]<-templike[,1]/tempsum
condlikedown[[d]][,k]<-templike[,3]/tempsum
}
clustlike[,k]<-clustlike[,k]+log(pai[k])
}
tempmax<-apply(clustlike,1,max)
for(k in 1:(K+1))
{
clustlike[,k]<-exp(clustlike[,k]-tempmax)
}
tempsum<-apply(clustlike,1,sum)
#update motif occurrence rate
for(k in 1:(K+1))
{
clustlike[,k]<-clustlike[,k]/tempsum
}
pai.new<-(apply(clustlike,2,sum)+1)/(sitenum+K+1)
#update motifs
qup.new<-matrix(0,K,studynum)
qdown.new<-matrix(0,K,studynum)
for(k in 1:K)
{
clustpsum<-sum(clustlike[,k+1])
for(d in 1:studynum)
{
qup.new[k,d]<-(sum(clustlike[,k+1]*condlikeup[[d]][,k+1])+1)/(clustpsum+3)
qdown.new[k,d]<-(sum(clustlike[,k+1]*condlikedown[[d]][,k+1])+1)/(clustpsum+3)
}
}
loglike.new<-(sum(tempmax+log(tempsum))+sum(log(pai.new))
+sum(log(qup.new)+log(qdown.new)+log(1-qup.new-qdown.new)))/sitenum
## evaluate convergence
err.pai<-max(abs(pai.new-pai)/pai)
err.qup<-max(abs(qup.new-qup)/qup)
err.qdown<-max(abs(qdown.new-qdown)/qdown)
err<-max(err.pai, err.qup,err.qdown)
pai<-pai.new
qup<-qup.new
qdown<-qdown.new
loglike.old<-loglike.new
if(err<tol)
{
break
}
}
for(d in 1:studynum)
{
condlikeup[[d]]<-matrix(0,sitenum,K+1)
condlikedown[[d]]<-matrix(0,sitenum,K+1)
}
#for ai=0 case
for(d in 1:studynum)
{
clustlike[,1]=clustlike[,1]+loglike0[[d]]
}
clustlike[,1]=clustlike[,1]+log(pai[1])
#for ai!=0 case
templike<-matrix(0,sitenum,3)
for(k in 2:(K+1))
{
for(d in 1:studynum)
{
templike[,1]<-log(qup[k-1,d])+loglikeup[[d]]
templike[,2]<-log(1-qup[k-1,d]-qdown[k-1,d])+loglike0[[d]]
templike[,3]<-log(qdown[k-1,d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
clustlike[,k]<-clustlike[,k]+tempmax+log(tempsum)
condlikeup[[d]][,k]<-templike[,1]/tempsum
condlikedown[[d]][,k]<-templike[,3]/tempsum
}
clustlike[,k]<-clustlike[,k]+log(pai[k])
}
tempmax<-apply(clustlike,1,max)
for(k in 1:(K+1))
{
clustlike[,k]<-exp(clustlike[,k]-tempmax)
}
tempsum<-apply(clustlike,1,sum)
for(k in 1:(K+1))
{
clustlike[,k]<-clustlike[,k]/tempsum
}
p.post=1-(clustlike[,1])
minus=sum(log(pai))+sum(log(qup)+log(qdown)+log(1-qup-qdown))
result<-list(p.post=p.post, motif.prior=pai,motif.qup=qup,motif.qdown=qdown,clustlike=clustlike,
c0j=c0j,d0j=d0j,loglike=loglike.new*sitenum-minus)
}
ASErawfit<-function(exprs,studyid,repid,refid)
{
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
#result<-list(xij=xij,nij=nij,design=experiment.unique)
########################################################################
#parameter preparation
########################################################################
p00d=rep(0,ncol(xij))
p0d=rep(0,studynum)
p0dz=rep(0,studynum)
z<-matrix(0,sitenum,studynum)
b<-matrix(0,sitenum,studynum)
B<-matrix(0,sitenum,studynum)
c0d<-rep(1,studynum)
d0d<-rep(1,studynum)
xid<-matrix(0,nrow=sitenum,ncol=studynum)
nid<-matrix(0,nrow=sitenum,ncol=studynum)
#pool counts
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]+tempx
tempn=nij[,sampleid]+tempn
}
xid[,d]=tempx
nid[,d]=tempn
prob=(tempx)/(tempn)
p00d[d]=mean(prob[which(tempn!=0)])
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]+tempx
tempn=nij[,sampleid]+tempn
}
#naive z statistic
p0dz[d]=mean(tempx/tempn,na.rm=TRUE)
z[,d]=abs(tempx/tempn-p0dz[d])
#naive Bayes statistic
addterm=p00d[d]
prob=(tempx+2*addterm)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0d[d]=mu*(mu*(1-mu)/v-1)
d0d[d]=(1-mu)*(mu*(1-mu)/v-1)
p0d[d]=mean(prob[which(tempn!=0)])
b[,d]=abs((tempx+2*addterm)/(tempn+2)-p0d[d])
B[,d]=abs((tempx+c0d[d])/(tempx+c0d[d]+d0d[d])-p0d[d])
}
result<-list(z=z,b=b,B=B,c0d=c0d,d0d=d0d,p0d=p0d,p0dz=p0dz)
}
iASeqmotif<-function(exprs,studyid,repid,refid,K,iter.max=100,tol=1e-3)
{
fitresult<-list()
for(i in 1:length(K))
{
fitresult[[i]]<-iASeqmotiffit(exprs,studyid,repid,refid,K[i],iter.max=iter.max,tol=tol)
}
if(K[1]>1)
{
studynum=ncol(fitresult[[1]]$motif.qup)
}
else
{
studynum=length(fitresult[[1]]$motif.qup)
}
sitenum=nrow(exprs)
bic<-rep(0,length(K))
loglike<-rep(0,length(K))
for(i in 1:length(K))
{
loglike[i]<-fitresult[[i]]$loglike
}
for(i in 1:length(K))
{
bic[i]<--2*fitresult[[i]]$loglike+(K[i]+K[i]*studynum*2+2*ncol(exprs)/2)*log(sitenum)
}
bestflag=which(bic==min(bic))
result<-list(bestmotif=fitresult[[bestflag]],bic=cbind(K,bic),loglike=cbind(K,loglike))
}
###########################################################
#Motif plot
###########################################################
plotMotif<-function(bestmotif,title="",cutoff)
{
par(mai=c(0.7,1,0.7,0.2))
layout(matrix(1:4,ncol=4))
u<-1:dim(bestmotif$motif.qup)[2]
v<-1:dim(bestmotif$motif.qup)[1]
image(u,v,t(bestmotif$motif.qup),
col=gray(seq(from=1,to=0,by=-0.1)),xlab="StudyID",yaxt = "n",
ylab="Motifs",main=paste(title,"V",sep=""),
cex.axis=2,cex.lab=2,cex.main=2)
axis(2,at=1:length(v),,cex.axis=2,cex.lab=3)
for(i in 1:(length(u)+1))
{
abline(v=(i-0.5))
}
for(i in 1:(length(v)+1))
{
abline(h=(i-0.5))
}
image(u,v,t(bestmotif$motif.qdown),
col=gray(seq(from=1,to=0,by=-0.1)),xlab="StudyID",yaxt = "n",
ylab="Motifs",main=paste(title,"W",sep=""),
cex.axis=2,cex.lab=2,cex.main=2)
axis(2,at=1:length(v),cex.axis=2,cex.lab=3)
for(i in 1:(length(u)+1))
{
abline(v=(i-0.5))
}
for(i in 1:(length(v)+1))
{
abline(h=(i-0.5))
}
Ng=nrow(bestmotif$clustlike)
genecount=floor(bestmotif$motif.p[-1]*Ng)
NK=nrow(bestmotif$motif.qup)
plot(0,0.7,pch=".",xlim=c(0,1.2),ylim=c(0.75,NK+0.25),frame.plot=FALSE,axes=FALSE,xlab="No. of SNPs",ylab=""
,cex.axis=1.5,cex.lab=2,cex.main=3,main=expression(pi))
segments(0,0.7,bestmotif$motif.p[2],0.7)
rect(0,1:NK-0.3,bestmotif$motif.p[-1],1:NK+0.3,col="dark grey")
mtext(1:NK,at=1:NK,side=2,cex=1.5)
text(bestmotif$motif.p[-1]+0.55,1:NK,labels=floor(bestmotif$motif.p[-1]*Ng),cex=1)
clustnum=sapply(1:NK, function(k)
length(which(sapply(1:nrow(bestmotif$clustlike), function(i) which.max(bestmotif$clustlike[i,]))==(k+1) &
bestmotif$clustlike[,(k+1)]>cutoff)))
plot(0,0.7,pch=".",xlim=c(0,1.2),ylim=c(0.75,NK+0.25),frame.plot=FALSE,axes=FALSE,xlab="No. of SNPs",ylab=""
,cex.axis=1.5,cex.lab=2,cex.main=3,main=expression(a[i]))
segments(0,0.7,clustnum[1]/Ng,0.7)
rect(0,1:NK-0.3,clustnum/Ng,1:NK+0.3,col="dark grey")
mtext(1:NK,at=1:NK,side=2,cex=1.5)
text(clustnum/Ng+0.55,1:NK,labels=clustnum,cex=1)
}
###########################################################
#BIC plot
###########################################################
plotBIC<-function(fitted_cormotif)
{
plot(fitted_cormotif$bic[,1], fitted_cormotif$bic[,2], type="b", xlab="Non-null Motif Number", ylab="BIC", main="BIC")
lines(fitted_cormotif$bic[,1], fitted_cormotif$bic[,2])
}
singleEMfit<-function(exprs,studyid,repid,refid,iter.max=100,tol=1e-3)
{
set.seed(1)
########################################################################
#data preparation
########################################################################
sitenum<-nrow(exprs)
studynum<-length(unique(studyid))
experiment.unique<-unique(cbind(studyid,repid))
colnames(experiment.unique)<-c("study","rep")
xij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
nij<-matrix(0,nrow=sitenum,ncol=nrow(experiment.unique))
#extract reference allele frequencies
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab & refid==0)
xij[,flag]=exprs[,sampleid]
flag=flag+1
}
}
#extract total count for each allele in each sample
flag=1
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
replab=experiment.unique[repseq,2][j]
sampleid=which(studyid==studylab & repid==replab)
nij[,flag]=rowSums(exprs[,sampleid])
flag=flag+1
}
}
#result<-list(xij=xij,nij=nij,design=experiment.unique)
########################################################################
#parameter preparation
########################################################################
p0j=rep(0,ncol(xij))
p00j=rep(0,ncol(xij))
#need fitting later
c0j<-rep(1,ncol(xij))
d0j<-rep(1,ncol(xij))
cd<-rep(1,studynum)
dd<-rep(1,studynum)
#initialization
qup=rep(0.1,studynum)
qdown=rep(0.1,studynum)
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
prob=(tempx)/(tempn)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
p00j[sampleid]=mean(prob[which(tempn!=0)])
}
}
for(d in 1:studynum)
{
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
tempx=rep(0,sitenum)
tempn=rep(0,sitenum)
for(j in 1:dn)
{
sampleid=repseq[j]
tempx=xij[,sampleid]
tempn=nij[,sampleid]
addterm=p00j[sampleid]
prob=(tempx+addterm*2)/(tempn+2)
mu=mean(prob[which(tempn!=0)])
v=var(prob[which(tempn!=0)])
c0j[sampleid]=mu*(mu*(1-mu)/v-1)
d0j[sampleid]=(1-mu)*(mu*(1-mu)/v-1)
p0j[sampleid]=mean(prob[which(tempn!=0)])
}
}
c0j=abs(c0j)
d0j=abs(d0j)
loglike.old<--1e10
########################################################################
#model fitting
########################################################################
loglike0<-list()
loglikeup<-list()
loglikedown<-list()
for(i.iter in 1:iter.max)
{
err<-tol+1
for(d in 1:studynum)
{
temploglike0=rep(0,sitenum)
temploglikeup=rep(0,sitenum)
temploglikedown=rep(0,sitenum)
studylab=unique(studyid)[d]
repseq=which(experiment.unique[,1]==studylab)
dn<-length(repseq)
for(j in 1:dn)
{
sampleid=repseq[j]
temploglike0<-temploglike0+f0.loglike(xij[,sampleid],nij[,sampleid],c0j[sampleid],d0j[sampleid])
temploglikeup<-temploglikeup+fup.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
temploglikedown<-temploglikedown+fdown.loglike(xij[,sampleid],nij[,sampleid],cd[d],dd[d],p0j[sampleid])
}
loglike0[[d]]<-temploglike0
loglikeup[[d]]<-temploglikeup
loglikedown[[d]]<-temploglikedown
}
#compute posterior cluster membership
condlike<-list()
for(d in 1:studynum)
{
condlike[[d]]<-matrix(0,sitenum,3)
}
templike<-matrix(0,sitenum,3)
for(d in 1:studynum)
{
templike[,1]<-log(qup[d])+loglikeup[[d]]
templike[,2]<-log(1-qup[d]-qdown[d])+loglike0[[d]]
templike[,3]<-log(qdown[d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
condlike[[d]][,1]<-templike[,1]/tempsum
condlike[[d]][,2]<-templike[,2]/tempsum
condlike[[d]][,3]<-templike[,3]/tempsum
}
qup.new<-rep(0,studynum)
qdown.new<-rep(0,studynum)
for(d in 1:studynum)
{
qup.new[d]<-(sum(condlike[[d]][,1])+1)/(sitenum+3)
qdown.new[d]<-(sum(condlike[[d]][,3])+1)/(sitenum+3)
}
loglike.new<-(sum(tempmax+log(tempsum))
+sum(log(qup.new)+log(qdown.new)+log(1-qup.new-qdown.new)))/sitenum
## evaluate convergence
err.qup<-max(abs(qup.new-qup)/qup)
err.qdown<-max(abs(qdown.new-qdown)/qdown)
err<-max(err.qup,err.qdown)
qup<-qup.new
qdown<-qdown.new
loglike.old<-loglike.new
if(err<tol)
{
break
}
}
for(d in 1:studynum)
{
condlike[[d]]<-matrix(0,sitenum,3)
}
templike<-matrix(0,sitenum,3)
for(d in 1:studynum)
{
templike[,1]<-log(qup[d])+loglikeup[[d]]
templike[,2]<-log(1-qup[d]-qdown[d])+loglike0[[d]]
templike[,3]<-log(qdown[d])+loglikedown[[d]]
tempmax<-pmax(templike[,1],templike[,2],templike[,3])
for(z in 1:3)
{
templike[,z]<-exp(templike[,z]-tempmax)
}
tempsum<-templike[,1]+templike[,2]+templike[,3]
condlike[[d]][,1]<-templike[,1]/tempsum
condlike[[d]][,2]<-templike[,2]/tempsum
condlike[[d]][,3]<-templike[,3]/tempsum
}
p.post2=matrix(0,sitenum,studynum)
for(d in 1:studynum)
{
p.post2[,d]=condlike[[d]][,2]
}
result<-list(p.study=p.post2,motif.qup=qup,motif.qdown=qdown,condlike=condlike)
}
Try the iASeq package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.