# R/kaksCodon.R In CorMut: Detect the correlated mutations based on selection pressure

kaksCodon <-
function(seq_formated){
seq_char=tolower(as.character(seq_formated))
refseq=seq_char[1]
seq_list=seq_char[2:length(seq_char)]
mat=sapply(as.list(seq_list),function(x)translate_cm(s2c(x)))
refseq=tolower(refseq)
refaa=translate_cm(s2c(refseq))
fvft=fvft(seq_list,refseq=refseq)
#final=list()
kaks_result=c()
LOD_result=c()
q_result=c()
KA_result=c()
KS_result=c()
for(i in 1:dim(mat)[1]){
NS=0
specific_row=mat[i,]
NY=sum(specific_row!=refaa[i])
KA_result=append(KA_result,NY)
for(k in 1:dim(mat)[2]){
codon_nu=substr(seq_list[[k]],(i-1)*3+1,i*3)
#cat(codon_nu,substr(refseq,(i-1)*3+1,i*3),"\n")
#cat("\n")
#cat(codon_nu!=substr(refseq,(i-1)*3+1,i*3))
if(codon_nu!=substr(refseq,(i-1)*3+1,i*3)&specific_row[k]==refaa[i]){
#cat(specific_row[k],refaa[i],specific_row[k]==refaa[i],"\n")
NS=NS+1
}
}
KS_result=append(KS_result,NS)
#compute NtNv
targetbp=s2c(refseq)[((i-1)*3+1):(i*3)]
#print(targetbp)
NtNv=Ntnv(targetbp)
model_fenzi=NtNv\$Nat*fvft\$ft+NtNv\$Nav*fvft\$fv
model_fenmu=NtNv\$Nst*fvft\$ft+NtNv\$Nsv*fvft\$fv
if(model_fenmu==0|model_fenzi==0){
random_model=((NtNv\$Nat+1)*fvft\$ft+(NtNv\$Nav+1)*fvft\$fv)/((NtNv\$Nst+1)*fvft\$ft+(NtNv\$Nsv+1)*fvft\$fv)
}else{
random_model=model_fenzi/model_fenmu
}
#random_model=(NtNv\$Nat*fvft\$ft+NtNv\$Nav*fvft\$fv)/(NtNv\$Nst*fvft\$ft+NtNv\$Nsv*fvft\$fv)
#if((random_model==Inf)|(random_model=="NaN")|(random_model==0))random_model=1
if(NY==0&NS==0){
KaKs=1
}else if(NY!=0&NS==0){
KaKs=((NY+1)/(NS+1))/random_model
}else{
KaKs=(NY/NS)/random_model
}
kaks_result=append(kaks_result,KaKs)
q=(NtNv\$Nat*fvft\$ft+NtNv\$Nav*fvft\$fv)/(3*fvft\$ft+6*fvft\$fv)
q_result=append(q_result,q)
LOD=0
#print(NY)
Nnum=NS+NY
for(ii in NY:Nnum){
#LOD=choose(Nnum,ii)*(q^ii)*((1-q)^(Nnum-ii))+LOD
LOD=dbinom(ii,Nnum,q)+LOD
#print(choose(dim(mat)[2],i)*(q^i)*((1-q)^(dim(mat)[2]-i)))
}
#print(LOD)
LOD_result=append(LOD_result,-log10(LOD))
}
final=new("kaksCodon",seq_num=length(seq_list),kaks=kaks_result,lod=LOD_result,q=q_result,ka=KA_result,ks=KS_result)
return(final)
}

## Try the CorMut package in your browser

Any scripts or data that you put into this service are public.

CorMut documentation built on May 6, 2019, 4:56 a.m.