###########################################################################
# R codes for the Poisson-PDNN model
# Manuscript: Modeling non-uniform read distribution based on the PDNN model in RNA-Seq data
# Authors:    Yuchao Xia”č, Fugui Wang”č, Minping Qian, Zhaohui Qin, MinghuaDeng
# Date:       2012.12.13
###########################################################################

###########################################################################
# This function makes use of Newton algorithm for cross-validation. 
# And the input of this function are both the data after using the 
# function of expData which is provided in mseq Package and the initial 
# parameters. The output of this function is R squared and the coefficient.
###########################################################################
`NewtonCV` <-
function(data,b,mu,fold = 5, seed = 281142, eps=0.01,iter=1000)
{
	cat("Cross validation using iterative glm...\n")
	train_index <- CVTrainIndexGene(data, fold, seed)
	n<-length(b)-1
	dev <- rep(0, fold)
	null_dev <- rep(0, fold)
	coef<-rep(0,fold)
	for (j in 1 : fold)
	{
		
            cat("\nFold", j, "...\n")
            log_pref <- getOriLogPref(data[train_index[j, ], ], 0.5)
            X<-as.matrix(data[train_index[j,],3:dim(data)[2]]/1000000)
     
            log_expr<-setOriOffset(data[train_index[j,],])
            log_pref1<-log_pref-log_expr
           
            #for(i in 1:max_iter)
            #{
            #cat("max_iter:",i,"\n")
            out<-Newton(b,mu,X,log_pref1,eps,iter)
            coef<-out$coef
            #predCount<-exp(coef[1]+1/(1+exp(X%*%coef[-1])))
            #log_expr<-updateOffset(data[train_index[j,],],predCount,log_expr)
            #log_pref1<-log_pref1-log_expr
            pred_count<-getPredCount_pdnn(data[train_index[j,],],X,log_expr,coef)
            nowdev<-getDev(pred_count,data[train_index[j,],]$count)
            cat("deviance=",nowdev,"\n")
            #b<-coef
            #if(i>1&&abs((nowdev-fordev)/fordev<thrd))
            #{
            #cat("dev done...\n")
            #break
            #}
            #fordev<-nowdev
            #}
            b<-out$coef
            X1<-as.matrix(data[!train_index[j,],3:dim(data)[2]]/1000000)
            log_expr_test<-setOriOffset(data[!train_index[j,],])
            pred_count<-getPredCount_pdnn(data[!train_index[j, ], ], X1,log_expr_test,out$coef)	 	
	      dev[j] <- getDev(pred_count, data[!train_index[j, ], ]$count)
       	null_dev[j] <- getDev(getNullCount_pdnn(data[!train_index[j, ], ],X1,n), data[!train_index[j, ], ]$count)
		cat("\nFold ", j, ": null deviance = ", null_dev[j], "; deviance = ", dev[j], "; R squared = ", 1 - dev[j] / null_dev[j], "\n", sep = '')
	}
	
	cat("\n")
	cat("deviance =", dev, "\n")
	cat("null deviance =", null_dev, "\n")
	coef<-out$coef
	R_squared <- 1 - sum(dev) / sum(null_dev)
	cat("R_squared =", R_squared, "\n")
	ls<-list(R_squared=R_squared,coef=coef)
	return(ls)
}


#########################################################################################
# This function is used to optimize the least squares method with penalty function to find optimal solution. 
# The inputs of this function are both weight parameters and the matrix of the data after
# using the function of expData which is provided in mseq Package.
#########################################################################################
Newton<-function(b,mu,X,Y,eps=0.01,iter=1000)
{
n<-length(b)
E<-matrix(0,n,n)
diag(E)<-1
Hessien<-matrix(0,n,n)
for(i in 1:iter)
{
cat("iter:...",i,"\n")
#hession
g1<-sum(-2*(Y-b[1]-1/(1+exp(X%*%cbind(b[-1])))))
g2<-t(cbind(2*(Y-b[1]-1/(1+exp(X%*%cbind(b[-1]))))*exp(X%*%cbind(b[-1]))/(1+exp(X%*%cbind(b[-1])))^2))%*%X-mu*(b[-1])^(-1)
g<-t(cbind(g1,g2))
g11<-2*length(Y)
g12<-t((-2)*exp(X%*%cbind(b[-1]))/(1+exp(X%*%cbind(b[-1])))^2)%*%X
tmp<-as.vector((2*(Y-b[1]-1/(1+exp(X%*%b[-1])))*((exp(X%*%b[-1]))^2-exp(X%*%b[-1]))/(1+exp(X%*%b[-1]))^3))
tmp_matrix<-matrix(rep(tmp,ncol(X)),nrow=ncol(X),byrow=T)
tmp1<-(t(X)*tmp_matrix)%*%X
g22<-2*t(t(exp(X%*%b[-1])/(1+exp(X%*%b[-1]))^2)%*%X)%*%(t(exp(X%*%b[-1])/(1+exp(X%*%b[-1]))^2)%*%X)+tmp1
g22<-g22+mu*diag(b[-1]^(-2))
Hessien[1,1]<-g11
Hessien[1,2:n]<-g12
Hessien[2:n,1]<-g12
Hessien[2:n,2:n]<-g22
H<-chol(Hessien+E,pivot=TRUE)
HH<-H%*%t(H)
b1<-b
#d<-(-solve(HH)%*%g)
b<-b1-solve(HH)%*%g
mu=mu*0.8
#result<-search(f,1,4,b,X,Y,d)
#lemda<-ercichazhi(f,result,b,X,Y,d)
#b<-b1+lemda*d
res<-(Y-b[1]-1/(1+exp(X%*%b[-1])))
dev<-sum(res^2)
cat("dev=",dev,"\n")
if(sum((b-b1)^2)<=eps)
{
cat("Done...\n")
break
}
if(sum(g^2)<=0.05)
{
cat("gDone...\n")
break
}
}
result<-list(coef=b,iter=iter,res=dev)
return(result)
}


##########################################################################################
# Get the counts of reads for each position under the null hypothesis (each positions have 
# the same weights). The output is null count.
###########################################################################################
`getNullCount_pdnn` <-
function(data,X,n)
{
if (n == 16)
{
b=rep(1/n,n)
}
else
{
b=rep(1/n,n)
}
Null_count<-exp(setOriOffset(data))*exp(1+1/(1+exp(X%*%b)))
return(Null_count)
}

#######################################################################################
# Get the predicted counts using the results of coefficient after using the function of Newton.
# The output is the predicted count
#######################################################################################
getPredCount_pdnn <-function(data,X,log_expr,coef)
{
    #pred_pref<-exp(log_expr)*exp(1/(1+exp(X%*%coef[2:length(coef)])))
	pred_pref<-exp(1/(1+exp(X%*%coef[2:length(coef)])))
	uniq_cat <- unique(data$index)
	pred_count <- rep(0, length(data$count))
	for (k in uniq_cat)
	{
		pred_count[data$index == k]<- sum((data$count[data$index == k])) / sum(pred_pref[data$index == k]) * pred_pref[data$index == k]
	}
	return(pred_count)
}



##################################################################################################
# Calculate the deviance under Poisson-PDNN model. The output is deviance.
##################################################################################################
getDev_pdnn <-function(pred_count, real_count)
{
dev<-2*sum(real_count[real_count!=0]*log(real_count[real_count!=0]/pred_count)+pred_count-real_count)	
return(dev)
}



####################################################################################################
# This function is used of calculating the energy between two nucleotides and gets a matrix of energy.
# The initial values of the energy between two nucleotides come from the PDNN model.
####################################################################################################

energy<-function(data,e=c(0.337130,-0.291760,0.779382,-0.222710,0.335910,-0.372490,0.827972,-0.143170,-0.394560,-0.948602,-0.020240,-0.989592,0.561440,-0.043490,0.969072,0.122600))
{
data1<-as.matrix(data[,3:dim(data)[2]])
energy_para<-matrix(0,dim(data1)[1],(dim(data1)[2]-1))
for(i in 1:dim(data1)[1])
{
for(j in 1:(dim(data1)[2]-1))
{
if(data1[i,j]=="A")
{ 
energy_para[i,j]<-switch(data1[i,j+1],"A"=e[1],"C"=e[2],"G"=e[3],"T"=e[4])
}
else if(data1[i,j]=="C")
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=e[5],"C"=e[6],"G"=e[7],"T"=e[8])
}
else if(data1[i,j]=="G")
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=e[9],"C"=e[10],"G"=e[11],"T"=e[12])
}
else
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=e[13],"C"=e[14],"G"=e[15],"T"=e[16])
}

}
}
energy_para_new<-cbind(data[,1:2],energy_para)
return(data.frame(energy_para_new))
}

###################################################################################################
# This function is from the package of mseq.
###################################################################################################

`setOriOffset` <-
function(data)
{
	uniq_cat <- unique(data$index)
	log_expr <- rep(0, length(data$index))
	
	for (k in uniq_cat)
	{
		log_expr[data$index == k] <- log(sum(data$count[data$index == k]) / sum(data$index == k))
	}

	return(log_expr)
}

#################################################################################################
# This function is used to calculate the weight of each poisitions when iterate the energy parameters
#################################################################################################

getWeight<-function(data,w)
{
data1<-as.matrix(data[,-(1:2)])
k<-length(w)
weight_matrix<-matrix(0,dim(data)[1],16)
for(i in 1:dim(data)[1])
{
#cat("iterRow:",i,"\n")
temp<-as.vector(data1[i,])
for(j in 1:16)
{
c<-sum(temp==j)
d<-sort.list(temp==j)
t<-d[-(1:(k-c))]
weight_matrix[i,j]<-sum(w[t])
}
}
new_data<-cbind(data[,1:2],weight_matrix)
return(data.frame(new_data))
}
###########################################################################################
# Use the data which comes from the result of the function of expData and get the weight 
# energy matrix in order to use the function of Energy calculating the energy.
###########################################################################################

energy_weight<-function(data)
{
data1<-as.matrix(data[,3:dim(data)[2]])
energy_para<-matrix(0,dim(data1)[1],(dim(data1)[2]-1))
for(i in 1:dim(data1)[1])
{
for(j in 1:(dim(data1)[2]-1))
{
if(data1[i,j]=="A")
{ 
energy_para[i,j]<-switch(data1[i,j+1],"A"=1,"C"=2,"G"=3,"T"=4)
}
else if(data1[i,j]=="C")
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=5,"C"=6,"G"=7,"T"=8)
}
else if(data1[i,j]=="G")
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=9,"C"=10,"G"=11,"T"=12)
}
else
{
energy_para[i,j]<-switch(data1[i,j+1],"A"=13,"C"=14,"G"=15,"T"=16)
}
}
}
energy_para_new<-cbind(data[,1:2],energy_para)
return(data.frame(energy_para_new))
}
##############################################################################################
# This function is from the package of mseq
##############################################################################################

`expData` <-
function(oriData, llen, rlen)
{
	cat("\nExpanding the surrounding sequences to get the data frame for formulas to use...\n")
	
	seq <- as.character(oriData$seq)

	run_i <- 0
	num_i <- sum(oriData$tag == 0)
	index <- numeric(num_i)
	sseq <- character(num_i * (llen + rlen))
	count <- numeric(num_i)
	for (i in 1 : length(seq))
	{
		if (oriData$tag[i] == 0)
		{
			run_i <- run_i + 1
			index[run_i] <- oriData$index[i]
			sseq[((run_i - 1) * (llen + rlen) + 1) : (run_i * (llen + rlen))] <- seq[(i - llen) : (i + rlen - 1)]
			count[run_i] <- oriData$count[i]
		}
	}

	# the line below is commented for Version 1.2
#	sseq <- factor(sseq, levels = c('T', 'A', 'C', 'G'))

	sseq <- matrix(sseq, ncol = llen + rlen, byrow = TRUE)
	data <- data.frame(index = index, count = count, sseq)

	## the three lines below are added for Version 1.2	
	for (i in 3 : (llen + rlen + 2))
	{
		data[, i] <- factor(as.character(data[, i]), levels=c('T', 'A', 'C', 'G'))
	}
	
	cname <- character(2 + llen + rlen)
	cname[1] <- "index"
	cname[2] <- "count"
	for (i in 3 : (2 + llen + rlen))
	{
		j <- i - 3 - llen
		if (j < 0)
		{
			cname[i] <- paste("pM", -j, sep = '')
		}
		else
		{
			cname[i] <- paste("p", j, sep = '')
		}
	}
	colnames(data) <- cname

	cat("number of genes =", length(unique(data$index)), "\n")
	cat("length of surrounding sequences =", dim(data)[2] - 2, "\n")
	cat("number of counts (positions) =", dim(data)[1], "\n")
	cat("total number of reads =", sum(data$count), "\n")

	return(data)
}

############################################################################################
# This function is from the package of mseq
############################################################################################
`getOriLogPref` <-
function(data, small_count = 0.5)
{
	cat("Calculating the original log preferences...\n")
	
	modi_count <- data$count
	modi_count[data$count == 0] <- small_count
	
	uniq_cat <- unique(data$index)
	log_pref <- rep(0, length(data$index))
	
	for (k in uniq_cat)
	{
		mean_k <- sum(data$count[data$index == k]) / sum(data$index == k)
		log_pref[data$index == k] <- log(modi_count[data$index == k] / mean_k)
	}
	
	return(log_pref)
}
######################################################################################################
# This function is from the package of mseq
######################################################################################################
`CVTrainIndexGene` <-
function(data, fold, seed)
{
	uniq_cat <- unique(data$index)
	train_index <- matrix(TRUE, fold, dim(data)[1])
	
	set.seed(seed)
	cv_cat <- sample(uniq_cat)
	for (j in 1 : length(uniq_cat))
	{
		row_num <- ceiling(j / floor(length(uniq_cat) / fold))
		if (row_num > fold)
		{
			row_num <- fold
		}
		train_index[row_num, data$index == cv_cat[j]] <- FALSE
	}
	
	return(train_index)
}
################################################################################################
#This function is from the package of mseq
################################################################################################
`getDev` <-
function(pred_count, real_count)
{
	dev <- sum(real_count[real_count != 0] * log(real_count[real_count != 0] / pred_count[real_count != 0]))
	dev <- dev - sum(real_count) + sum(pred_count)
	dev <- 2 * dev
	
	return(dev)
}

####################################################################################################
# This function is used to compute R-squared and train the weight and energy parameters. It first read 
# the data file. Then it uses the expData function get a csv file named data1 as an input of the energy 
# function. After using energy function, it gets an energy matrix called data11 which is a input file of 
# NewtonCV function. It gets R-squared and the coefficients after executing the NewtonCV function. In this
# paper we consider the twenty nucleotides before each poistion and twenty after it. The initial value of 
# weight is same for each poisiton and the initial value of energy parameter is from the PDNN model. The 
# output is a list which consists of R-squared, weight parameters and energy parameters.
###################################################################################################
compute_exe<-function(data,b=c(10,rep(1/39,39)),e=rep(1/16,16),mu=0.1,llen=20,rlen=20,fold = 5, seed = 281142, eps=0.1,max_iter=10,iter=1000)
{
data1<-expData(data,llen,rlen)
for(i in 1:max_iter)
{
cat("max_iter...",i,"\n")
cat("compute energy...\n")
data11<-energy(data1,e)
cat("compute R_squared...\n")
out<-NewtonCV(data11,b,mu,fold,seed,eps,iter)
coef<-out$coef
coef1<-c(coef[1],coef[-1])
cat("train energy...\n")
data12<-energy_weight(data1)
data12<-getWeight(data12,coef1[-1])
#Y<-data12$count
#X<-as.matrix(data12[,-(1:2)])
e1<-c(1,e)
#out2<-Newton(e1,X,Y,eps,iter=3)
out2<-NewtonCV_energy(data12,e1,fold,seed,eps,iter)
e1<-out2$coef
#e<-e1[-1]/max(abs(e1[-1]))
e<-e1[-1]
b<-coef
}
result<-list(R_squared=out$R_squared,weight=out$coef,energy=e1)
return(result)
}




#####################################################################################################
# This function is used to calculate the expressions. The inputs are the test data, the weight and 
# energy parameters which come from the result of compute_exe function.
#####################################################################################################

geneExpression<-function(data,b,e,llen=20,rlen=20)
{
data1<-expData(data,llen,rlen)
data11<-energy(data1,e)
unique_cat=unique(data1$index)
gene_expression=rep(0,length(unique_cat))
for (i in unique_cat)
{
temp_data=data11[data11$index==i,]
X=as.matrix(temp_data[,3:dim(temp_data)[2]]/1000000)
pred_count=rep(0,length(temp_data$count))
count=sum(temp_data$count[temp_data$index==i])
pref=exp(b[1]+1/(1+exp(X%*%b[-1])))
pred_count=count/sum(pref)*pref
gene_expression[i]=count/sum(pref)
}
return(gene_expression)
}


###########################################################################
# This function makes use of Newton algorithm for cross-validation. 
# And the input of this function are both the data after using the 
# function of expData which is provided in mseq Package and the initial 
# parameters. The output of this function is R squared and the coefficient.
###########################################################################
`NewtonCV_energy` <-
function(data,b, fold = 5, seed = 281142, eps=0.01,iter=1000)
{
	cat("Cross validation using iterative glm...\n")
	train_index <- CVTrainIndexGene(data, fold, seed)
	n<-length(b)-1
	dev <- rep(0, fold)
	null_dev <- rep(0, fold)
	coef<-rep(0,fold)
	for (j in 1 : fold)
	{
		
            cat("\nFold", j, "...\n")
            log_pref <- getOriLogPref(data[train_index[j, ], ], 0.5)
            X<-as.matrix(data[train_index[j,],3:dim(data)[2]]/1000000)
     
            log_expr<-setOriOffset(data[train_index[j,],])
            log_pref1<-log_pref-log_expr
           
            #for(i in 1:max_iter)
            #{
            #cat("max_iter:",i,"\n")
            out<-Newton_energy(b,X,log_pref1,eps,iter)
            coef<-out$coef
            #predCount<-exp(coef[1]+1/(1+exp(X%*%coef[-1])))
            #log_expr<-updateOffset(data[train_index[j,],],predCount,log_expr)
            #log_pref1<-log_pref1-log_expr
            pred_count<-getPredCount_pdnn(data[train_index[j,],],X,log_expr,coef)
            nowdev<-getDev(pred_count,data[train_index[j,],]$count)
            cat("deviance=",nowdev,"\n")
            #b<-coef
            #if(i>1&&abs((nowdev-fordev)/fordev<thrd))
            #{
            #cat("dev done...\n")
            #break
            #}
            #fordev<-nowdev
            #}
            b<-out$coef
            X1<-as.matrix(data[!train_index[j,],3:dim(data)[2]]/1000000)
            log_expr_test<-setOriOffset(data[!train_index[j,],])
            pred_count<-getPredCount_pdnn(data[!train_index[j, ], ], X1,log_expr_test,out$coef)	 	
	      dev[j] <- getDev(pred_count, data[!train_index[j, ], ]$count)
       	null_dev[j] <- getDev(getNullCount_pdnn(data[!train_index[j, ], ],X1,n), data[!train_index[j, ], ]$count)
		cat("\nFold ", j, ": null deviance = ", null_dev[j], "; deviance = ", dev[j], "; R squared = ", 1 - dev[j] / null_dev[j], "\n", sep = '')
	}
	
	cat("\n")
	cat("deviance =", dev, "\n")
	cat("null deviance =", null_dev, "\n")
	coef<-out$coef
	R_squared <- 1 - sum(dev) / sum(null_dev)
	cat("R_squared =", R_squared, "\n")
	ls<-list(R_squared=R_squared,coef=coef)
	return(ls)
}


#########################################################################################
# This function is used to optimize the least squares method to find optimal solution. 
# The inputs of this function are both weight parameters and the matrix of the data after
# using the function of expData which is provided in mseq Package.
#########################################################################################
Newton_energy<-function(b,X,Y,eps=0.01,iter=1000)
{
n<-length(b)
E<-matrix(0,n,n)
diag(E)<-1
Hessien<-matrix(0,n,n)
for(i in 1:iter)
{
cat("iter:...",i,"\n")
#hession
g1<-sum(-2*(Y-b[1]-1/(1+exp(X%*%cbind(b[-1])))))
g2<-t(cbind(2*(Y-b[1]-1/(1+exp(X%*%cbind(b[-1]))))*exp(X%*%cbind(b[-1]))/(1+exp(X%*%cbind(b[-1])))^2))%*%X
g<-t(cbind(g1,g2))
g11<-2*length(Y)
g12<-t((-2)*exp(X%*%cbind(b[-1]))/(1+exp(X%*%cbind(b[-1])))^2)%*%X
tmp<-as.vector((2*(Y-b[1]-1/(1+exp(X%*%b[-1])))*((exp(X%*%b[-1]))^2-exp(X%*%b[-1]))/(1+exp(X%*%b[-1]))^3))
tmp_matrix<-matrix(rep(tmp,ncol(X)),nrow=ncol(X),byrow=T)
tmp1<-(t(X)*tmp_matrix)%*%X
g22<-2*t(t(exp(X%*%b[-1])/(1+exp(X%*%b[-1]))^2)%*%X)%*%(t(exp(X%*%b[-1])/(1+exp(X%*%b[-1]))^2)%*%X)+tmp1
Hessien[1,1]<-g11
Hessien[1,2:n]<-g12
Hessien[2:n,1]<-g12
Hessien[-1,-1]<-g22
H<-chol(Hessien+E,pivot=TRUE)
HH<-H%*%t(H)
b1<-b
#d<-(-solve(HH)%*%g)
b<-b1-solve(HH)%*%g
#result<-search(f,1,4,b,X,Y,d)
#lemda<-ercichazhi(f,result,b,X,Y,d)
#b<-b1+lemda*d
res<-(Y-b[1]-1/(1+exp(X%*%b[-1])))
dev<-sum(res^2)
cat("dev=",dev,"\n")
if(sum((b-b1)^2)<=eps)
{
cat("Done...\n")
break
}
}
result<-list(coef=b,iter=iter,res=dev)
return(result)
}
#using coefficients to predict count
pred_count<-function(data,b,e,index,llen=20,rlen=20)
{
data1<-expData(data,llen,rlen)
data11<-energy(data1,e)
print(data11[data11$index==index,][,1:4])
log_pref <- getOriLogPref(data11[data11$index==index,], 0.5)
X<-as.matrix(data11[data11$index==index,3:dim(data11)[2]]/1000000)     
log_expr<-setOriOffset(data11[data11$index==index,])
log_pref1<-log_pref-log_expr
print(log_pref1)
print(dim(X))
log_count<-log_pref1+b[1]+1/(1+exp(sum(X%*%b[-1])))
count<-exp(log_count)*10^8
}



