Fit a kinetic model using the degradation rate transformed NTR posterior distribution.

Fit the standard mass action kinetics model of gene expression by maximum a posteriori on a model based on the NTR posterior. The fit takes only the NTRs into account and is completely independent on normalization, but it cannot be performed without assuming steady state. The parameters are fit per Condition.

Usage

FitKineticsGeneNtr(
  data,
  gene,
  slot = DefaultSlot(data),
  time = Design$dur.4sU,
  CI.size = 0.95,
  transformed.NTR.MAP = TRUE,
  exact.ci = FALSE,
  total.fun = median
)

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
CI.size: A number between 0 and 1 representing the size of the credible interval
transformed.NTR.MAP: Use the transformed NTR MAP estimator instead of the MAP of the transformed posterior
exact.ci: compute exact credible intervals (see details)
total.fun: use this function to summarize the expression values (only relevant for computing the synthesis rate s)

Value

A named list containing the model fit:

data: a data frame containing the observed value used for fitting
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
total: the total sum of all new and old RNA values used for fitting
type: always "ntr"

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.

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.

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. Further assuming steady state allows to derive the function transforming from the NTR to the degradation rate d as $d(ntr)=-1/t log(1-ntr)$. Furthermore, if the ntr is (approximately) beta distributed, it is possible to derive the distribution of the transformed random variable for the degradation rate (see Juerges et al., Bioinformatics 2018).

This function primarily finds d by maximizing the degradation rate posterior distribution. For that, data does not have to be normalized, but this only works under steady-state conditions. The synthesis rate is then computed (under the assumption of steady state) as $s=f0 \cdot d$

The maximum-a-posteriori estimator is biased. Bias can be removed by a correction factor (which is done by default).

By default the chi-squared approximation of the log-posterior function is used to compute credible intervals. If exact.ci is used, the posterior is integrated numerically.

Examples

sars <- ReadGRAND(system.file("extdata", "sars.tsv.gz", package = "grandR"),
                  design=c("Condition",Design$dur.4sU,Design$Replicate))
#> Warning: Duplicate gene symbols (n=1, e.g. MATR3) present, making unique!
sars <- Normalize(sars)
sars <- subset(sars,columns=Condition=="Mock")
FitKineticsGeneNtr(sars,"SRSF6")
#> $Mock
#> $Mock$data
#>        Name     alpha      beta t     Quantile Expected.NTR Observed.NTR
#> 1 Mock.1h.A  371.4187 1277.9872 1 0.9739303846    0.2454541    0.2248497
#> 2 Mock.2h.A  604.6440  833.4801 2 0.7841455206    0.4306605    0.4203286
#> 3 Mock.2h.B  759.8269  987.6798 2 0.3640291596    0.4306605    0.4347316
#> 4 Mock.3h.A 1297.4445 1053.6619 3 0.9651349807    0.5704073    0.5518884
#> 5 Mock.4h.A 1154.1171  453.4019 4 0.0001271576    0.6758526    0.7182208
#>   Absolute.Residual Relative.Residual
#> 1       -0.02060446      -0.083944217
#> 2       -0.01033197      -0.023990985
#> 3        0.00407104       0.009453014
#> 4       -0.01851891      -0.032466111
#> 5        0.04236819       0.062688502
#> 
#> $Mock$Synthesis
#> [1] 344.0347
#> 
#> $Mock$Degradation
#> [1] 0.2816392
#> 
#> $Mock$`Half-life`
#> [1] 2.461117
#> 
#> $Mock$conf.lower
#>   Synthesis Degradation   Half-life 
#> 330.8044664   0.2708084   2.3677358 
#> 
#> $Mock$conf.upper
#>   Synthesis Degradation   Half-life 
#> 357.6031805   0.2927468   2.5595479 
#> 
#> $Mock$f0
#> [1] 1221.544
#> 
#> $Mock$logLik
#> [1] -5631.111
#> 
#> $Mock$rmse
#> [1] 0.01918177
#> 
#> $Mock$total
#> [1] 1221.544
#> 
#> $Mock$type
#> [1] "ntr"
#> 
#>