FitKineticsGeneLeastSquares {grandR}

Fit a kinetic model according to non-linear least squares.

Description

Fit the standard mass action kinetics model of gene expression using least squares (i.e. assuming gaussian homoscedastic errors) for the given gene. The fit takes both old and new RNA into account and requires proper normalization, but can be performed without assuming steady state. The parameters are fit per Condition.

Usage

FitKineticsGeneLeastSquares(
  data,
  gene,
  slot = DefaultSlot(data),
  time = Design$dur.4sU,
  chase = FALSE,
  CI.size = 0.95,
  steady.state = NULL,
  use.old = TRUE,
  use.new = TRUE,
  maxiter = 250,
  compute.residuals = TRUE
)

Arguments

data

A grandR object

gene

The gene for which to fit the model

slot

The data slot to take expression values from

time

The column in the column annotation table representing the labeling duration

chase

is this a pulse-chase experiment? (see details)

CI.size

A number between 0 and 1 representing the size of the confidence interval

steady.state

either a named list of logical values representing conditions in steady state or not, or a single logical value for all conditions

use.old

a logical vector to exclude old RNA from specific time points

use.new

a logical vector to exclude new RNA from specific time points

maxiter

the maximal number of iterations for the Levenberg-Marquardt algorithm used to minimize the least squares

compute.residuals

set this to TRUE to compute the residual matrix

Details

The start of labeling for all samples should be the same experimental time point. The fit gets more precise with multiple samples from multiple labeling durations. In particular (but not only) without assuming steady state, also a sample without 4sU (representing time 0) is useful.

The standard mass action kinetics model of gene expression arises from the following differential equation:

df/dt = s - d f(t)

This model assumes constant synthesis and degradation rates (but not necessarily that the system is in steady state at time 0). From the solution of this differential equation, it is straight forward to derive the expected abundance of old and new RNA at time t for given parameters s (synthesis rate), d (degradation rate) and f0=f(0) (the abundance at time 0). These equations are implemented in f.old.equi (old RNA assuming steady state gene expression, i.e. f0=s/d), f.old.nonequi (old RNA without assuming steady state gene expression) and f.new (new RNA; whether or not it is steady state does not matter).

This function finds s and d such that the squared error between the observed values of old and new RNA and their corresponding functions is minimized. For that to work, data has to be properly normalized.

For pulse-chase designs, only the drop of the labeled RNA is considered. Note that in this case the notion "new" / "old" RNA is misleading, since labeled RNA corresponds to pre-existing RNA!

Value

A named list containing the model fit:

• data: a data frame containing the observed value used for fitting

• residuals: the computed residuals if compute.residuals=TRUE, otherwise NA

• Synthesis: the synthesis rate (in U/h, where U is the unit of the slot)

• Degradation: the degradation rate (in 1/h)

• Half-life: the RNA half-life (in h, always equal to log(2)/degradation-rate

• conf.lower: a vector containing the lower confidence bounds for Synthesis, Degradation and Half-life

• conf.upper: a vector containing the lower confidence bounds for Synthesis, Degradation and Half-life

• f0: The abundance at time 0 (in U)

• logLik: the log likelihood of the model

• rmse: the total root mean square error

• rmse.new: the total root mean square error for all new RNA values used for fitting

• rmse.old: the total root mean square error for all old RNA values used for fitting

• total: the total sum of all new and old RNA values used for fitting

• type: non-equi or equi

If Condition(data) is not NULL, the return value is a named list (named according to the levels of Condition(data)), each element containing such a structure.

See Also

FitKinetics, FitKineticsGeneLogSpaceLinear, FitKineticsGeneNtr

Examples

sars <- ReadGRAND(system.file("extdata", "sars.tsv.gz", package = "grandR"),
                  design=c("Condition",Design$dur.4sU,Design$Replicate))
sars <- Normalize(sars)
FitKineticsGeneLeastSquares(sars,"SRSF6",steady.state=list(Mock=TRUE,SARS=FALSE))