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

Source:`R/modeling.R`

`FitKineticsGeneNtr.Rd`

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"
#>
#>
```