R/analysis.R
In ffeng23/ADASPR: An R package to characterize the distributed affinities in a multi-ligand system
Defines functions
FitSteadyStateSPR
logFactorial
logPoly
preparePolynomialTsWithFactorial
preparePolynomialVariables
getSprPolynomialCoef
getPolynomialApproximatedSPR
predictPolynomialFitNoFactorial
predictPolynomialApproximatedSPR
predictPolynomialApproximatedSprKoff
fitPolySPR
plotFitSPR.kon
plotFitSPR.koff
fitPolySPRsOffsetted
fitPolySPRs
fitPolySPRsKoff
fitCoefficientPolySPRs
fitCoefficientPolySPRsKoff
fitSPR.kon
fitSPR.koff
fitSPR.KD
###this file encodes the data analysis function to
##characterize the affinity/rate constant distributions
###two parts: one for steady state, the other for non-steady state
###
#####===================================>
####now we can do the fitting.
#mode, indicating what model will be used
# 2, two state steady state fitting by linear regression
# 1, two state fitting by nonlinear
#'@title fit the sensorgrame steady state data based on the MultiLigand model
#'
#'@description estimate the distributin for equilibrium dissociation constant based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param x SensorgramData containing the data to be fitted
#'@param steadyStateStart numeric the starting time for steady state
#' optional and if provided, will overwrite the one in sensorgramdata
#'@param steadyStateEnd numeric the ending time for steady state
#' optional and if provided, will overwrite the ones in sensorgramdata
#'
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (\code{\link{fitSPR.kon}})
#'
#'@param fix.ligand NOT IMPLEMENT YET \cr a boolean indicating which ligand immoblization
#' model is using. \cr
#' TRUE, the fixed ligand immobilzation model. Rmax is used\cr
#' FALSE, the variable ligand immbolization model. Rligand
#' and efficiency are used. See also \code{\link{LangmuirModel-class}}
#'@param Rligand NOT IMPLEMENT YET \cr the input for the variable immobilization levels of
#' ligands on the chip. This is used by the variable ligand
#' immbolization model.See also \code{\link{LangmuirModel-class}}
#'@param auto bool used to control the plotting. TRUE by default,
#' so to wait for user in order to continue; FALSE if otherwise for
#' batch running
#'@note 1. update on 4/30/2017, add code to smoothly plot the fitting;
#' add parameter to control the plotting in order not to wait for user
#' so as to run in batch\cr
#'@return a list of parameters estimated
#'@examples
#' #------>for steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm)
#'
#' #plotting R_AB
#' plot(sData[[1]])
#'
#' #calculate the expected moments of kD distribution
# e_k<-sum(koff/kon*(Rmax/sum(Rmax)))
#' e_k2<-sum((koff/kon)^2*(Rmax/sum(Rmax)))
#' e_k3<-sum((koff/kon)^3*(Rmax/sum(Rmax)))
#' e_k4<-sum((koff/kon)^4*(Rmax/sum(Rmax)))
#'
#' #fit the steady state data
#' FitSteadyStateSPR(sData[[1]], degree=5, steadyStateStart=950,steadyStateEnd=1000)
#'
#'@export
#input:
####setMethod("FitTwoStateSPR", c("x"="SensorgramData"),#, "sampleFreq"="numeric"),
FitSteadyStateSPR<-function(x, degree=4, steadyStateStart,steadyStateEnd,
weights.matrix, weights.scale=1, weights.type="uniform",
#init.association=NULL, #init.dissociation=NULL,
#control=list(maxiter = 500,tol = 1e-2, minFactor=1/1E10)
#,trace=F,
fix.ligand=TRUE, Rligand, auto=T
)
{
#we first get data ready
if(class(x)!="SensorgramData")
{
stop("the input data is not the correct SensorgramData\n");
# return(FALSE);
}
#get the steady state window ready
ssStart<-x@steadyStateStart
ssEnd<-x@steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(!missing(steadyStateStart))
{
#cat("in hereJ:", steadyStateStart,"\n")
#cat(
ssStart<-steadyStateStart
}
if(!missing(steadyStateEnd))
ssEnd<-steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(ssStart<0||ssStart<min(x@associationData[,1]))
{
cat("******ERROR:steady state start point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<0||ssEnd>max(x@associationData[,1]))
{
cat("******ERROR:steady state end point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<=ssStart)
{
stop("*******ERROR: steady state end point is not set correctly...\n")
}
#if(mode!=1&&mode!=2)
#{
# cat("******WARNING: unknown value for fitting mode, use \"mode=1\" instead...\n");
# mode<-1;
# #stop("\n");
#}
if(!fix.ligand&&missing(Rligand))
{
stop("ERROR: Rligand can not be null when the variable ligand level is selected");
}
if(!fix.ligand&&(length(Rligand)!=length(x@analyteConcentrations)))
{
stop("ERROR: Rligand doesn't have the correct number of elements compared with analyteConcentrations");
}
#by now, we have set the window correctly
#next, let's get mean of the values to do the fitting
op<-par(ask=auto, #need to confirm for showing picture
#mfrow=c(2,2),#2x2 pictures on one plot
pty="s" #square plotting
)
plot(x)
meanSPR<-rep(0,length(x@associationData[1,])/2);
for(i in 1:length(meanSPR))
{
meanSPR[i]<-mean(x@associationData[x@associationData[,i*2-1]>=ssStart&x@associationData[,i*2-1]<=ssEnd,i*2])
lines(c(ssStart,ssEnd),c(meanSPR[i], meanSPR[i]), col=2,lwd=2)
}
#now we got the means, need to do polynomial fitting
inv_concs<-1/x@analyteConcentrations
#efficiency<-1
#now need to determine the weight matrix
scaledM<-(inv_concs)/min(inv_concs)/weights.scale
#steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(scaledM),
factorial=1/factorial(scaledM),
uniform=rep(1,length(scaledM)),
1/(scaledM+1e-5)
);
}
polyFit<-lm(formula=meanSPR~poly(inv_concs,degree=degree, raw=TRUE)#,
,weights=weights.matrix
)
mainStr<-expression(paste ( "Polynomial Regression:Req=Rmax-Rmax*", frac(1,group("[","B","]")),"E(KD)+ Rmax*",
frac(1,group("[","B","]")^2),"E(",KD^2,")-Rmax*",frac(1,group("[","B","]")^3),"+...",sep=""))
plot(inv_concs, meanSPR,main=mainStr #sub=subStr
,#frac(1,group("[","B","]")+frac(1,Rmax)"),
xlab=expression(frac(1,group("[","B","]"))),
,ylab="RUs", type="p", cex.main=0.8)
#prepare the concpoly nomial
xconcPoly<-seq(min(inv_concs),max(inv_concs),by=(max(inv_concs)-min(inv_concs))/100)
concPoly<-poly(xconcPoly,degree=degree, raw=TRUE)
#add interceptor term
concPoly<-cbind(rep(1,dim(concPoly)[1]),concPoly)
predPoly<-rep(0,dim(concPoly)[1])
for(i in 1:dim(concPoly)[1])
{
temp<-concPoly[i,]*polyFit$coefficients
temp[is.na(temp)]<-0
predPoly[i]<-sum(temp)
}
#e_RUs<-KD*(1/Rmax)*concs_inv+1/Rmax
lines(xconcPoly,predPoly, lty=2, col=2, lwd=2)
legend(max(inv_concs)/2,max(meanSPR),legend=c("observed","fitted"),lty=c(1,2),col=c(1,2),lwd=c(1,2))
#if(fix.ligand)
#{
# e_RUs<-Rmax*e_concs/(KD+e_concs)
#}else
#{
# cat("length: Rligand:",length(Rligand), ";efficiency:",length(efficiency), ";e_concs:",length(e_concs),"\n");
# e_RUs<-e_Rligand*efficiency*e_concs/(KD+e_concs)
#}
#lines(e_concs,e_RUs, lty=2, col=2, lwd=2)
#legend(min(concs),max(meanSPR),legend=c("fitted"),lty=c(2),col=c(2),lwd=c(2))
par(op);
Rmax<-polyFit$coefficients[1]
output<-polyFit$coefficients/Rmax*(-1)^c(0:(degree))
output[1]<-Rmax;
nameStr<-rep("", degree+1);
nameStr[1]<-"Rmax";
for(i in 1:degree)
{
nameStr[i+1]<-paste("moment",i,sep=" ")
}
names(output)<-nameStr;
return(output)
}#end of the function
#log-transformed factorial
#implementation:
# noninteger : get floor(n)
# negative and zero: factorial(0)=1, factorial(-n)=1
#'@title logarithem of factorial of some integer
#'
#'@description it calculate the factorial of integer in a logarithm format to avoid overflow
#'
#'@details it calculate as \cr
#' sum(log(1)+log(2)+...+log(n))
#'
#'@param n numeric or vector to be calculated for factorial values
#'
#'
#'@return logarithm transformed factorial values
#'@examples
#' logFactorial(c(1:10))
#' log(factorial(c(1:10)))
#'@export
logFactorial<-function(n)
{
if(missing(n))
{
stop("Error: missing input n!")
}
lfac<-n
for(i in 1:length(n))
{
if(n[i]<=0)
{
lfac[i]<-0
}else
{
#nInt<-floor(n)
nIntArr<-c(1:n[i])
lfac[i]<-sum(log(nIntArr))
}
}
return(lfac)
}
#equivalent to poly, but log-transformed and always with raw=TRUE
#'@title raw polynomials with log transformed
#'
#'@description prepare the raw polynomials for the input in log transformed format to avoid number overflow
#'
#'@details This is implemented as \cr
#' t'=1*log(t)+2log(t)+3log(t)...+4log(4) \cr
#' we implement this internal with logarithm to avoid numeric overflow
#'
#'@param x numeric vector the input to be expressed as polynomial
#'
#'
#'@return a matrix of raw polynomials log-transformed
#'@examples
#' logPoly(10,degree=4)
#' log(poly(10,degree=4, raw=T))
#'@export
logPoly<-function(x, degree)
{
pl<-matrix(nrow=length(x), ncol=degree)
nArr<-1:degree
for(i in 1:length(x))
{
pl[i,]<-nArr*log(x[i])
}
pl
}
#transform time t to become a polynomial with factorial n, t=1+t/1!+t^2/2!+..+t^n/n!
#this will be new input based on ts and factorial
#'@title transform the time t arrays with polynomial factors as a new input
#'
#'@description prepare the time t by scaled with factorials as a new input for polynomial approximation of exponential
#'
#'@details as an approximation of exponential, we transform the time t arrays by dividing
#' each polynomial term with factorial of degree \cr
#' t'=t/1!+t^2/2!+...+t^d/d! \cr
#' we implement this internal with logarithm to avoid numeric overflow
#'
#'@param t numeric vector time t array
#'@param degree numeric degree of the polynomial
#'
#'@return a matrix of raw polynomial transformed/divided by factorial
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
#input:
preparePolynomialTsWithFactorial<-function(ts, degree)
{
ptf<-logPoly(ts,degree)- matrix(rep(logFactorial(1:degree),length(ts)),nrow=length(ts), ncol=degree, byrow=T)
eptf<-exp(ptf)
eptf
}
preparePolynomialVariables<-function(v, degree, x)
{
xs<-matrix(rep(x),nrow=length(v),ncol=degree,byrow=T)
ptf<-poly(v,degree,raw=T)*xs# matrix(logFactorial(1:degree)),nrow=length(ts), ncol=degree, byrow=T)
#eptf<-exp(ptf)
#eptf
}
#this is the SPR re
#simluate SPR response by polynomail on time based on multiligand model
#return a array of polynomial coefficient over time t with degree
# sum(-1)^d*(Rmax*kon*c)*(kon*c+koff)^(d-1)
#for each individual conc
# kon, koff, Rmax are vectors
# conc is single numeric
# degree is total degree going up to.
# pc returns a vector of poly coef (without scaled by factorial) up to degree
#===============>
#'@title generate polynomial coefficients based on the input
#'
#'@description based on input kon, koff, and Rmax, generate the polynomial coefficients.
#'
#'@details this function is not used by data analysis, but simply used to testing and debugging
#' to compare the expected polynomial approximation value vs. real exponential simulations
#' under the multiligand model.
#'
#'@param kon numeric vector on rate array
#'@param koff numeric vector off rate array
#'@param Rmax numeric vector Rmax array
#'@param conc numeric single conc of analyte
#'@param degree numeric maximum degree of polynomial approximation
#'
#'@return a vector of raw polynomial coefficients for a multiliang model SPR
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
getSprPolynomialCoef<-function(kon, koff, Rmax, conc, degree)
{
n<-length(kon)
pc<-rep(0,degree)
#cat(pc)
for(i in 1:degree)
{
for(j in 1:n)
{
#cat("i:",i,";j:",j)
pc[i]<-pc[i]+(-1)^(i+1)*(Rmax[j]*kon[j]*conc)*(kon[j]*conc+koff[j])^(i-1)
#cat("pc:",pc,"\n")
}
}
pc
}
#get single time point, sinle conc, polynomial approxmated multiligand SPR value
#'@title get polynomial approximated multiligand SPR value as simulation
#'
#'@description based on input kon, koff, and Rmax, generate the polynomial approximated SPR values.
#'
#'@details this function is not used by data analysis, but simply used to testing and debugging
#' to compare the expected polynomial approximation value vs. real exponential simulations
#'
#'@param kon numeric vector on rate array
#'@param koff numeric vector off rate array
#'@param Rmax numeric vector Rmax array
#'@param conc numeric single conc of analyte
#'@param degree numeric maximum degree of polynomial approximation
#'
#'@return a vector of raw polynomial coefficients for a multiliang model SPR
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
getPolynomialApproximatedSPR<-function(kon, koff, Rmax, Ts, conc, degree)
{
#cat("degree:", degree, "\n")
pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree)
spptf<-Ts
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc)
}
spptf
}
#'@title predict the response based on fitted polynomail coefficients
#'
#'@description By the input x and fitted polyCoefficients, predict the responses
#' so to check the fitting quality.
#'
#'@details The polynomial approximation is prepared WITHOUT factorials scaling. Also,
#' we assuming the fitting has no intercept term. Only the coefficients for degree 1
#' and up are used to predict the response.
#'@param x numeric array the input parameter.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial
#'@return a numeric vector of predicted responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#'
#' predictFitSPR<-predictPolynomialFitNoFactorial(Ts,sfitPoly$coefficient)
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'
#'@export
predictPolynomialFitNoFactorial<-function(x, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient);
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-poly(x, degree, raw=T);
pc<-polyCoefficient;
pc[is.na(polyCoefficient)]<-0;
spptf<-x;
for(i in 1:length(x))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc);
}
spptf
}
#'@title predict the SPR response based on fitting polynomail coefficients
#'
#'@description By the input time and fitted polyCoefficients, predict the SPR responses
#' so to check the fitting quality
#'
#'@details The polynomial approximation is prepared with factorials scaling.
#'@param Ts time t array.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial based on the observed SPR responses assuming a multiligand model
#'@return a numeric vector of predicted SPR responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySPR(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'@export
predictPolynomialApproximatedSPR<-function(Ts, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient);
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree);
pc<-polyCoefficient;
pc[is.na(polyCoefficient)]<-0;
spptf<-Ts;
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc);
}
spptf
}
#'@title predict the SPR dissociation response based on fitting polynomail coefficients
#'
#'@description By the input time and fitted polyCoefficients, predict the SPR dissociation responses
#' so to check the fitting quality
#'
#'@details The polynomial approximation is prepared with scaled factorials.
#'@param Ts time t array.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial based on the observed SPR responses assuming a multiligand model
#'@return a numeric vector of predicted SPR responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySprKoff(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts+max(sData[[1]]@associationData[,1]), predictFitSPR, col=3, lty=3)
#'@export
predictPolynomialApproximatedSprKoff<-function(Ts, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient)-1;
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree);
pc<-polyCoefficient[-1];
pc[is.na(pc)]<-0;
spptf<-Ts;
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc)+polyCoefficient[1];
}
spptf
}
#the work horse for estimating on rate constant
#this one fit a polynomial with certain degree for RUs over Time
#to estimate the coefficience for each concentration of analyte
#Ts, time array
#RUs, RU array
#degree, degree
#'@title fit the non-steady state SPR response with polynomial approximation (single)
#'
#'@description This function only takes one set of SPR response data. Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on-rate constant
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param Ts numeric vector time t array
#'@param RUs numeric vector observed SPR responses
#'@param degree numeric degree of the polynomial
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'@return a vector of polynomial coefficients estimated
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@seealso \code{\link{fitPolySPRs}}
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySPR(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'@export
fitPolySPR<-function(Ts, RUs, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(Ts))
{
stop("Ts has not been specified!")
}
if(missing(RUs))
{
stop("RUs has not been specified!")
}
if(degree<1)
{
stop("degree can not be negative or zero!")
}
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
#else case meaning we using user-define the weight matrix
tf<-preparePolynomialTsWithFactorial(Ts, degree);
#get polynomial coefficient
#spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
fp<-lm(formula = RUs ~ tf+0,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
cf<-fp$coefficient;
names(cf)<-nstr;
#now do some debugging plot
pRUs<-predictPolynomialApproximatedSPR(Ts, cf)
if(debug)
{
plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts, pRUs, col=2, lty=2, lwd=2)
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
}
lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
###plot fitted spr, in this case, we want to control to plot/show only
### the good fitted part, but not the region where the fitting is crazy
##data, the sensorgram data
##fit, the polynomial fitted results
##
###To define the good part, we want to do the following
###1. at least 10% of data will be plotted
###2. difference is bigger than 2 folder (more or less)
##
#'@export
plotFitSPR.kon<-function(data, fit )
{
#now do some debugging plot
plot(data)
offFactor<-1.2 #2 fold difference
streak_to_search<-10 #ten points to check for difference
##know which analyte concentration is the biggest one.
#conc_in<-which.max(data@analyteConcentrations)
#if(debug)
#{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
#dataPoints<-data@associationData[, (i-1)*2+2]
RUs<-data@associationData[,(i-1)*2+2];
if(length(data@analyteConcentrations)==1){
pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[1]])#fit[[i]][[1]])
}else
{
pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[i]][[1]])
}
###looking for the points to stop
break_point_in<-floor(length(Ts)*0.1)
start_pt<-break_point_in
flag_break<-FALSE
#starting at least 10% of data points for plotting
for(i in start_pt:(length(Ts)-streak_to_search))
{
#flag_break<-TRUE
sub_flag_break<-TRUE
#search for a streak of points for breaking
for(j in 1:streak_to_search)
{
##ad a small value to make it non-zero
if(abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))<offFactor&&abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))>1/offFactor)
{
sub_flag_break<-FALSE
#break_point_in<-i+j
break;
}
}
if(sub_flag_break) #in this case, we have found a streak of no-good fitting points, break now
{
flag_break<-TRUE
break_point_in<-i
if(break_point_in>2)
{
break_point_in<-break_point_in-1
}
break;
}
}
if(!flag_break)
{
break_point_in<-length(Ts)
}
Ts.good<-Ts[c(1:break_point_in)]
pRUs.good<-pRUs[c(1:break_point_in)]
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts.good, pRUs.good, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#}
}
#'@export
plotFitSPR.koff<-function(data, fit )
{
#now do some debugging plot
plot(data)
offFactor<-2 #2 fold difference
streak_to_search<-10 #ten points to check for difference
##know which analyte concentration is the biggest one.
#conc_in<-which.max(data@analyteConcentrations)
#if(debug)
#{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@dissociationData[,(i-1)*2+1];
#dataPoints<-data@associationData[, (i-1)*2+2]
RUs<-data@dissociationData[,(i-1)*2+2];
#cat("after entering predict.....\n")
#flush.console();
associationT<-max(data@associationData[,(i-1)*2+1])
if(length(data@analyteConcentrations)==1){
#pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[1]])#fit[[i]][[1]])
pRUs<-predictPolynomialApproximatedSprKoff(Ts, fit[[2]]$coefficients)
}else
{
#pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[i]][[1]])
pRUs<-predictPolynomialApproximatedSprKoff(Ts, fit[[i]][[2]]$coefficients)
}
###looking for the points to stop
break_point_in<-floor(length(Ts)*0.1)
start_pt<-break_point_in
flag_break<-FALSE
#starting at least 10% of data points for plotting
for(i in start_pt:(length(Ts)-streak_to_search))
{
#flag_break<-FALSE
sub_flag_break<-TRUE
#search for a streak of points for breaking
for(j in 1:streak_to_search)
{
##ad a small value to make it non-zero
if(abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))<offFactor&&abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))>1/offFactor)
{
sub_flag_break<-FALSE
#break_point_in<-i+j
break;
}
}
if(sub_flag_break) #in this case, we have found a streak of no-good fitting points, break now
{
break_point_in<-i
if(break_point_in>2)
{
break_point_in<-break_point_in-1
}
flag_break<-TRUE
break;
}
#break_point_in
}
if(!flag_break)
{
break_point_in<-length(Ts)
}
Ts.good<-Ts[c(1:break_point_in)]
pRUs.good<-pRUs[c(1:break_point_in)]
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(associationT+Ts.good, pRUs.good, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#}
}
#'@title fit the non-steady state SPR responses with polynomial approximation （flipped)
#'
#'@description Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on rate
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.In this implementation, we also try offset/flip the response
#' units with Rmax. This way of flipping should weight the large RU points less similar to
#' off-rate constant estimation, hoping stabalize the estimation.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated, with first to be
#' Rmax (which is added to the front of the output).
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#'
#'@export
fitPolySPRsOffsetted<-function(data, Rmax, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
if(missing(Rmax))
{
stop("Rmax can not be found!!")
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
RUs<-data@associationData[,(i-1)*2+2];
RUsFlipped<-Rmax-RUs
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
fp<-lm(formula = RUsFlipped ~ tf+1,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
nstr<-c("Rmax",nstr)
cf<-fp$coefficient;
names(cf)<-nstr;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot, flipped and offsetted one.
plot(data)
if(debug)
{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
pRUs<-predictPolynomialApproximatedSPR(Ts, lst[[i]][[1]][-1])*(-1)
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts, pRUs, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#do not use plot fit spr good function yet. need to modify!!!
#plotFitSPR.kon(data, lst);
}
#lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
#'@title fit the non-steady state SPR responses with polynomial approximation
#'
#'@description Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on rate
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#'
#'@export
fitPolySPRs<-function(data, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
RUs<-data@associationData[,(i-1)*2+2];
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
fp<-lm(formula = RUs ~ tf+0,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
cf<-fp$coefficient;
names(cf)<-nstr;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot
#plot(data)
if(debug)
{
#for(i in 1:length(data@analyteConcentrations))
#{
# Ts<-data@associationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
# pRUs<-predictPolynomialApproximatedSPR(Ts, lst[[i]][[1]])
#=====>>>>>>>>>>>>check for the fitting quality and only plot the good fitting part.....
##what is a good criteria for a "good" fitting???
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
# lines(Ts, pRUs, col=2, lty=2, lwd=2)
#}
#legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
plotFitSPR.kon(data, lst);
}
#lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
#fitting the polynomial approximated dissociation phase for off rate constant
#
#'@title fit the non-steady state SPR responses for off rate constant
#'
#'@description Analyze the non-steady state SPR dissociation responses assuming the multiligand model
#' with the polynomial approximation. The distribution of off rate constant is
#' related to the levels of complex formed at the beginning of the dissocation phase.
#' In this sense, it depends on the concentration and time length of the previous
#' association phase.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better.
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@seealso \code{\link{fitPolySPR}}
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#'
#'###=======> now following the previous code, we start testing the koffs
#' m<-fitPolySPRsKoff(sData[[1]])
#'
#' #start to get the distribution from koffs
#' r0<-matrix(Rmax,nrow=length(analyteConcentrations),ncol=length(Rmax), byrow=T)
#' for(i in 1:length(analyteConcentrations))
#' {
#' r0[i,]<-Rmax*kon*analyteConcentrations[i]/(kon*analyteConcentrations[i]+koff)*(1-exp(-(kon*analyteConcentrations[i]+koff)*associationLength))
#' }
#' #expected koff moments
#' E_R0<-c(sum(r0[1,]),sum(r0[2,]),sum(r0[3,]),sum(r0[4,]),sum(r0[5,]))
#' E_koff1<-rep(0,length(analyteConcentrations))
#' E_koff2<-rep(0,length(analyteConcentrations))
#' E_koff3<-rep(0,length(analyteConcentrations))
#' E_koff4<-rep(0,length(analyteConcentrations))
#'
#' for(i in 1:length(analyteConcentrations))
#' {
#' E_koff1[i]<-sum(koff*r0[i,]/sum(r0[i,]))
#' E_koff2[i]<-sum(koff^2*r0[i,]/sum(r0[i,]))
#' E_koff3[i]<-sum(koff^3*r0[i,]/sum(r0[i,]))
#' E_koff4[i]<-sum(koff^4*r0[i,]/sum(r0[i,]))
#'
#' }
#' #fitted koff moments
#' f_E_koff1<-rep(0,length(analyteConcentrations))
#' f_E_koff2<-rep(0,length(analyteConcentrations))
#' f_E_koff3<-rep(0,length(analyteConcentrations))
#' f_E_koff4<-rep(0,length(analyteConcentrations))
#' for(i in 1:length(analyteConcentrations))
#' {
#' f_E_koff1[i]<-m[[i]][[1]][2]
#' f_E_koff2[i]<-m[[i]][[1]][3]
#' f_E_koff3[i]<-m[[i]][[1]][4]
#' f_E_koff4[i]<-m[[i]][[1]][5]
#'
#' }
#'
#'@export
fitPolySPRsKoff<-function(data, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@dissociationData[,(i-1)*2+1];
RUs<-data@dissociationData[,(i-1)*2+2];
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
#weights.matrix<-rep(1,length(steppedM))
fp<-lm(formula = RUs ~ tf,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:(degree)), sep="");
cf<-fp$coefficient*(-1)^(c(1:(degree+1))-1);
names(cf)<-c("R0",nstr);
R0<-cf[1]
cf[-1]<-cf[-1]/R0;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot
#plot(data)
if(debug)
{
#for(i in 1:length(data@analyteConcentrations))
#{
# Ts<-data@dissociationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
#cat("before entering predict..i...",i,"\n")
#flush.console();
# pRUs<-predictPolynomialApproximatedSprKoff(Ts, lst[[i]][[2]]$coefficients)
#cat("after entering predict.....\n")
#flush.console();
# associationT<-max(data@associationData[,(i-1)*2+1])
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
# lines(Ts+associationT, pRUs, col=2, lty=2, lwd=2)
#}
#legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
plotFitSPR.koff(data,lst)
}
lst
}
#this function takes the output from the function fitPolySPR and then fit a
#polynomials over conc and estimated the distribution moments of rate constants.
#
#'@title Esitmate the on rate distribution on non-steady state SPR data
#'
#'@description Characterize the on rate distribution using non-steady state SPR data. The moments
#' of on rate distribution will be estimated.
#'
#'@details The moments of the on rate distributions will be estimated based on the
#' polynominal approximation of SPR responses. The multiligand model is assumed
#' and if total Rmax is known, the true moments of the distribution will be computed.
#' Otherwise, the results are Rmax scaled moments. Practically, Rmax can be estimated
#' based on the levels of immobilized ligand assuming the binding efficiency of 100%.
#'
#'@param coefficient list object holding the estimated polynomial coefficients outputed
#' by the function \code{\link{fitPolySPRs}}
#'@param mode type of coefficient of esitmated. mode=1, regular estimation and
#' mode=2, flipped/offsetted estimation. In mode 2 estimation, the output
#' coefficient contains RMax in the beginning, but mode=1 doesn't.
#'@param analyteConcentrations numeric vector the analyte concentrations for SPR responses
#'@param Rmax numeric optional the total Rmax for ligands. Practically, this could be
#' levels of immobilized ligand when the true Rmax is not known.
#'@param debug bool used to indicate whether to plot and save the debugging
#' figures for fitting. Also write the fitting information to file.
#'
#'@return a vector of polynomial coefficients estimated
#'
#'@seealso \code{\link{fitPolySPRs}} \code{\link{fitPolySPR}}
#'@examples
#'
#' ####--------->
#' #testing fitPolySPRs
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' #5E-4, 4E-4, 3E-4, 2E-4,
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#'
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#' ###======>
#' ###following the above code, starting testing the fitting coefficient of polySPR
#' #first get the distribution of Rmax
#' E_kon<-rep(0,length(kon))
#' E_kon[1]<-sum(Rmax/sum(Rmax)*kon)
#' E_kon[2]<-sum(Rmax/sum(Rmax)*kon^2)
#' E_kon[3]<-sum(Rmax/sum(Rmax)*kon^3)
#' E_kon[4]<-sum(Rmax/sum(Rmax)*kon^4)
#' E_kon[5]<-sum(Rmax/sum(Rmax)*kon^5)
#'
#' m<-fitCoefficientPolySPRs(fpc, analyteConcentrations)
#'
#'@export
fitCoefficientPolySPRs<-function(coefficient, analyteConcentrations, Rmax,debug=F, mode=1)
{
#first go through the list from the fitted polynomial coefficients
if(class(coefficient)!="list")
{
stop("the coefficient is not set up correctly")
}
degree<-length(coefficient)
if(degree<1)
{
stop("coefficient has not set up correctly with zero degree")
}
if(mode!=1&&mode!=2)
{
stop("unknown fitting mode, please specify again!!")
}
if(missing(Rmax))
{
if(mode==1)
{
Rmax<-1;
}
#else if mode ==2 ,we don't care, we don't need
#{
# stop("rmax has not been specified.")
#}
}
#degree<-6
m<-rep(0,degree)
#for degree j
#fit lm with degree j
timeStep<-as.numeric(Sys.time())
if(debug){
fileConn<-file(paste("debug_",timeStep,".txt",sep=""), "a")
}
Rmax.fitted<-Rmax
if(mode==2){ ##mode==2 we already have Rmax, mode=1 add Rmax to the beginning
Y<-rep(0, length(coefficient))
for(i in 1:degree)
{
Y[i]<-coefficient[[i]][[1]][1]
coefficient[[i]][[1]]<-coefficient[[i]][[1]][-1];
}
lm.f<-lm(formula=Y~1)
Rmax.fitted<-lm.f$coefficients[1]
}
for(j in 1:degree)
{
Y<-rep(0, degree)
curr_degree<-j
for(i in 1:degree)
{
Y[i]<-coefficient[[i]][[1]][curr_degree]
}
conc_t<-poly(analyteConcentrations,degree=curr_degree, raw=T)
if(mode==1){
Y<-(-1)^(curr_degree-1)*Y
}
else{
Y<-(-1)^(curr_degree)*Y
}
lm.f<-lm(formula = Y ~ conc_t+0
# weights = 1/factorial(floor(x + 0.1))
#weights=weights.matrix
)
m[j]<-lm.f$coefficients[curr_degree]
#plotting
if(debug)
{
jpeg(paste("fitting_degree_",j,"_",timeStep,".jpg",sep=""))
plot(analyteConcentrations, Y, main=paste("degree=",j, sep=""),xlab="analyte conc", ylab="coef fitted")
#pYs<-predictPolynomialFitNoFactorial(analyteConcentrations, lm.f$coefficients)
x<-seq(min(analyteConcentrations), max(analyteConcentrations), by=(max(analyteConcentrations)-min(analyteConcentrations))/100)
pxYs<-predictPolynomialFitNoFactorial(x, lm.f$coefficients)
lines(x, pxYs, col=2, lty=2)
dev.off()
#save debugging information for fitting
cat("***********",format(Sys.time(), "%a %b %d %X %Y"),"\n",file=fileConn,sep="")
cat("fitting summaries for fitted coefficient SPR kon rates\n",file=fileConn,sep="")
cat("\t===>degree ",j,":\n",file=fileConn,sep="")
sink(file=fileConn, append=T)
#cat("hello\n")
print(summary(lm.f))
sink();
}
}
if(debug)
{
cat(degree," debugging plots written to ", getwd(),"\n");
cat("fitting information has been written to file at\n\t ");
cat(getwd(),"debug.txt\n")
close(fileConn)
}
if(length(Rmax.fitted)<1)
{
warning("Rmax has not been set correctly, size of zero")
Rmax.fitted<-1
}else if(length(Rmax.fitted)>degree)
{
warning("Rmax has not been set correctly, size larger than expected")
Rmax.fitted<-Rmax.fitted[1:degree]
}
cat("rmax:", Rmax.fitted,"\n")
m<-c(Rmax.fitted,m)
nstr<-c(1:degree);
nstr<-paste("moment ", nstr, sep="");
nstr<-c("Rmax", nstr)
names(m)<-nstr;
m<-m/Rmax.fitted
m
}
#this is for steady state data off rate constant.
#'@title Esitmate the off rate distribution on steady state SPR data
#'
#'@description Characterize the off rate distribution using steady state SPR data
#'(assuming all the series reach the steady state. Is this realistic/possible??).
#'The moments
#' of off rate distribution will be estimated.
#'
#'@details The moments of the off rate distributions will be estimated based on the
#' polynominal approximation of SPR responses. The multiligand model is assumed
#' and assuming total Rmax is known, the true moments of the distribution (based on Rmax) will be computed.
#'
#'@param coefficient list object holding the estimated polynomial coefficients outputed
# by the function \code{\link{fitPolySPRs}}
#'@param analyteConcentrations numeric vector the analyte concentrations for SPR responses
#'@param Rmax numeric total level of ligands
#'@param moments numeric number of moments to estimate
#'@param debug bool used to indication whether to plot and save the
#' fitting figures as well as the fitting results
#'@param degree numeric indicates the number of degree of the polynomial
#' for fitting the moments of the distribution. Without noise, a degree
#' of 2 or 3 might be good, but with noise, first degree of polynomial
#' might give out the better results since otherwise the problem of
#' of overfitting might be very pronounced.
#'@param moments numeric number of moments to be estimated. Moment up to
#' seven might be possible. otherwise, the bias might affect the results
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of moments estimated
#'
#'@seealso \code{\link{fitPolySPRsKoff}} \code{\link{fitCoefficientPolySPRs}}
#'@examples
#'
#' ####--------->
#' #testing fitPolySPRs
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' #5E-4, 4E-4, 3E-4, 2E-4,
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#'
#'###=======> now following the previous code, we start testing the koffs
#' m<-fitPolySPRsKoff(sData[[1]])
#'
#' #start to get the distribution from koffs
#' r0<-matrix(Rmax,nrow=length(analyteConcentrations),ncol=length(Rmax), byrow=T)
#' for(i in 1:length(analyteConcentrations))
#' {
#' r0[i,]<-Rmax*kon*analyteConcentrations[i]/(kon*analyteConcentrations[i]+koff)*(1-exp(-(kon*analyteConcentrations[i]+koff)*associationLength))
#' }
#' #expected koff moments
#' E_R0<-c(sum(r0[1,]),sum(r0[2,]),sum(r0[3,]),sum(r0[4,]),sum(r0[5,]))
#' E_koff1<-rep(0,length(analyteConcentrations))
#' E_koff2<-rep(0,length(analyteConcentrations))
#' E_koff3<-rep(0,length(analyteConcentrations))
#' E_koff4<-rep(0,length(analyteConcentrations))
#'
#' for(i in 1:length(analyteConcentrations))
#' {
#' E_koff1[i]<-sum(koff*r0[i,]/sum(r0[i,]))
#' E_koff2[i]<-sum(koff^2*r0[i,]/sum(r0[i,]))
#' E_koff3[i]<-sum(koff^3*r0[i,]/sum(r0[i,]))
#' E_koff4[i]<-sum(koff^4*r0[i,]/sum(r0[i,]))
#'
#' }
#' #fitted koff moments
#' f_E_koff1<-rep(0,length(analyteConcentrations))
#' f_E_koff2<-rep(0,length(analyteConcentrations))
#' f_E_koff3<-rep(0,length(analyteConcentrations))
#' f_E_koff4<-rep(0,length(analyteConcentrations))
#' for(i in 1:length(analyteConcentrations))
#' {
#' f_E_koff1[i]<-m[[i]][[1]][2]
#' f_E_koff2[i]<-m[[i]][[1]][3]
#' f_E_koff3[i]<-m[[i]][[1]][4]
#' f_E_koff4[i]<-m[[i]][[1]][5]
#' }
#'
#' ###now testing the steady state koff
#' ##based on the previous fitPolySPRsKoff
#' mss<-fitCoefficientPolySPRsKoff(m, analyteConcentrations, sum(Rmax))
#'
#' #expected koff moments
#' E_koff<-rep(0,7)
#' E_koff[1]<-sum(Rmax)
#' E_koff[2]<-sum(koff*Rmax/sum(Rmax))
#' E_koff[3]<-sum(koff^2*Rmax/sum(Rmax))
#' E_koff[4]<-sum(koff^3*Rmax/sum(Rmax))
#' E_koff[5]<-sum(koff^4*Rmax/sum(Rmax))
#' E_koff[6]<-sum(koff^5*Rmax/sum(Rmax))
#' E_koff[7]<-sum(koff^6*Rmax/sum(Rmax))
#'
#'@export
fitCoefficientPolySPRsKoff<-function(coefficient, analyteConcentrations, Rmax, debug=F, degree=3, moments=7
, weights.type="exp", weights.scale=1, weights.matrix
)
{
#first go through the list from the fitted polynomial coefficients
if(class(coefficient)!="list")
{
stop("the coefficient is not set up correctly")
}
if(degree>length(analyteConcentrations))
degree<-length(analyteConcentrations)
#mNum<-length(coefficient)
#if(degree<1)
#{
# stop("coefficient has not set up correctly with zero degree")
#}
#moments<-degree
if(length(coefficient)!=length(analyteConcentrations))
{
stop("somehting wrong, the analyte concentration and input coefficient are different in length!")
}
#degree<-6
m<-rep(0,moments)
lenConc<-length(analyteConcentrations)
#for degree 1
#fit lm with degree 1
timeStep<-as.numeric(Sys.time())
if(debug){
fileConn<-file(paste("debug_",timeStep,".txt",sep=""), "a")
}
scaledM<-1/analyteConcentrations/min(1/analyteConcentrations)/weights.scale
#steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(scaledM),
factorial=1/factorial(scaledM),
uniform=rep(1,length(scaledM)),
1/(scaledM+1e-5)
);
}
for(j in 1:(moments)) #moments is same as degree!!! so, that is how many diff conc and also how many coef can be estimated
{
Y<-rep(0, lenConc)
curr_moments<-j
for(i in 1:lenConc)
{
Y[i]<-coefficient[[i]][[2]]$coefficient[curr_moments] #here using the fit data, since the coefficient data has been change (divided by R0)
}
conc_t<-poly(1/analyteConcentrations,degree=degree, raw=T)
Y<-(-1)^(curr_moments-1)*Y
lm.f<-lm(formula = Y ~ conc_t
,weights = weights.matrix
#weights=weights.matrix
)
m[j]<-lm.f$coefficient[1]
#plotting
if(debug)
{
jpeg(paste("koff_fitting_degree_",j,"_",timeStep,".jpg",sep=""))
plot(1/analyteConcentrations, Y, main=paste("degree=",j, sep=""),xlab="1/[analyte conc]", ylab="coef fitted")
#pYs<-predictPolynomialFitNoFactorial(1/analyteConcentrations, lm.f$coefficients[-1])+lm.f$coefficients[1]
x<-seq(min(1/analyteConcentrations), max(1/analyteConcentrations), by=(max(1/analyteConcentrations)-min(1/analyteConcentrations))/100)
pxYs<-predictPolynomialFitNoFactorial(x, lm.f$coefficients[-1])+lm.f$coefficients[1]
#lines(1/analyteConcentrations, pYs, col=2, lty=2)
lines(x, pxYs, col=2, lty=2)
dev.off()
#save debugging information for fitting
cat("***********",format(Sys.time(), "%a %b %d %X %Y"),"\n",file=fileConn,sep="")
cat("fitting summaries for fitted coefficient SPR koff rates\n",file=fileConn,sep="")
cat("\t===>degree ",j,":\n",file=fileConn,sep="")
sink(file=fileConn, append=T)
#cat("hello\n")
print(summary(lm.f))
sink();
}
}
if(debug)
{
cat(degree," debugging plots written to ", getwd(),"\n");
cat("fitting information has been written to file at\n\t ");
cat(getwd(),"debug.txt\n")
close(fileConn)
}
nstr<-c(1:(moments-1));
nstr<-paste("moment ", nstr, sep="");
names(m)<-c("Rmax", nstr);
#if(!missing(Rmax))
#{
# if(length(Rmax)<1)
# {
# warning("Rmax has not been set correctly, size of zero")
# Rmax<-1
# }else if(length(Rmax)>degree)
# {
# warning("Rmax has not been set correctly, size larger than expected")
# Rmax<-Rmax[1:degree]
# }
if(missing(Rmax)){
Rmax<-m[1]
}
m[-1]<-m[-1]/Rmax
#}
m
}
####-->this is a wrapper function (or not???) to get all the workering functions together
#### to do fitting by assuming anyway the steady state of reaction.
#'@title fit the sensorgrame to estimate the on-rate distribution information
#'
#'@description estimate the raw moments of the on-rate constant distributin based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param Rmax numeric optional the total Rmax for ligands. When not provided,
#' the output moments is Rmax-scaled. An alternative way practically is to
#' specify the levels of immobilized ligand when the true Rmax is not known.
#' In this case, the output is constant scaled moments. The constant is so-called
#' efficiency of the ligand determined by the orientation as well as activity
#' of the ligand, etc.
#'@param degree numeric degree of the polynomial. Practically the number of data points
#' is big and more than the degree we can estimated. Usually a degree of 100 is probably
#' the best we can do.
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of moments estimated
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitSPR.kon(sData[[1]], Rmax=sum(Rmax),debug=TRUE);
#'
#' #first get the distribution of Rmax
#' E_kon<-rep(0,length(kon))
#' E_kon[1]<-sum(Rmax/sum(Rmax)*kon)
#' E_kon[2]<-sum(Rmax/sum(Rmax)*kon^2)
#' E_kon[3]<-sum(Rmax/sum(Rmax)*kon^3)
#' E_kon[4]<-sum(Rmax/sum(Rmax)*kon^4)
#' E_kon[5]<-sum(Rmax/sum(Rmax)*kon^5)
#'
#'
#'@return a list of parameters estimated
#'@export
# mode ==1 regular fitting mode; ==2 flipped and offsetted mode.
fitSPR.kon<-function(data, Rmax, mode=1, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
if(mode!=1&&mode!=2)
{
stop("unkown fitting mode, please specify again")
}
if(missing(Rmax))
{
if(mode==1)
{
Rmax<-1
}
else{
stop("please specify Rmax")
}
}
SPR.fpc<-NULL
if(missing(weights.matrix))
{
cat("doing fitting the SPR sensorgrams......\n")
if(mode==1){
SPR.fpc<-fitPolySPRs(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale
)
}
else{ #flip &offset mode
SPR.fpc<-fitPolySPRsOffsetted(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale,
Rmax=Rmax
)
}
}
else{
if(class(weights.matrix)!="numeric"||class(weights.matrix)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
cat("doing fitting the SPR sensorgrams......\n")
if(mode==1){
SPR.fpc<-fitPolySPRs(data, degree=degree,
weights.matrix=weights.matrix,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale
)
}
else{ #flip &offset mode
SPR.fpc<-fitPolySPRsOffsetted(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale,
Rmax=Rmax
)
}
}
#call to fit
cat("doing fitting for moments of kon.......")
m<-fitCoefficientPolySPRs(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, mode=mode)
return(m)
}#end of the function
###---->wrapper function to put together all the code to estimate the off-rate constant distribution
#'@title fit the sensorgrame to estimate the off-rate distribution information
#'
#'@description estimate the raw moments of the off-rate constant distributin based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the individual data set, which might have different quality
#' of the polynomial approximation of the exponential function. The more the quality of approximation
#' the more the number of the moments to be estimated
#' . This is different from kon estimation. Check detail of the implementation here
#' \url{http://}
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param Rmax numeric optional the total Rmax for ligands. When not provided,
#' the output moments is Rmax-scaled. An alternative way practically is to
#' specify the levels of immobilized ligand when the true Rmax is not known.
#' In this case, the output is constant scaled moments. The constant is so-called
#' efficiency of the ligand determined by the orientation as well as activity
#' of the ligand, etc.
#'@param degree.fitSPR numeric degree of the polynomial for fitting the SPR responses.
#" Practically the number of data points
#' is big and more than the degree we can estimated. Usually a degree of 100 is probably
#' the best we can do.
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'@param degree.fitMoments numeric the number of degree for the polynomial fitting of
#' moments for off-rate constants.
#'@param moments.num numeric the number of moments to estimate. Practically 6 or 7 is the best we can do.
#'
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitSPR.koff(sData[[1]], Rmax=sum(Rmax),debug=TRUE,degree.fitMoments=length(analyteConcentrations));
#'
#' #expected koff moments
#' E_koff<-rep(0,7)
#' E_koff[1]<-sum(Rmax)
#' E_koff[2]<-sum(koff*Rmax/sum(Rmax))
#' E_koff[3]<-sum(koff^2*Rmax/sum(Rmax))
#' E_koff[4]<-sum(koff^3*Rmax/sum(Rmax))
#' E_koff[5]<-sum(koff^4*Rmax/sum(Rmax))
#' E_koff[6]<-sum(koff^5*Rmax/sum(Rmax))
#' E_koff[7]<-sum(koff^6*Rmax/sum(Rmax))
#'@return A list contains two pieces of results, R0 distributed and Rmax distributed koff.
#' R0 distributed koff contains several distributions based on how many different analyte
#' concentrations are used in the experiment. Each analyte concentration generates
#' one distribution of koff. Rmax distributed koff only has one distribution of koff.
#'@export
fitSPR.koff<-function(data, Rmax, degree.fitSPR=100, weightsType.fitSPR="exp", weightsMatrix.fitSPR, debug=FALSE,
weightsStep.fitSPR=1, weightsScale.fitSPR=1,
degree.fitMoments=1, moments.num=7,
weightsType.fitMoments="exp", weightsMatrix.fitMoments, weightsScale.fitMoments=1
)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree.fitSPR<1||degree.fitMoments<1||moments.num<1)
{
stop("degree can not be negative or zero!");
}
if(missing(Rmax))
{
Rmax<-1
}
SPR.fpc<-NULL
SPR.lst<-vector("list", length=2)
if(missing(weightsMatrix.fitSPR))
{
cat("doing fitting the SPR sensorgrams......\n")
SPR.fpc<-fitPolySPRsKoff(data, degree=degree.fitSPR,
weights.type=weightsType.fitSPR,
debug=debug,
weights.step=weightsStep.fitSPR, weights.scale=weightsScale.fitSPR
)
}
else{
if(class(weightsMatrix.fitSPR)!="numeric"||class(weightsMatrix.fitSPR)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
cat("doing fitting the SPR sensorgrams......\n")
SPR.fpc<-fitPolySPRsKoff(data, degree=degree.fitSPR,
weights.matrix=weightsMatrix.fitSPR,
debug=debug#,
#weights.step=weights.step,weights.scale=weights.scale
)
}
#call to fit
cat("doing fitting for moments of koff.......")
SPR.fpm=NULL;
if(moments.num>degree.fitSPR){
moments.num<-degree.fitSPR
}
if(missing(weightsMatrix.fitSPR))
{
SPR.fpm<-fitCoefficientPolySPRsKoff(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, degree=degree.fitMoments, moments=moments.num
,weights.type=weightsType.fitMoments, weights.scale=weightsScale.fitMoments
)
}
else{
if(class(weightsMatrix.fitSPR)!="numeric"||class(weightsMatrix.fitSPR)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
SPR.fpm<-fitCoefficientPolySPRsKoff(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, degree=degree.fitMoments, moments=moments.num
,weights.matrix=weightsMatrix.fitMoments
)
}
SPR.lst<-list("R0.distributed"=SPR.fpc, "Rmax.distributed"=SPR.fpm)
return(SPR.lst)
}#end of the function
###--> this is a wrapper function to call the worker function to do the fitting
#'@title fit the sensorgrame steady state data based on the MultiLigand model
#'
#'@description estimate the distributin for equilibrium dissociation constant based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param x SensorgramData containing the data to be fitted
#'
#'@param steadyStateStart numeric the starting time for steady state
#' optional and if provided, will overwrite the one in sensorgramdata
#'@param steadyStateEnd numeric the ending time for steady state
#' optional and if provided, will overwrite the ones in sensorgramdata
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (\code{\link{fitSPR.kon}})
#'
#'@param fix.ligand NOT IMPLEMENT YET \cr a boolean indicating which ligand immoblization
#' model is using. \cr
#' TRUE, the fixed ligand immobilzation model. Rmax is used\cr
#' FALSE, the variable ligand immbolization model. Rligand
#' and efficiency are used. See also \code{\link{LangmuirModel-class}}
#'@param Rligand NOT IMPLEMENT YET \cr the input for the variable immobilization levels of
#' ligands on the chip. This is used by the variable ligand
#' immbolization model.See also \code{\link{LangmuirModel-class}}
#'@param auto bool used to control the plotting. TRUE by default,
#' so to wait for user in order to continue; FALSE if otherwise for
#' batch running
#'@return a list of parameters estimated
#'@examples
#' #------>for steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm)
#'
#' #plotting R_AB
#' plot(sData[[1]])
#'
#' #calculate the expected moments of kD distribution
#' e_k<-rep(0,length(analyteConcentrations))
#' e_k[1]<-sum(koff/kon*(Rmax/sum(Rmax)))
#' e_k[2]<-sum((koff/kon)^2*(Rmax/sum(Rmax)))
#' e_k[3]<-sum((koff/kon)^3*(Rmax/sum(Rmax)))
#' e_k[4]<-sum((koff/kon)^4*(Rmax/sum(Rmax)))
#'
#' #fit the steady state data
#' fitSPR.KD(sData[[1]], degree=length(analyteConcentrations), steadyStateStart=950,steadyStateEnd=1000)
#'
#'@export
#input:
####setMethod("FitTwoStateSPR", c("x"="SensorgramData"),#, "sampleFreq"="numeric"),
fitSPR.KD<-function(x, degree=4, steadyStateStart,steadyStateEnd,
weights.matrix, weights.scale=1, weights.type="uniform",
#init.association=NULL, #init.dissociation=NULL,
#control=list(maxiter = 500,tol = 1e-2, minFactor=1/1E10)
#,trace=F,
fix.ligand=TRUE, Rligand, auto=T
)
{
#we first get data ready
if(class(x)!="SensorgramData")
{
stop("the input data is not the correct SensorgramData\n");
# return(FALSE);
}
#get the steady state window ready
ssStart<-x@steadyStateStart
ssEnd<-x@steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(!missing(steadyStateStart))
{
#cat("in hereJ:", steadyStateStart,"\n")
#cat(
ssStart<-steadyStateStart
}
if(!missing(steadyStateEnd))
ssEnd<-steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(ssStart<0||ssStart<min(x@associationData[,1]))
{
cat("******ERROR:steady state start point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<0||ssEnd>max(x@associationData[,1]))
{
cat("******ERROR:steady state end point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<=ssStart)
{
stop("*******ERROR: steady state end point is not set correctly...\n")
}
#if(mode!=1&&mode!=2)
#{
# cat("******WARNING: unknown value for fitting mode, use \"mode=1\" instead...\n");
# mode<-1;
# #stop("\n");
#}
if(!fix.ligand&&missing(Rligand))
{
stop("ERROR: Rligand can not be null when the variable ligand level is selected");
}
if(!fix.ligand&&(length(Rligand)!=length(x@analyteConcentrations)))
{
stop("ERROR: Rligand doesn't have the correct number of elements compared with analyteConcentrations");
}
if(degree>length(x@analyteConcentrations))
{
warnings("the degree set from the input is larger than the number of analyteConcentrations!!!Automatically set it to the biggest one!",
call.=TRUE, immediate.=TRUE
)
degree<-length(x@analyteConcentrations)
}
output<-NULL
if(missing(weights.matrix))
{
output<-FitSteadyStateSPR(x,degree=degree,
steadyStateStart=steadyStateStart, steadyStateEnd=steadyStateEnd,
weights.scale=weights.scale, weights.type=weights.type,
fix.ligand=fix.ligand, Rligand=Rligand,
auto=auto
)
}
else
{
output<-FitSteadyStateSPR(x,degree=degree,
steadyStateStart=steadyStateStart, steadyStateEnd=steadyStateEnd,
weights.matrix=weights.matrix,
fix.ligand=fix.ligand, Rligand=Rligand,
auto=auto
)
}
return(output)
}#end of the function
ffeng23/ADASPR documentation built on July 13, 2019, 1:15 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 FitSteadyStateSPR logFactorial logPoly preparePolynomialTsWithFactorial preparePolynomialVariables getSprPolynomialCoef getPolynomialApproximatedSPR predictPolynomialFitNoFactorial predictPolynomialApproximatedSPR predictPolynomialApproximatedSprKoff fitPolySPR plotFitSPR.kon plotFitSPR.koff fitPolySPRsOffsetted fitPolySPRs fitPolySPRsKoff fitCoefficientPolySPRs fitCoefficientPolySPRsKoff fitSPR.kon fitSPR.koff fitSPR.KD
###this file encodes the data analysis function to
##characterize the affinity/rate constant distributions
###two parts: one for steady state, the other for non-steady state
###
#####===================================>
####now we can do the fitting.
#mode, indicating what model will be used
# 2, two state steady state fitting by linear regression
# 1, two state fitting by nonlinear
#'@title fit the sensorgrame steady state data based on the MultiLigand model
#'
#'@description estimate the distributin for equilibrium dissociation constant based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param x SensorgramData containing the data to be fitted
#'@param steadyStateStart numeric the starting time for steady state
#' optional and if provided, will overwrite the one in sensorgramdata
#'@param steadyStateEnd numeric the ending time for steady state
#' optional and if provided, will overwrite the ones in sensorgramdata
#'
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (\code{\link{fitSPR.kon}})
#'
#'@param fix.ligand NOT IMPLEMENT YET \cr a boolean indicating which ligand immoblization
#' model is using. \cr
#' TRUE, the fixed ligand immobilzation model. Rmax is used\cr
#' FALSE, the variable ligand immbolization model. Rligand
#' and efficiency are used. See also \code{\link{LangmuirModel-class}}
#'@param Rligand NOT IMPLEMENT YET \cr the input for the variable immobilization levels of
#' ligands on the chip. This is used by the variable ligand
#' immbolization model.See also \code{\link{LangmuirModel-class}}
#'@param auto bool used to control the plotting. TRUE by default,
#' so to wait for user in order to continue; FALSE if otherwise for
#' batch running
#'@note 1. update on 4/30/2017, add code to smoothly plot the fitting;
#' add parameter to control the plotting in order not to wait for user
#' so as to run in batch\cr
#'@return a list of parameters estimated
#'@examples
#' #------>for steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm)
#'
#' #plotting R_AB
#' plot(sData[[1]])
#'
#' #calculate the expected moments of kD distribution
# e_k<-sum(koff/kon*(Rmax/sum(Rmax)))
#' e_k2<-sum((koff/kon)^2*(Rmax/sum(Rmax)))
#' e_k3<-sum((koff/kon)^3*(Rmax/sum(Rmax)))
#' e_k4<-sum((koff/kon)^4*(Rmax/sum(Rmax)))
#'
#' #fit the steady state data
#' FitSteadyStateSPR(sData[[1]], degree=5, steadyStateStart=950,steadyStateEnd=1000)
#'
#'@export
#input:
####setMethod("FitTwoStateSPR", c("x"="SensorgramData"),#, "sampleFreq"="numeric"),
FitSteadyStateSPR<-function(x, degree=4, steadyStateStart,steadyStateEnd,
weights.matrix, weights.scale=1, weights.type="uniform",
#init.association=NULL, #init.dissociation=NULL,
#control=list(maxiter = 500,tol = 1e-2, minFactor=1/1E10)
#,trace=F,
fix.ligand=TRUE, Rligand, auto=T
)
{
#we first get data ready
if(class(x)!="SensorgramData")
{
stop("the input data is not the correct SensorgramData\n");
# return(FALSE);
}
#get the steady state window ready
ssStart<-x@steadyStateStart
ssEnd<-x@steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(!missing(steadyStateStart))
{
#cat("in hereJ:", steadyStateStart,"\n")
#cat(
ssStart<-steadyStateStart
}
if(!missing(steadyStateEnd))
ssEnd<-steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(ssStart<0||ssStart<min(x@associationData[,1]))
{
cat("******ERROR:steady state start point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<0||ssEnd>max(x@associationData[,1]))
{
cat("******ERROR:steady state end point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<=ssStart)
{
stop("*******ERROR: steady state end point is not set correctly...\n")
}
#if(mode!=1&&mode!=2)
#{
# cat("******WARNING: unknown value for fitting mode, use \"mode=1\" instead...\n");
# mode<-1;
# #stop("\n");
#}
if(!fix.ligand&&missing(Rligand))
{
stop("ERROR: Rligand can not be null when the variable ligand level is selected");
}
if(!fix.ligand&&(length(Rligand)!=length(x@analyteConcentrations)))
{
stop("ERROR: Rligand doesn't have the correct number of elements compared with analyteConcentrations");
}
#by now, we have set the window correctly
#next, let's get mean of the values to do the fitting
op<-par(ask=auto, #need to confirm for showing picture
#mfrow=c(2,2),#2x2 pictures on one plot
pty="s" #square plotting
)
plot(x)
meanSPR<-rep(0,length(x@associationData[1,])/2);
for(i in 1:length(meanSPR))
{
meanSPR[i]<-mean(x@associationData[x@associationData[,i*2-1]>=ssStart&x@associationData[,i*2-1]<=ssEnd,i*2])
lines(c(ssStart,ssEnd),c(meanSPR[i], meanSPR[i]), col=2,lwd=2)
}
#now we got the means, need to do polynomial fitting
inv_concs<-1/x@analyteConcentrations
#efficiency<-1
#now need to determine the weight matrix
scaledM<-(inv_concs)/min(inv_concs)/weights.scale
#steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(scaledM),
factorial=1/factorial(scaledM),
uniform=rep(1,length(scaledM)),
1/(scaledM+1e-5)
);
}
polyFit<-lm(formula=meanSPR~poly(inv_concs,degree=degree, raw=TRUE)#,
,weights=weights.matrix
)
mainStr<-expression(paste ( "Polynomial Regression:Req=Rmax-Rmax*", frac(1,group("[","B","]")),"E(KD)+ Rmax*",
frac(1,group("[","B","]")^2),"E(",KD^2,")-Rmax*",frac(1,group("[","B","]")^3),"+...",sep=""))
plot(inv_concs, meanSPR,main=mainStr #sub=subStr
,#frac(1,group("[","B","]")+frac(1,Rmax)"),
xlab=expression(frac(1,group("[","B","]"))),
,ylab="RUs", type="p", cex.main=0.8)
#prepare the concpoly nomial
xconcPoly<-seq(min(inv_concs),max(inv_concs),by=(max(inv_concs)-min(inv_concs))/100)
concPoly<-poly(xconcPoly,degree=degree, raw=TRUE)
#add interceptor term
concPoly<-cbind(rep(1,dim(concPoly)[1]),concPoly)
predPoly<-rep(0,dim(concPoly)[1])
for(i in 1:dim(concPoly)[1])
{
temp<-concPoly[i,]*polyFit$coefficients
temp[is.na(temp)]<-0
predPoly[i]<-sum(temp)
}
#e_RUs<-KD*(1/Rmax)*concs_inv+1/Rmax
lines(xconcPoly,predPoly, lty=2, col=2, lwd=2)
legend(max(inv_concs)/2,max(meanSPR),legend=c("observed","fitted"),lty=c(1,2),col=c(1,2),lwd=c(1,2))
#if(fix.ligand)
#{
# e_RUs<-Rmax*e_concs/(KD+e_concs)
#}else
#{
# cat("length: Rligand:",length(Rligand), ";efficiency:",length(efficiency), ";e_concs:",length(e_concs),"\n");
# e_RUs<-e_Rligand*efficiency*e_concs/(KD+e_concs)
#}
#lines(e_concs,e_RUs, lty=2, col=2, lwd=2)
#legend(min(concs),max(meanSPR),legend=c("fitted"),lty=c(2),col=c(2),lwd=c(2))
par(op);
Rmax<-polyFit$coefficients[1]
output<-polyFit$coefficients/Rmax*(-1)^c(0:(degree))
output[1]<-Rmax;
nameStr<-rep("", degree+1);
nameStr[1]<-"Rmax";
for(i in 1:degree)
{
nameStr[i+1]<-paste("moment",i,sep=" ")
}
names(output)<-nameStr;
return(output)
}#end of the function
#log-transformed factorial
#implementation:
# noninteger : get floor(n)
# negative and zero: factorial(0)=1, factorial(-n)=1
#'@title logarithem of factorial of some integer
#'
#'@description it calculate the factorial of integer in a logarithm format to avoid overflow
#'
#'@details it calculate as \cr
#' sum(log(1)+log(2)+...+log(n))
#'
#'@param n numeric or vector to be calculated for factorial values
#'
#'
#'@return logarithm transformed factorial values
#'@examples
#' logFactorial(c(1:10))
#' log(factorial(c(1:10)))
#'@export
logFactorial<-function(n)
{
if(missing(n))
{
stop("Error: missing input n!")
}
lfac<-n
for(i in 1:length(n))
{
if(n[i]<=0)
{
lfac[i]<-0
}else
{
#nInt<-floor(n)
nIntArr<-c(1:n[i])
lfac[i]<-sum(log(nIntArr))
}
}
return(lfac)
}
#equivalent to poly, but log-transformed and always with raw=TRUE
#'@title raw polynomials with log transformed
#'
#'@description prepare the raw polynomials for the input in log transformed format to avoid number overflow
#'
#'@details This is implemented as \cr
#' t'=1*log(t)+2log(t)+3log(t)...+4log(4) \cr
#' we implement this internal with logarithm to avoid numeric overflow
#'
#'@param x numeric vector the input to be expressed as polynomial
#'
#'
#'@return a matrix of raw polynomials log-transformed
#'@examples
#' logPoly(10,degree=4)
#' log(poly(10,degree=4, raw=T))
#'@export
logPoly<-function(x, degree)
{
pl<-matrix(nrow=length(x), ncol=degree)
nArr<-1:degree
for(i in 1:length(x))
{
pl[i,]<-nArr*log(x[i])
}
pl
}
#transform time t to become a polynomial with factorial n, t=1+t/1!+t^2/2!+..+t^n/n!
#this will be new input based on ts and factorial
#'@title transform the time t arrays with polynomial factors as a new input
#'
#'@description prepare the time t by scaled with factorials as a new input for polynomial approximation of exponential
#'
#'@details as an approximation of exponential, we transform the time t arrays by dividing
#' each polynomial term with factorial of degree \cr
#' t'=t/1!+t^2/2!+...+t^d/d! \cr
#' we implement this internal with logarithm to avoid numeric overflow
#'
#'@param t numeric vector time t array
#'@param degree numeric degree of the polynomial
#'
#'@return a matrix of raw polynomial transformed/divided by factorial
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
#input:
preparePolynomialTsWithFactorial<-function(ts, degree)
{
ptf<-logPoly(ts,degree)- matrix(rep(logFactorial(1:degree),length(ts)),nrow=length(ts), ncol=degree, byrow=T)
eptf<-exp(ptf)
eptf
}
preparePolynomialVariables<-function(v, degree, x)
{
xs<-matrix(rep(x),nrow=length(v),ncol=degree,byrow=T)
ptf<-poly(v,degree,raw=T)*xs# matrix(logFactorial(1:degree)),nrow=length(ts), ncol=degree, byrow=T)
#eptf<-exp(ptf)
#eptf
}
#this is the SPR re
#simluate SPR response by polynomail on time based on multiligand model
#return a array of polynomial coefficient over time t with degree
# sum(-1)^d*(Rmax*kon*c)*(kon*c+koff)^(d-1)
#for each individual conc
# kon, koff, Rmax are vectors
# conc is single numeric
# degree is total degree going up to.
# pc returns a vector of poly coef (without scaled by factorial) up to degree
#===============>
#'@title generate polynomial coefficients based on the input
#'
#'@description based on input kon, koff, and Rmax, generate the polynomial coefficients.
#'
#'@details this function is not used by data analysis, but simply used to testing and debugging
#' to compare the expected polynomial approximation value vs. real exponential simulations
#' under the multiligand model.
#'
#'@param kon numeric vector on rate array
#'@param koff numeric vector off rate array
#'@param Rmax numeric vector Rmax array
#'@param conc numeric single conc of analyte
#'@param degree numeric maximum degree of polynomial approximation
#'
#'@return a vector of raw polynomial coefficients for a multiliang model SPR
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
getSprPolynomialCoef<-function(kon, koff, Rmax, conc, degree)
{
n<-length(kon)
pc<-rep(0,degree)
#cat(pc)
for(i in 1:degree)
{
for(j in 1:n)
{
#cat("i:",i,";j:",j)
pc[i]<-pc[i]+(-1)^(i+1)*(Rmax[j]*kon[j]*conc)*(kon[j]*conc+koff[j])^(i-1)
#cat("pc:",pc,"\n")
}
}
pc
}
#get single time point, sinle conc, polynomial approxmated multiligand SPR value
#'@title get polynomial approximated multiligand SPR value as simulation
#'
#'@description based on input kon, koff, and Rmax, generate the polynomial approximated SPR values.
#'
#'@details this function is not used by data analysis, but simply used to testing and debugging
#' to compare the expected polynomial approximation value vs. real exponential simulations
#'
#'@param kon numeric vector on rate array
#'@param koff numeric vector off rate array
#'@param Rmax numeric vector Rmax array
#'@param conc numeric single conc of analyte
#'@param degree numeric maximum degree of polynomial approximation
#'
#'@return a vector of raw polynomial coefficients for a multiliang model SPR
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' Ts<-seq(0,100,by=0.1)
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #first testing transform ts with factorial
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#'
#' #get polynomial coefficient
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
#'
#' #get SPR polynomial approximated
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#'@export
getPolynomialApproximatedSPR<-function(kon, koff, Rmax, Ts, conc, degree)
{
#cat("degree:", degree, "\n")
pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree)
spptf<-Ts
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc)
}
spptf
}
#'@title predict the response based on fitted polynomail coefficients
#'
#'@description By the input x and fitted polyCoefficients, predict the responses
#' so to check the fitting quality.
#'
#'@details The polynomial approximation is prepared WITHOUT factorials scaling. Also,
#' we assuming the fitting has no intercept term. Only the coefficients for degree 1
#' and up are used to predict the response.
#'@param x numeric array the input parameter.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial
#'@return a numeric vector of predicted responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#'
#' predictFitSPR<-predictPolynomialFitNoFactorial(Ts,sfitPoly$coefficient)
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'
#'@export
predictPolynomialFitNoFactorial<-function(x, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient);
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-poly(x, degree, raw=T);
pc<-polyCoefficient;
pc[is.na(polyCoefficient)]<-0;
spptf<-x;
for(i in 1:length(x))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc);
}
spptf
}
#'@title predict the SPR response based on fitting polynomail coefficients
#'
#'@description By the input time and fitted polyCoefficients, predict the SPR responses
#' so to check the fitting quality
#'
#'@details The polynomial approximation is prepared with factorials scaling.
#'@param Ts time t array.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial based on the observed SPR responses assuming a multiligand model
#'@return a numeric vector of predicted SPR responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySPR(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'@export
predictPolynomialApproximatedSPR<-function(Ts, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient);
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree);
pc<-polyCoefficient;
pc[is.na(polyCoefficient)]<-0;
spptf<-Ts;
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc);
}
spptf
}
#'@title predict the SPR dissociation response based on fitting polynomail coefficients
#'
#'@description By the input time and fitted polyCoefficients, predict the SPR dissociation responses
#' so to check the fitting quality
#'
#'@details The polynomial approximation is prepared with scaled factorials.
#'@param Ts time t array.
#'@param polyCoefficient numeric vector that output from the linear regression
#' of polynomial based on the observed SPR responses assuming a multiligand model
#'@return a numeric vector of predicted SPR responses
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySprKoff(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts+max(sData[[1]]@associationData[,1]), predictFitSPR, col=3, lty=3)
#'@export
predictPolynomialApproximatedSprKoff<-function(Ts, polyCoefficient)
{
#cat("degree:", degree, "\n")
degree<-length(polyCoefficient)-1;
#pc<-getSprPolynomialCoef(kon, koff, Rmax, conc, degree)
ptf<-preparePolynomialTsWithFactorial(Ts, degree);
pc<-polyCoefficient[-1];
pc[is.na(pc)]<-0;
spptf<-Ts;
for(i in 1:length(Ts))
{#cat("i:pc:", pc, "\n")
#cat("\tptf[i]:", ptf[i,]*pc, "\n")
#cat("\tptf*pc", ptf[i,]*pc, "\n")
spptf[i]<-sum(ptf[i,]*pc)+polyCoefficient[1];
}
spptf
}
#the work horse for estimating on rate constant
#this one fit a polynomial with certain degree for RUs over Time
#to estimate the coefficience for each concentration of analyte
#Ts, time array
#RUs, RU array
#degree, degree
#'@title fit the non-steady state SPR response with polynomial approximation (single)
#'
#'@description This function only takes one set of SPR response data. Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on-rate constant
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param Ts numeric vector time t array
#'@param RUs numeric vector observed SPR responses
#'@param degree numeric degree of the polynomial
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'@return a vector of polynomial coefficients estimated
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@seealso \code{\link{fitPolySPRs}}
#'@examples
#' #start testing the polynomial
#' kon<-c(3e3)
#' koff<-c(5E-2)#,2E-2)
#' Rmax<-c(60)
#' conc<-1E-4
#' degree<-110
#' associationLength<-100
#' dissociationLength<-100
#'
#' #now start doing simulation with multiligand model, with a single component
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=conc,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' Ts<-sData[[1]]@associationData[,1]
#' tf<-preparePolynomialTsWithFactorial(Ts, degree);
#' paSPR<-getPolynomialApproximatedSPR(kon, koff, Rmax,Ts, conc, degree)
#'
#' #plot them together
#' plot(sData[[1]])
#'
#' lines(Ts,paSPR, col=2, lty=2)
#' legend(50,40, c("simulation by multiLigand model","polynomial with degree of 110"), col=c(1,2), pch=c(-1,-1), lty=c(1,2))
#'
#' sfitPoly<-lm(formula = sData[[1]]@associationData[,2] ~ tf +0
#' , weights=exp(-1*Ts)
#' #, weights = 1/factorial(floor(x + 0.1))
#' )
#' #fitting
#' fpc<-fitPolySPR(Ts, sData[[1]]@associationData[,2])
#' predictFitSPR<-predictPolynomialApproximatedSPR(Ts,fpc[["lm.coefficient"]])
#' lines(Ts, predictFitSPR, col=3, lty=3)
#'@export
fitPolySPR<-function(Ts, RUs, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(Ts))
{
stop("Ts has not been specified!")
}
if(missing(RUs))
{
stop("RUs has not been specified!")
}
if(degree<1)
{
stop("degree can not be negative or zero!")
}
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
#else case meaning we using user-define the weight matrix
tf<-preparePolynomialTsWithFactorial(Ts, degree);
#get polynomial coefficient
#spc<-getSprPolynomialCoef(kon,koff, Rmax, conc, degree)
fp<-lm(formula = RUs ~ tf+0,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
cf<-fp$coefficient;
names(cf)<-nstr;
#now do some debugging plot
pRUs<-predictPolynomialApproximatedSPR(Ts, cf)
if(debug)
{
plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts, pRUs, col=2, lty=2, lwd=2)
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
}
lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
###plot fitted spr, in this case, we want to control to plot/show only
### the good fitted part, but not the region where the fitting is crazy
##data, the sensorgram data
##fit, the polynomial fitted results
##
###To define the good part, we want to do the following
###1. at least 10% of data will be plotted
###2. difference is bigger than 2 folder (more or less)
##
#'@export
plotFitSPR.kon<-function(data, fit )
{
#now do some debugging plot
plot(data)
offFactor<-1.2 #2 fold difference
streak_to_search<-10 #ten points to check for difference
##know which analyte concentration is the biggest one.
#conc_in<-which.max(data@analyteConcentrations)
#if(debug)
#{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
#dataPoints<-data@associationData[, (i-1)*2+2]
RUs<-data@associationData[,(i-1)*2+2];
if(length(data@analyteConcentrations)==1){
pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[1]])#fit[[i]][[1]])
}else
{
pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[i]][[1]])
}
###looking for the points to stop
break_point_in<-floor(length(Ts)*0.1)
start_pt<-break_point_in
flag_break<-FALSE
#starting at least 10% of data points for plotting
for(i in start_pt:(length(Ts)-streak_to_search))
{
#flag_break<-TRUE
sub_flag_break<-TRUE
#search for a streak of points for breaking
for(j in 1:streak_to_search)
{
##ad a small value to make it non-zero
if(abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))<offFactor&&abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))>1/offFactor)
{
sub_flag_break<-FALSE
#break_point_in<-i+j
break;
}
}
if(sub_flag_break) #in this case, we have found a streak of no-good fitting points, break now
{
flag_break<-TRUE
break_point_in<-i
if(break_point_in>2)
{
break_point_in<-break_point_in-1
}
break;
}
}
if(!flag_break)
{
break_point_in<-length(Ts)
}
Ts.good<-Ts[c(1:break_point_in)]
pRUs.good<-pRUs[c(1:break_point_in)]
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts.good, pRUs.good, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#}
}
#'@export
plotFitSPR.koff<-function(data, fit )
{
#now do some debugging plot
plot(data)
offFactor<-2 #2 fold difference
streak_to_search<-10 #ten points to check for difference
##know which analyte concentration is the biggest one.
#conc_in<-which.max(data@analyteConcentrations)
#if(debug)
#{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@dissociationData[,(i-1)*2+1];
#dataPoints<-data@associationData[, (i-1)*2+2]
RUs<-data@dissociationData[,(i-1)*2+2];
#cat("after entering predict.....\n")
#flush.console();
associationT<-max(data@associationData[,(i-1)*2+1])
if(length(data@analyteConcentrations)==1){
#pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[1]])#fit[[i]][[1]])
pRUs<-predictPolynomialApproximatedSprKoff(Ts, fit[[2]]$coefficients)
}else
{
#pRUs<-predictPolynomialApproximatedSPR(Ts, fit[[i]][[1]])
pRUs<-predictPolynomialApproximatedSprKoff(Ts, fit[[i]][[2]]$coefficients)
}
###looking for the points to stop
break_point_in<-floor(length(Ts)*0.1)
start_pt<-break_point_in
flag_break<-FALSE
#starting at least 10% of data points for plotting
for(i in start_pt:(length(Ts)-streak_to_search))
{
#flag_break<-FALSE
sub_flag_break<-TRUE
#search for a streak of points for breaking
for(j in 1:streak_to_search)
{
##ad a small value to make it non-zero
if(abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))<offFactor&&abs((pRUs[(i+j)]+1E-20)/(RUs[i+j]+1E-20))>1/offFactor)
{
sub_flag_break<-FALSE
#break_point_in<-i+j
break;
}
}
if(sub_flag_break) #in this case, we have found a streak of no-good fitting points, break now
{
break_point_in<-i
if(break_point_in>2)
{
break_point_in<-break_point_in-1
}
flag_break<-TRUE
break;
}
#break_point_in
}
if(!flag_break)
{
break_point_in<-length(Ts)
}
Ts.good<-Ts[c(1:break_point_in)]
pRUs.good<-pRUs[c(1:break_point_in)]
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(associationT+Ts.good, pRUs.good, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#}
}
#'@title fit the non-steady state SPR responses with polynomial approximation （flipped)
#'
#'@description Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on rate
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.In this implementation, we also try offset/flip the response
#' units with Rmax. This way of flipping should weight the large RU points less similar to
#' off-rate constant estimation, hoping stabalize the estimation.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated, with first to be
#' Rmax (which is added to the front of the output).
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#'
#'@export
fitPolySPRsOffsetted<-function(data, Rmax, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
if(missing(Rmax))
{
stop("Rmax can not be found!!")
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
RUs<-data@associationData[,(i-1)*2+2];
RUsFlipped<-Rmax-RUs
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
fp<-lm(formula = RUsFlipped ~ tf+1,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
nstr<-c("Rmax",nstr)
cf<-fp$coefficient;
names(cf)<-nstr;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot, flipped and offsetted one.
plot(data)
if(debug)
{
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
pRUs<-predictPolynomialApproximatedSPR(Ts, lst[[i]][[1]][-1])*(-1)
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
lines(Ts, pRUs, col=2, lty=2, lwd=2)
}
legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
#do not use plot fit spr good function yet. need to modify!!!
#plotFitSPR.kon(data, lst);
}
#lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
#'@title fit the non-steady state SPR responses with polynomial approximation
#'
#'@description Analyze the non-steady state SPR responses assuming the multiligand model
#' with the polynomial approximation. This is the first step in order to get on rate
#' distribution.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#'
#'@export
fitPolySPRs<-function(data, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@associationData[,(i-1)*2+1];
RUs<-data@associationData[,(i-1)*2+2];
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
fp<-lm(formula = RUs ~ tf+0,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:degree), sep="");
cf<-fp$coefficient;
names(cf)<-nstr;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot
#plot(data)
if(debug)
{
#for(i in 1:length(data@analyteConcentrations))
#{
# Ts<-data@associationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
# pRUs<-predictPolynomialApproximatedSPR(Ts, lst[[i]][[1]])
#=====>>>>>>>>>>>>check for the fitting quality and only plot the good fitting part.....
##what is a good criteria for a "good" fitting???
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
# lines(Ts, pRUs, col=2, lty=2, lwd=2)
#}
#legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
plotFitSPR.kon(data, lst);
}
#lst<-list(lm.coefficient=cf,lm.summary=summary(fp));
lst
}
#fitting the polynomial approximated dissociation phase for off rate constant
#
#'@title fit the non-steady state SPR responses for off rate constant
#'
#'@description Analyze the non-steady state SPR dissociation responses assuming the multiligand model
#' with the polynomial approximation. The distribution of off rate constant is
#' related to the levels of complex formed at the beginning of the dissocation phase.
#' In this sense, it depends on the concentration and time length of the previous
#' association phase.
#'
#'@details The SPR responses that is assuming a sum of exponentials is approximated
#' by polynomials. In this step, we fit the responses to estimate the
#' coefficients of approximated polynomial.
#' The implementation is to prepare the time t arrays as polymonials scaled
#' with factorials and then we fit linearly the data to estimate the coefficients
#' for each degree of polynomials. In order to make the estimation accurate,
#' an exponential weight matrix is used to lower the contribution of large time t values
#' , at which the approximation could potentially goes off.
#' A factorial weight matrix could be used too. But empiricially, an exponential one is better.
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param degree numeric degree of the polynomial
#'
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of polynomial coefficients estimated
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@seealso \code{\link{fitPolySPR}}
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#'
#'###=======> now following the previous code, we start testing the koffs
#' m<-fitPolySPRsKoff(sData[[1]])
#'
#' #start to get the distribution from koffs
#' r0<-matrix(Rmax,nrow=length(analyteConcentrations),ncol=length(Rmax), byrow=T)
#' for(i in 1:length(analyteConcentrations))
#' {
#' r0[i,]<-Rmax*kon*analyteConcentrations[i]/(kon*analyteConcentrations[i]+koff)*(1-exp(-(kon*analyteConcentrations[i]+koff)*associationLength))
#' }
#' #expected koff moments
#' E_R0<-c(sum(r0[1,]),sum(r0[2,]),sum(r0[3,]),sum(r0[4,]),sum(r0[5,]))
#' E_koff1<-rep(0,length(analyteConcentrations))
#' E_koff2<-rep(0,length(analyteConcentrations))
#' E_koff3<-rep(0,length(analyteConcentrations))
#' E_koff4<-rep(0,length(analyteConcentrations))
#'
#' for(i in 1:length(analyteConcentrations))
#' {
#' E_koff1[i]<-sum(koff*r0[i,]/sum(r0[i,]))
#' E_koff2[i]<-sum(koff^2*r0[i,]/sum(r0[i,]))
#' E_koff3[i]<-sum(koff^3*r0[i,]/sum(r0[i,]))
#' E_koff4[i]<-sum(koff^4*r0[i,]/sum(r0[i,]))
#'
#' }
#' #fitted koff moments
#' f_E_koff1<-rep(0,length(analyteConcentrations))
#' f_E_koff2<-rep(0,length(analyteConcentrations))
#' f_E_koff3<-rep(0,length(analyteConcentrations))
#' f_E_koff4<-rep(0,length(analyteConcentrations))
#' for(i in 1:length(analyteConcentrations))
#' {
#' f_E_koff1[i]<-m[[i]][[1]][2]
#' f_E_koff2[i]<-m[[i]][[1]][3]
#' f_E_koff3[i]<-m[[i]][[1]][4]
#' f_E_koff4[i]<-m[[i]][[1]][5]
#'
#' }
#'
#'@export
fitPolySPRsKoff<-function(data, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
#else case meaning we using user-define the weight matrix
lst<-vector("list", length(data@analyteConcentrations))
#now, we need to go through each analyte conc and fit for polynomial coefficients
for(i in 1:length(data@analyteConcentrations))
{
Ts<-data@dissociationData[,(i-1)*2+1];
RUs<-data@dissociationData[,(i-1)*2+2];
tf<-preparePolynomialTsWithFactorial(Ts, degree);
scaledM<-Ts/weights.scale
steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(steppedM),
factorial=1/factorial(steppedM),
uniform=rep(1,length(steppedM)),
1/(steppedM+1e-5)
);
}
#weights.matrix<-rep(1,length(steppedM))
fp<-lm(formula = RUs ~ tf,
# weights = 1/factorial(floor(x + 0.1))
weights=weights.matrix
);
nstr<-rep("coef", degree);
nstr<-paste(nstr, c(1:(degree)), sep="");
cf<-fp$coefficient*(-1)^(c(1:(degree+1))-1);
names(cf)<-c("R0",nstr);
R0<-cf[1]
cf[-1]<-cf[-1]/R0;
lst[[i]]<-list(lm.coef=cf, lm.fit=fp)
}
#now do some debugging plot
#plot(data)
if(debug)
{
#for(i in 1:length(data@analyteConcentrations))
#{
# Ts<-data@dissociationData[,(i-1)*2+1];
#RUs<-data@associationData[,(i-1)*2+2];
#cat("before entering predict..i...",i,"\n")
#flush.console();
# pRUs<-predictPolynomialApproximatedSprKoff(Ts, lst[[i]][[2]]$coefficients)
#cat("after entering predict.....\n")
#flush.console();
# associationT<-max(data@associationData[,(i-1)*2+1])
#plot(Ts, RUs, main="fitted SPR response by polynomials",type="l",xlab="time(sec)", ylab="RUs")
# lines(Ts+associationT, pRUs, col=2, lty=2, lwd=2)
#}
#legend(Ts[length(Ts)]*0.5, RUs[length(RUs)]*0.8, c("fitting SPR", "Observed SPR"), lty=c(2,1),col=c(2,1), lwd=c(2,1))
plotFitSPR.koff(data,lst)
}
lst
}
#this function takes the output from the function fitPolySPR and then fit a
#polynomials over conc and estimated the distribution moments of rate constants.
#
#'@title Esitmate the on rate distribution on non-steady state SPR data
#'
#'@description Characterize the on rate distribution using non-steady state SPR data. The moments
#' of on rate distribution will be estimated.
#'
#'@details The moments of the on rate distributions will be estimated based on the
#' polynominal approximation of SPR responses. The multiligand model is assumed
#' and if total Rmax is known, the true moments of the distribution will be computed.
#' Otherwise, the results are Rmax scaled moments. Practically, Rmax can be estimated
#' based on the levels of immobilized ligand assuming the binding efficiency of 100%.
#'
#'@param coefficient list object holding the estimated polynomial coefficients outputed
#' by the function \code{\link{fitPolySPRs}}
#'@param mode type of coefficient of esitmated. mode=1, regular estimation and
#' mode=2, flipped/offsetted estimation. In mode 2 estimation, the output
#' coefficient contains RMax in the beginning, but mode=1 doesn't.
#'@param analyteConcentrations numeric vector the analyte concentrations for SPR responses
#'@param Rmax numeric optional the total Rmax for ligands. Practically, this could be
#' levels of immobilized ligand when the true Rmax is not known.
#'@param debug bool used to indicate whether to plot and save the debugging
#' figures for fitting. Also write the fitting information to file.
#'
#'@return a vector of polynomial coefficients estimated
#'
#'@seealso \code{\link{fitPolySPRs}} \code{\link{fitPolySPR}}
#'@examples
#'
#' ####--------->
#' #testing fitPolySPRs
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' #5E-4, 4E-4, 3E-4, 2E-4,
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitPolySPRs(sData[[1]]);
#'
#' spc<-getSprPolynomialCoef(kon,koff, Rmax, analyteConcentrations[1], degree=100)
#'
#' ###======>
#' ###following the above code, starting testing the fitting coefficient of polySPR
#' #first get the distribution of Rmax
#' E_kon<-rep(0,length(kon))
#' E_kon[1]<-sum(Rmax/sum(Rmax)*kon)
#' E_kon[2]<-sum(Rmax/sum(Rmax)*kon^2)
#' E_kon[3]<-sum(Rmax/sum(Rmax)*kon^3)
#' E_kon[4]<-sum(Rmax/sum(Rmax)*kon^4)
#' E_kon[5]<-sum(Rmax/sum(Rmax)*kon^5)
#'
#' m<-fitCoefficientPolySPRs(fpc, analyteConcentrations)
#'
#'@export
fitCoefficientPolySPRs<-function(coefficient, analyteConcentrations, Rmax,debug=F, mode=1)
{
#first go through the list from the fitted polynomial coefficients
if(class(coefficient)!="list")
{
stop("the coefficient is not set up correctly")
}
degree<-length(coefficient)
if(degree<1)
{
stop("coefficient has not set up correctly with zero degree")
}
if(mode!=1&&mode!=2)
{
stop("unknown fitting mode, please specify again!!")
}
if(missing(Rmax))
{
if(mode==1)
{
Rmax<-1;
}
#else if mode ==2 ,we don't care, we don't need
#{
# stop("rmax has not been specified.")
#}
}
#degree<-6
m<-rep(0,degree)
#for degree j
#fit lm with degree j
timeStep<-as.numeric(Sys.time())
if(debug){
fileConn<-file(paste("debug_",timeStep,".txt",sep=""), "a")
}
Rmax.fitted<-Rmax
if(mode==2){ ##mode==2 we already have Rmax, mode=1 add Rmax to the beginning
Y<-rep(0, length(coefficient))
for(i in 1:degree)
{
Y[i]<-coefficient[[i]][[1]][1]
coefficient[[i]][[1]]<-coefficient[[i]][[1]][-1];
}
lm.f<-lm(formula=Y~1)
Rmax.fitted<-lm.f$coefficients[1]
}
for(j in 1:degree)
{
Y<-rep(0, degree)
curr_degree<-j
for(i in 1:degree)
{
Y[i]<-coefficient[[i]][[1]][curr_degree]
}
conc_t<-poly(analyteConcentrations,degree=curr_degree, raw=T)
if(mode==1){
Y<-(-1)^(curr_degree-1)*Y
}
else{
Y<-(-1)^(curr_degree)*Y
}
lm.f<-lm(formula = Y ~ conc_t+0
# weights = 1/factorial(floor(x + 0.1))
#weights=weights.matrix
)
m[j]<-lm.f$coefficients[curr_degree]
#plotting
if(debug)
{
jpeg(paste("fitting_degree_",j,"_",timeStep,".jpg",sep=""))
plot(analyteConcentrations, Y, main=paste("degree=",j, sep=""),xlab="analyte conc", ylab="coef fitted")
#pYs<-predictPolynomialFitNoFactorial(analyteConcentrations, lm.f$coefficients)
x<-seq(min(analyteConcentrations), max(analyteConcentrations), by=(max(analyteConcentrations)-min(analyteConcentrations))/100)
pxYs<-predictPolynomialFitNoFactorial(x, lm.f$coefficients)
lines(x, pxYs, col=2, lty=2)
dev.off()
#save debugging information for fitting
cat("***********",format(Sys.time(), "%a %b %d %X %Y"),"\n",file=fileConn,sep="")
cat("fitting summaries for fitted coefficient SPR kon rates\n",file=fileConn,sep="")
cat("\t===>degree ",j,":\n",file=fileConn,sep="")
sink(file=fileConn, append=T)
#cat("hello\n")
print(summary(lm.f))
sink();
}
}
if(debug)
{
cat(degree," debugging plots written to ", getwd(),"\n");
cat("fitting information has been written to file at\n\t ");
cat(getwd(),"debug.txt\n")
close(fileConn)
}
if(length(Rmax.fitted)<1)
{
warning("Rmax has not been set correctly, size of zero")
Rmax.fitted<-1
}else if(length(Rmax.fitted)>degree)
{
warning("Rmax has not been set correctly, size larger than expected")
Rmax.fitted<-Rmax.fitted[1:degree]
}
cat("rmax:", Rmax.fitted,"\n")
m<-c(Rmax.fitted,m)
nstr<-c(1:degree);
nstr<-paste("moment ", nstr, sep="");
nstr<-c("Rmax", nstr)
names(m)<-nstr;
m<-m/Rmax.fitted
m
}
#this is for steady state data off rate constant.
#'@title Esitmate the off rate distribution on steady state SPR data
#'
#'@description Characterize the off rate distribution using steady state SPR data
#'(assuming all the series reach the steady state. Is this realistic/possible??).
#'The moments
#' of off rate distribution will be estimated.
#'
#'@details The moments of the off rate distributions will be estimated based on the
#' polynominal approximation of SPR responses. The multiligand model is assumed
#' and assuming total Rmax is known, the true moments of the distribution (based on Rmax) will be computed.
#'
#'@param coefficient list object holding the estimated polynomial coefficients outputed
# by the function \code{\link{fitPolySPRs}}
#'@param analyteConcentrations numeric vector the analyte concentrations for SPR responses
#'@param Rmax numeric total level of ligands
#'@param moments numeric number of moments to estimate
#'@param debug bool used to indication whether to plot and save the
#' fitting figures as well as the fitting results
#'@param degree numeric indicates the number of degree of the polynomial
#' for fitting the moments of the distribution. Without noise, a degree
#' of 2 or 3 might be good, but with noise, first degree of polynomial
#' might give out the better results since otherwise the problem of
#' of overfitting might be very pronounced.
#'@param moments numeric number of moments to be estimated. Moment up to
#' seven might be possible. otherwise, the bias might affect the results
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of moments estimated
#'
#'@seealso \code{\link{fitPolySPRsKoff}} \code{\link{fitCoefficientPolySPRs}}
#'@examples
#'
#' ####--------->
#' #testing fitPolySPRs
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' #5E-4, 4E-4, 3E-4, 2E-4,
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#'
#'###=======> now following the previous code, we start testing the koffs
#' m<-fitPolySPRsKoff(sData[[1]])
#'
#' #start to get the distribution from koffs
#' r0<-matrix(Rmax,nrow=length(analyteConcentrations),ncol=length(Rmax), byrow=T)
#' for(i in 1:length(analyteConcentrations))
#' {
#' r0[i,]<-Rmax*kon*analyteConcentrations[i]/(kon*analyteConcentrations[i]+koff)*(1-exp(-(kon*analyteConcentrations[i]+koff)*associationLength))
#' }
#' #expected koff moments
#' E_R0<-c(sum(r0[1,]),sum(r0[2,]),sum(r0[3,]),sum(r0[4,]),sum(r0[5,]))
#' E_koff1<-rep(0,length(analyteConcentrations))
#' E_koff2<-rep(0,length(analyteConcentrations))
#' E_koff3<-rep(0,length(analyteConcentrations))
#' E_koff4<-rep(0,length(analyteConcentrations))
#'
#' for(i in 1:length(analyteConcentrations))
#' {
#' E_koff1[i]<-sum(koff*r0[i,]/sum(r0[i,]))
#' E_koff2[i]<-sum(koff^2*r0[i,]/sum(r0[i,]))
#' E_koff3[i]<-sum(koff^3*r0[i,]/sum(r0[i,]))
#' E_koff4[i]<-sum(koff^4*r0[i,]/sum(r0[i,]))
#'
#' }
#' #fitted koff moments
#' f_E_koff1<-rep(0,length(analyteConcentrations))
#' f_E_koff2<-rep(0,length(analyteConcentrations))
#' f_E_koff3<-rep(0,length(analyteConcentrations))
#' f_E_koff4<-rep(0,length(analyteConcentrations))
#' for(i in 1:length(analyteConcentrations))
#' {
#' f_E_koff1[i]<-m[[i]][[1]][2]
#' f_E_koff2[i]<-m[[i]][[1]][3]
#' f_E_koff3[i]<-m[[i]][[1]][4]
#' f_E_koff4[i]<-m[[i]][[1]][5]
#' }
#'
#' ###now testing the steady state koff
#' ##based on the previous fitPolySPRsKoff
#' mss<-fitCoefficientPolySPRsKoff(m, analyteConcentrations, sum(Rmax))
#'
#' #expected koff moments
#' E_koff<-rep(0,7)
#' E_koff[1]<-sum(Rmax)
#' E_koff[2]<-sum(koff*Rmax/sum(Rmax))
#' E_koff[3]<-sum(koff^2*Rmax/sum(Rmax))
#' E_koff[4]<-sum(koff^3*Rmax/sum(Rmax))
#' E_koff[5]<-sum(koff^4*Rmax/sum(Rmax))
#' E_koff[6]<-sum(koff^5*Rmax/sum(Rmax))
#' E_koff[7]<-sum(koff^6*Rmax/sum(Rmax))
#'
#'@export
fitCoefficientPolySPRsKoff<-function(coefficient, analyteConcentrations, Rmax, debug=F, degree=3, moments=7
, weights.type="exp", weights.scale=1, weights.matrix
)
{
#first go through the list from the fitted polynomial coefficients
if(class(coefficient)!="list")
{
stop("the coefficient is not set up correctly")
}
if(degree>length(analyteConcentrations))
degree<-length(analyteConcentrations)
#mNum<-length(coefficient)
#if(degree<1)
#{
# stop("coefficient has not set up correctly with zero degree")
#}
#moments<-degree
if(length(coefficient)!=length(analyteConcentrations))
{
stop("somehting wrong, the analyte concentration and input coefficient are different in length!")
}
#degree<-6
m<-rep(0,moments)
lenConc<-length(analyteConcentrations)
#for degree 1
#fit lm with degree 1
timeStep<-as.numeric(Sys.time())
if(debug){
fileConn<-file(paste("debug_",timeStep,".txt",sep=""), "a")
}
scaledM<-1/analyteConcentrations/min(1/analyteConcentrations)/weights.scale
#steppedM<-floor(scaledM/weights.step)*weights.step
if(missing(weights.matrix))
{
weights.matrix=
switch(weights.type,
exp=1/exp(scaledM),
factorial=1/factorial(scaledM),
uniform=rep(1,length(scaledM)),
1/(scaledM+1e-5)
);
}
for(j in 1:(moments)) #moments is same as degree!!! so, that is how many diff conc and also how many coef can be estimated
{
Y<-rep(0, lenConc)
curr_moments<-j
for(i in 1:lenConc)
{
Y[i]<-coefficient[[i]][[2]]$coefficient[curr_moments] #here using the fit data, since the coefficient data has been change (divided by R0)
}
conc_t<-poly(1/analyteConcentrations,degree=degree, raw=T)
Y<-(-1)^(curr_moments-1)*Y
lm.f<-lm(formula = Y ~ conc_t
,weights = weights.matrix
#weights=weights.matrix
)
m[j]<-lm.f$coefficient[1]
#plotting
if(debug)
{
jpeg(paste("koff_fitting_degree_",j,"_",timeStep,".jpg",sep=""))
plot(1/analyteConcentrations, Y, main=paste("degree=",j, sep=""),xlab="1/[analyte conc]", ylab="coef fitted")
#pYs<-predictPolynomialFitNoFactorial(1/analyteConcentrations, lm.f$coefficients[-1])+lm.f$coefficients[1]
x<-seq(min(1/analyteConcentrations), max(1/analyteConcentrations), by=(max(1/analyteConcentrations)-min(1/analyteConcentrations))/100)
pxYs<-predictPolynomialFitNoFactorial(x, lm.f$coefficients[-1])+lm.f$coefficients[1]
#lines(1/analyteConcentrations, pYs, col=2, lty=2)
lines(x, pxYs, col=2, lty=2)
dev.off()
#save debugging information for fitting
cat("***********",format(Sys.time(), "%a %b %d %X %Y"),"\n",file=fileConn,sep="")
cat("fitting summaries for fitted coefficient SPR koff rates\n",file=fileConn,sep="")
cat("\t===>degree ",j,":\n",file=fileConn,sep="")
sink(file=fileConn, append=T)
#cat("hello\n")
print(summary(lm.f))
sink();
}
}
if(debug)
{
cat(degree," debugging plots written to ", getwd(),"\n");
cat("fitting information has been written to file at\n\t ");
cat(getwd(),"debug.txt\n")
close(fileConn)
}
nstr<-c(1:(moments-1));
nstr<-paste("moment ", nstr, sep="");
names(m)<-c("Rmax", nstr);
#if(!missing(Rmax))
#{
# if(length(Rmax)<1)
# {
# warning("Rmax has not been set correctly, size of zero")
# Rmax<-1
# }else if(length(Rmax)>degree)
# {
# warning("Rmax has not been set correctly, size larger than expected")
# Rmax<-Rmax[1:degree]
# }
if(missing(Rmax)){
Rmax<-m[1]
}
m[-1]<-m[-1]/Rmax
#}
m
}
####-->this is a wrapper function (or not???) to get all the workering functions together
#### to do fitting by assuming anyway the steady state of reaction.
#'@title fit the sensorgrame to estimate the on-rate distribution information
#'
#'@description estimate the raw moments of the on-rate constant distributin based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param Rmax numeric optional the total Rmax for ligands. When not provided,
#' the output moments is Rmax-scaled. An alternative way practically is to
#' specify the levels of immobilized ligand when the true Rmax is not known.
#' In this case, the output is constant scaled moments. The constant is so-called
#' efficiency of the ligand determined by the orientation as well as activity
#' of the ligand, etc.
#'@param degree numeric degree of the polynomial. Practically the number of data points
#' is big and more than the degree we can estimated. Usually a degree of 100 is probably
#' the best we can do.
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'
#'@return a vector of moments estimated
#'@seealso \code{\link{fitPolySPR}}
#'
#'@note The weights of the fitting is very important in this function. First of all we are doing
#' the polynomial approximation of either exponentional or fraction (Taylor expansion),
#' so we theoratically/ideally need infinite number of degrees for fully accuracy.
#' Practically we, however, could only have some degree of 100 or 200 maximally.
#' This leads to inaccuracy when the independent variable reaches a large value.
#' In order to avoid this inaccuracy, we try to fit the curve with much higher weights
#' for smaller x values than for the larger x values (see weights.type). This way it has
#' better fittings. To control this behavior, two more parameters are specified,
#' weights.step and weights.scale. Weights.step is used to control the local window so
#' as to have equal weights within this window. The bigger the value, the bigger the
#' the window size within which all the data points have equal weights for fitting.
#' This parameter intends to better estimate the noise levels. Parameter weights.scale
#' is used to control how many data points contontribute the fitting, since when value
#' of the independent variable is bigger enough, 1/exp(x) is close to zero, meaning x
#' data point or up contribute nothing to fitting. When weights.scale is small, it
#' increase x values and make less data points contribute to fitting. Otherwise,
#' more data points affect fitting.
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-100
#' dissociationLength<-100
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitSPR.kon(sData[[1]], Rmax=sum(Rmax),debug=TRUE);
#'
#' #first get the distribution of Rmax
#' E_kon<-rep(0,length(kon))
#' E_kon[1]<-sum(Rmax/sum(Rmax)*kon)
#' E_kon[2]<-sum(Rmax/sum(Rmax)*kon^2)
#' E_kon[3]<-sum(Rmax/sum(Rmax)*kon^3)
#' E_kon[4]<-sum(Rmax/sum(Rmax)*kon^4)
#' E_kon[5]<-sum(Rmax/sum(Rmax)*kon^5)
#'
#'
#'@return a list of parameters estimated
#'@export
# mode ==1 regular fitting mode; ==2 flipped and offsetted mode.
fitSPR.kon<-function(data, Rmax, mode=1, degree=100, weights.type="exp", weights.matrix, debug=FALSE, weights.step=1, weights.scale=1)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree<1)
{
stop("degree can not be negative or zero!");
}
if(mode!=1&&mode!=2)
{
stop("unkown fitting mode, please specify again")
}
if(missing(Rmax))
{
if(mode==1)
{
Rmax<-1
}
else{
stop("please specify Rmax")
}
}
SPR.fpc<-NULL
if(missing(weights.matrix))
{
cat("doing fitting the SPR sensorgrams......\n")
if(mode==1){
SPR.fpc<-fitPolySPRs(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale
)
}
else{ #flip &offset mode
SPR.fpc<-fitPolySPRsOffsetted(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale,
Rmax=Rmax
)
}
}
else{
if(class(weights.matrix)!="numeric"||class(weights.matrix)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
cat("doing fitting the SPR sensorgrams......\n")
if(mode==1){
SPR.fpc<-fitPolySPRs(data, degree=degree,
weights.matrix=weights.matrix,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale
)
}
else{ #flip &offset mode
SPR.fpc<-fitPolySPRsOffsetted(data, degree=degree,
weights.type=weights.type,
debug=debug,
weights.step=weights.step,weights.scale=weights.scale,
Rmax=Rmax
)
}
}
#call to fit
cat("doing fitting for moments of kon.......")
m<-fitCoefficientPolySPRs(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, mode=mode)
return(m)
}#end of the function
###---->wrapper function to put together all the code to estimate the off-rate constant distribution
#'@title fit the sensorgrame to estimate the off-rate distribution information
#'
#'@description estimate the raw moments of the off-rate constant distributin based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the individual data set, which might have different quality
#' of the polynomial approximation of the exponential function. The more the quality of approximation
#' the more the number of the moments to be estimated
#' . This is different from kon estimation. Check detail of the implementation here
#' \url{http://}
#'
#'@param data sensorgramData object holding the spr data to be fitted
#'@param Rmax numeric optional the total Rmax for ligands. When not provided,
#' the output moments is Rmax-scaled. An alternative way practically is to
#' specify the levels of immobilized ligand when the true Rmax is not known.
#' In this case, the output is constant scaled moments. The constant is so-called
#' efficiency of the ligand determined by the orientation as well as activity
#' of the ligand, etc.
#'@param degree.fitSPR numeric degree of the polynomial for fitting the SPR responses.
#" Practically the number of data points
#' is big and more than the degree we can estimated. Usually a degree of 100 is probably
#' the best we can do.
#'@param weights.matrix matrix user defined matrix
#'@param debug bool indicate whether to plot the debugging figure
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.step numeric controlling the local window size of equal weights. (see notes)
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (see notes)
#'@param degree.fitMoments numeric the number of degree for the polynomial fitting of
#' moments for off-rate constants.
#'@param moments.num numeric the number of moments to estimate. Practically 6 or 7 is the best we can do.
#'
#'@examples
#' #for non-steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm, sampleFreq=1)
#' fpc<-fitSPR.koff(sData[[1]], Rmax=sum(Rmax),debug=TRUE,degree.fitMoments=length(analyteConcentrations));
#'
#' #expected koff moments
#' E_koff<-rep(0,7)
#' E_koff[1]<-sum(Rmax)
#' E_koff[2]<-sum(koff*Rmax/sum(Rmax))
#' E_koff[3]<-sum(koff^2*Rmax/sum(Rmax))
#' E_koff[4]<-sum(koff^3*Rmax/sum(Rmax))
#' E_koff[5]<-sum(koff^4*Rmax/sum(Rmax))
#' E_koff[6]<-sum(koff^5*Rmax/sum(Rmax))
#' E_koff[7]<-sum(koff^6*Rmax/sum(Rmax))
#'@return A list contains two pieces of results, R0 distributed and Rmax distributed koff.
#' R0 distributed koff contains several distributions based on how many different analyte
#' concentrations are used in the experiment. Each analyte concentration generates
#' one distribution of koff. Rmax distributed koff only has one distribution of koff.
#'@export
fitSPR.koff<-function(data, Rmax, degree.fitSPR=100, weightsType.fitSPR="exp", weightsMatrix.fitSPR, debug=FALSE,
weightsStep.fitSPR=1, weightsScale.fitSPR=1,
degree.fitMoments=1, moments.num=7,
weightsType.fitMoments="exp", weightsMatrix.fitMoments, weightsScale.fitMoments=1
)
{
if(missing(data))
{
stop("data has not been specified!");
}
if(class(data)!="SensorgramData")
{
stop("data is not a valid SensorgramData!");
}
if(degree.fitSPR<1||degree.fitMoments<1||moments.num<1)
{
stop("degree can not be negative or zero!");
}
if(missing(Rmax))
{
Rmax<-1
}
SPR.fpc<-NULL
SPR.lst<-vector("list", length=2)
if(missing(weightsMatrix.fitSPR))
{
cat("doing fitting the SPR sensorgrams......\n")
SPR.fpc<-fitPolySPRsKoff(data, degree=degree.fitSPR,
weights.type=weightsType.fitSPR,
debug=debug,
weights.step=weightsStep.fitSPR, weights.scale=weightsScale.fitSPR
)
}
else{
if(class(weightsMatrix.fitSPR)!="numeric"||class(weightsMatrix.fitSPR)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
cat("doing fitting the SPR sensorgrams......\n")
SPR.fpc<-fitPolySPRsKoff(data, degree=degree.fitSPR,
weights.matrix=weightsMatrix.fitSPR,
debug=debug#,
#weights.step=weights.step,weights.scale=weights.scale
)
}
#call to fit
cat("doing fitting for moments of koff.......")
SPR.fpm=NULL;
if(moments.num>degree.fitSPR){
moments.num<-degree.fitSPR
}
if(missing(weightsMatrix.fitSPR))
{
SPR.fpm<-fitCoefficientPolySPRsKoff(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, degree=degree.fitMoments, moments=moments.num
,weights.type=weightsType.fitMoments, weights.scale=weightsScale.fitMoments
)
}
else{
if(class(weightsMatrix.fitSPR)!="numeric"||class(weightsMatrix.fitSPR)!="matrix")
{
stop("ERROR: the input of weights.matrix is not correctly set!!")
}
SPR.fpm<-fitCoefficientPolySPRsKoff(SPR.fpc, data@analyteConcentrations, Rmax=Rmax, debug=debug, degree=degree.fitMoments, moments=moments.num
,weights.matrix=weightsMatrix.fitMoments
)
}
SPR.lst<-list("R0.distributed"=SPR.fpc, "Rmax.distributed"=SPR.fpm)
return(SPR.lst)
}#end of the function
###--> this is a wrapper function to call the worker function to do the fitting
#'@title fit the sensorgrame steady state data based on the MultiLigand model
#'
#'@description estimate the distributin for equilibrium dissociation constant based on the MultiLigand model
#'
#'@details It assumes the system reached steady state. The moments of distribution of
#' equilibrium dissocation constant are estimated. How many moments of the distribution
#' can be estimated depends on the number of data points. The more the number of
#' of different analyte concentrations are run in the experiment, the
#' more moments we can estimate. Check detail of the implementation here
#' \url{http://}
#'
#'@param x SensorgramData containing the data to be fitted
#'
#'@param steadyStateStart numeric the starting time for steady state
#' optional and if provided, will overwrite the one in sensorgramdata
#'@param steadyStateEnd numeric the ending time for steady state
#' optional and if provided, will overwrite the ones in sensorgramdata
#'@param weights.matrix matrix user defined matrix
#'@param weight.type choose what weight matrix will be used.
# #' \tabular{ll}{
# #' \tab\"exp\", exponential weights with the form 1/exp(x), x is the independent variable\c
# #' }
#' "exp", by default, exponential weights with the form 1/exp(x), x is the independent variable\cr
#' "factorial", factorial weights with the form 1/factorial(x), x is the independent variable\cr
#' "uniform", all data points are weighted equally, with weight being 1 uniformly\cr
#' any other input will be interprated as being weighted by x values of 1/x\cr
#'
#'@param weights.scale numeric controlling the number of data points contributing to the fitting. (\code{\link{fitSPR.kon}})
#'
#'@param fix.ligand NOT IMPLEMENT YET \cr a boolean indicating which ligand immoblization
#' model is using. \cr
#' TRUE, the fixed ligand immobilzation model. Rmax is used\cr
#' FALSE, the variable ligand immbolization model. Rligand
#' and efficiency are used. See also \code{\link{LangmuirModel-class}}
#'@param Rligand NOT IMPLEMENT YET \cr the input for the variable immobilization levels of
#' ligands on the chip. This is used by the variable ligand
#' immbolization model.See also \code{\link{LangmuirModel-class}}
#'@param auto bool used to control the plotting. TRUE by default,
#' so to wait for user in order to continue; FALSE if otherwise for
#' batch running
#'@return a list of parameters estimated
#'@examples
#' #------>for steady state spr
#' kon<-c(2E3, 1E3, 3E3)
#' koff<-c(0.001, 0.01,0.05)
#'
#' analyteConcentrations<-c(
#' 1E-4, 9E-5,6E-5, 3E-5, 2E-5)
#' associationLength<-1000
#' dissociationLength<-1000
#' Rmax<-c(50,80,60)
#'
#' mlgm<- new("MultiLigandModel", kon=kon, koff=koff, analyteConcentrations=analyteConcentrations,
#' associationLength=associationLength, dissociationLength=dissociationLength, Rmax=Rmax
#' )
#' sData<-Simulate(mlgm)
#'
#' #plotting R_AB
#' plot(sData[[1]])
#'
#' #calculate the expected moments of kD distribution
#' e_k<-rep(0,length(analyteConcentrations))
#' e_k[1]<-sum(koff/kon*(Rmax/sum(Rmax)))
#' e_k[2]<-sum((koff/kon)^2*(Rmax/sum(Rmax)))
#' e_k[3]<-sum((koff/kon)^3*(Rmax/sum(Rmax)))
#' e_k[4]<-sum((koff/kon)^4*(Rmax/sum(Rmax)))
#'
#' #fit the steady state data
#' fitSPR.KD(sData[[1]], degree=length(analyteConcentrations), steadyStateStart=950,steadyStateEnd=1000)
#'
#'@export
#input:
####setMethod("FitTwoStateSPR", c("x"="SensorgramData"),#, "sampleFreq"="numeric"),
fitSPR.KD<-function(x, degree=4, steadyStateStart,steadyStateEnd,
weights.matrix, weights.scale=1, weights.type="uniform",
#init.association=NULL, #init.dissociation=NULL,
#control=list(maxiter = 500,tol = 1e-2, minFactor=1/1E10)
#,trace=F,
fix.ligand=TRUE, Rligand, auto=T
)
{
#we first get data ready
if(class(x)!="SensorgramData")
{
stop("the input data is not the correct SensorgramData\n");
# return(FALSE);
}
#get the steady state window ready
ssStart<-x@steadyStateStart
ssEnd<-x@steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(!missing(steadyStateStart))
{
#cat("in hereJ:", steadyStateStart,"\n")
#cat(
ssStart<-steadyStateStart
}
if(!missing(steadyStateEnd))
ssEnd<-steadyStateEnd
#cat("ssStart:",ssStart, ";ssEnd:", ssEnd,"\n");
if(ssStart<0||ssStart<min(x@associationData[,1]))
{
cat("******ERROR:steady state start point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<0||ssEnd>max(x@associationData[,1]))
{
cat("******ERROR:steady state end point is not set correctly...\n");
return(FALSE)
}
if(ssEnd<=ssStart)
{
stop("*******ERROR: steady state end point is not set correctly...\n")
}
#if(mode!=1&&mode!=2)
#{
# cat("******WARNING: unknown value for fitting mode, use \"mode=1\" instead...\n");
# mode<-1;
# #stop("\n");
#}
if(!fix.ligand&&missing(Rligand))
{
stop("ERROR: Rligand can not be null when the variable ligand level is selected");
}
if(!fix.ligand&&(length(Rligand)!=length(x@analyteConcentrations)))
{
stop("ERROR: Rligand doesn't have the correct number of elements compared with analyteConcentrations");
}
if(degree>length(x@analyteConcentrations))
{
warnings("the degree set from the input is larger than the number of analyteConcentrations!!!Automatically set it to the biggest one!",
call.=TRUE, immediate.=TRUE
)
degree<-length(x@analyteConcentrations)
}
output<-NULL
if(missing(weights.matrix))
{
output<-FitSteadyStateSPR(x,degree=degree,
steadyStateStart=steadyStateStart, steadyStateEnd=steadyStateEnd,
weights.scale=weights.scale, weights.type=weights.type,
fix.ligand=fix.ligand, Rligand=Rligand,
auto=auto
)
}
else
{
output<-FitSteadyStateSPR(x,degree=degree,
steadyStateStart=steadyStateStart, steadyStateEnd=steadyStateEnd,
weights.matrix=weights.matrix,
fix.ligand=fix.ligand, Rligand=Rligand,
auto=auto
)
}
return(output)
}#end of the function
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.