---
title: "Standard kinetic parameter estimation: EstimateKinetics()"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{EstimateKinetics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.dim = c(6, 4)
)
```

## Introduction

The goal of most standard, single label, NR-seq experiments is to assess the
kinetics of RNA synthesis and degradation. In EZbakR, this is the task of
`EstimateKinetics()`. This vignette will discuss the basics of the standard
analysis implemented by `EstimateKinetics()`, and touch on two alternative
analysis strategies that `EstimateKinetics()` also implements: non-steady-state
and short-feed analyses.


```{r setup}
library(EZbakR)
```

## The basics

### Running the standard kdeg and ksyn estimation

Below I demonstrate how to estimate synthesis and degradation rate constants
using the standard analysis workflow and simulated data:

```{r}
# Simulate data
simdata <- EZSimulate(nfeatures = 300, nreps = 2)

# Create an `EZbakRData` object
ezbdo <- EZbakRData(simdata$cB, simdata$metadf)

# Estimate fractions
ezbdo <- EstimateFractions(ezbdo)

# Estimate synthesis and degradation rate constants
ezbdo <- EstimateKinetics(ezbdo,
                          strategy = "standard")
```

Setting `strategy = "standard"` is not strictly necessary, but I do it here
to highlight that this is an argument that you can alter. If you have multiple
`fractions` tables, then you will need to specify aspects of the metadata of
the table you want to use. Under the hood, EZbakR is passing these arguments to
`EZget()`, so check out the vignette on `EZget()` to better understand how this
works:

```{r}
# Create another fractions table
ezbdo2 <- EstimateFractions(ezbdo, overwrite = FALSE)

# Estimate synthesis and degradation rate constants
  # Use the 1st identical fractions table you created
ezbdo2 <- EstimateKinetics(ezbdo2,
                          repeatID = 1,
                          strategy = "standard")
```

In this case, `repeatID` has to be specified because the two fractions table
present are identical in all other metadata aspects. Other metadata that you
can specify are:

* `features`: Set of features analyzed, or a subset of features unique to the
table of interest.
* `populations`: Mutation populations analyzed. Relevant if one of the tables
is from a single-label analysis, and the other is from a multi-label analysis.
`EstimateKinetics()`'s standard analysis strategy is only compatible with
single-label analyses.
* `fraction_design`: The fraction design table used for the analysis of interest.
Rarely will be necessary to specify in this case.

### What is the "standard" strategy?

The standard strategy makes the following assumptions:

1. You have conducted a pulse-label experiment. That means, cells were labeled
for a period of time that I will refer to as the label time, and RNA was
extracted after this labeling. This is in contrast to a pulse-chase design, where
you label for a certain amount of time (the pulse time) and then chase with 
regular nucleotide for a certain amount of time (the chase time). Eventually,
EZbakR will support kinetic parameter inference from a pulse-chase design, but
I will say more about that at the end of this vignette.
1. All of the populations of RNA you analyzed were at steady-state during the 
course of metabolic labeling. This means that their levels were constant and 
not changing in response to some perturbation during the labeling.
1. RNA synthesis is a 0-th order, constant rate process, and RNA degradation is
a 1st-order, single rate constant process.

If all of these assumptions are met, then the dynamics of a given RNA population
are well described by the following differential equation:

$$
\begin{align}
\frac{\text{dR}}{\text{dt}} = k_{\text{syn}} - k_{\text{deg}}*\text{R} 
\end{align}
$$
where $\text{R}$ represents the concentration of RNA over time and 
$k_{\text{syn}}$ and $k_{\text{deg}}$ are the synthesis and degradation rate
constants, respectively. The general solution to this differential equation is:

$$
\begin{align}
\text{R(t)} = \frac{ k_{\text{syn}}}{k_{\text{deg}}} +  (\text{R}_\text{o} - \frac{ k_{\text{syn}}}{k_{\text{deg}}})* e^{-\text{tl}*k_{\text{deg}}}
\end{align}
$$

where $\text{R}_\text{o}$ is the initial concentration of the RNA and $\text{tl}$
is the label time. Thus the concentration of new RNA ($\text{R}_\text{o} = 0$) 
as a function of time is:

$$
\begin{align}
\text{R}_{\text{new}}\text{(t)} = \frac{ k_{\text{syn}}}{k_{\text{deg}}}*(1 - e^{-\text{tl}*k_{\text{deg}}})
\end{align}
$$
As we are assuming that the total RNA concentration is not changing during the
labeling, and since the steady-state concentration of RNA (solve for R setting 
$\frac{\text{dR}}{\text{dt}} = 0$) is $\frac{ k_{\text{syn}}}{k_{\text{deg}}}$,
we get that the fraction of RNA that is new (denoted $\theta$) is:

$$
\begin{align}
\theta = 1 - e^{-\text{tl}*k_{\text{deg}}}
\end{align}
$$
We can thus solve for the degradation rate constant, and use the normalized 
read counts to infer the synthesis rate constant. 

$$
\begin{align}
k_{\text{deg}} &= -\frac{\text{log}(1 - \theta)}{\text{tl}} \\
k_{\text{syn}} &= (\text{normalized read count})*k_{\text{deg}}
\end{align}
$$

**NOTE**: degradation rate constants estimated in this manner will be in 
absolute units (1/time), whereas synthesis rate constants will be in relative 
units (read counts/time). Thus, degradation rate constant estimation is 
"internally normalized" (the fraction new is a ratio of two unnormalized 
quantities that would have the same normalization scale factor), whereas 
synthesis rate constant estimates require normalization across libraries and 
samples. At this stage, EZbakR uses TMM normalization (as in packages like edgeR 
and DESeq2) to calculate the normalized read counts. In the near future, it will 
also support users providing external scale factors (e.g., those derived from 
spike-ins).

## Short feed analyses

Sometimes it is advantageous to use a fairly short label time. For example, 
NR-seq is commonly used as a replacement for enrichment based metabolic labeling
methods (e.g., TT-seq), and a short label time is used to ensure that there is
limited degradation of newly synthesized RNA during the feed. EZbakR formalizes
this strategy and provides a way by which to estimate synthesis and degradation
rate constants in this context.


### How to run it

Everything is the same as it was in the standard case, except now we set
`strategy = "shortfeed"`:

```{r}
# Estimate synthesis and degradation rate constants
ezbdo <- EstimateKinetics(ezbdo,
                          strategy = "shortfeed")
```


### What is the "shortfeed" strategy?

In the case of an extremely short label feed time (relative to the average
half-life of RNA being analyzed), we can assume that there is no degradation
newly synthesized RNA during the label time. Thus, the following differential
equation and solution to said ODE describes the dynamics of the new RNA:

$$
\begin{align}
\frac{\text{dR}}{\text{dt}} &= k_{\text{syn}} \\
\text{R}_{\text{new}}\text{(t)} &=  k_{\text{syn}}*\text{t}
\end{align}
$$
Assuming that the concentration of new RNA is proportional to the normalized 
number of new reads yields the following estimate for $k_{\text{syn}}$:

$$
k_{\text{syn}} =  \frac{\theta*(\text{normalized read count})}{\text{tl}}
$$
Assuming the RNA levels are at steady-state (as in the standard analysis),
yields the following estimate for the degradation rate constant:

$$
\begin{align}
\text{normalized read count} &= \frac{k_{\text{syn}}}{k_{\text{deg}}} \\
k_{\text{deg}} &= \frac{k_{\text{syn}}}{\text{normalized read count}} \\
k_{\text{deg}} &= \frac{\theta}{\text{tl}}
\end{align}
$$

The calculus enthusiasts among you may recognize that this estimate for
$k_{\text{deg}}$ is equivalent to a 1st order Taylor series approximation of
the "standard" estimate strategy. This is no coincidence, as the assumption of
a short feed is equivalent to assuming that $\theta$ is close to 0 (most of the
RNA is old). It may also seem odd that this estimate is bounded between 0 
and 1/tl, but this just reflects the assumption made by this strategy that RNA
half lives are much shorter than the label time.


## Non-steady-state (NSS) analyses

**NOTE:** The strategy described here is a slight modification of that discussed
in [Narain et al., 2021](https://www.sciencedirect.com/science/article/pii/S1097276521004962). If you use EZbakR's "NSS" strategy, please cite that
paper as well as EZbakR.

One of the key assumptions made by the "standard" analysis strategy is the 
assumption that RNA levels are at steady-state during the labeling. This 
assumption can be violated if a perturbation was applied to cells shortly
before the start of, or during, labeling. It can also be violated if you are
studying a process (e.g., development) where the state of the cells is changing
over the course of the experiment. `EstimateKinetic()` can adjust its kinetic
parameter estimation strategy accordingly, though it requires making some
assumptions about how the experiment was performed. See below for details

### Requirements of the "NSS" strategy

In this case, it is important to discuss some key aspects about how this
strategy works before showcasing it. Namely:

1. Instead of assuming the fraction new in a given sample is related to the
degradation rate constant, EZbakR will compare the old read levels in the 
labeled sample and a set of samples that provide information about RNA abundances
at the start of the labeling.
1. EZbakR assumes that your metadf contains two additional columns. One should
specify which samples belong to the same biological condition/analysis time point
and thus should assay similar RNA populations. The other should specify which
group of samples identified by the first column provide information about the
RNA population that existed at the start of labeling.

For example, suppose you have a time course with the following samples:

1. "sampleA": -s4U data collected at the start of the time course
1. "sampleB": +s4U data collected after a 2 hour labeling period that started
at the start of the time course.
1. "sampleC": +s4U data collected after a 2 hour labeling period that started
at a time point equivalent to the end of the labeling in "sampleB".

In this case, "sampleA" provides information about the initial RNA concentrations
for "sampleB", and "sampleB" provides information about the initial RNA concentrations
for "sampleC". Therefore, your metadf could look like:

```{r}
metadf <- data.frame(
  sample = c('sampleA',
             'sampleB',
             'sampleC'),
  tl = c(0, 2, 2),
  group = c("a", 
            "b", 
            "c"),
  reference = c("a",
                "a",
                "c")
)
```

One thing to note is that the value you enter for "sampleA"'s reference is
technically arbitrary, as it is a -s4U sample that you will not be estimating
degradation rate constants from. Thus, you could enter any string for its value
in the column named "reference" in this example.

### How to run it

The only difference in this case is that you have to tell EZbakR which metadf
columns correspond to the two additional columns explained above. The relevant
parameters are:

1. grouping_factor: Which column defines groups of samples which assay the same
RNA population?
1. reference_factor: Which column defines the groups of samples which assay the
RNA population that existed at the start of labeling?

```{r}
# new metadf
metadf <- data.frame(sample = paste0('sample', 1:6),
                 tl = c(2, 2, 2, 2, 0, 0),
                 group = c('b', 'b',
                           'c', 'c',
                           'a', 'b'),
                 ref = c('a', 'a',
                         'b', 'b',
                         'a', 'a'))

# Create an `EZbakRData` object
ezbdo <- EZbakRData(simdata$cB, metadf)

# Estimate fractions
ezbdo <- EstimateFractions(ezbdo)

# Estimate synthesis and degradation rate constants
ezbdo <- EstimateKinetics(ezbdo,
                          strategy = "NSS",
                          grouping_factor = "group",
                          reference_factor = "ref")
```

In this case, the only experimental details column of the metadf is called
"treatment", so this is specified in `grouping_factors`.


### What is the "NSS" strategy?

You can see [Narain et al.](https://www.sciencedirect.com/science/article/pii/S1097276521004962) for an alternative discussion of this strategy.

We would like to relax the assumption that the total RNA levels are unchanging
over the course of the labeling. We can do this by considering the dynamics
of the old RNA during the labeling. If we assume that a single, unchanging
degradation rate constant describes their turnover, then the following ODE
and solution describe its dynamics:

$$
\begin{align}
\frac{\text{dR}}{\text{dt}} &= -k_{\text{deg}}*\text{R} \\
\text{R}_{\text{old}}\text{(t)} &=  \text{R}_{\text{o}}*e^{-k_{\text{deg}}*\text{t}}
\end{align}
$$
In the steady-state case, $\text{R}_{\text{o}}$ would be equivalent to the 
steady-state RNA level ($\frac{ k_{\text{syn}}}{k_{\text{deg}}}$). In the away
from steady-state case though, this assumption is no longer valid. Despite this,
assume for a second that we could estimate $\text{R}_{\text{o}}$. Under this
assumption, we could then estimate $k_{\text{deg}}$ as follows:

$$
k_{\text{deg}} = -\frac{\text{log}(\text{R}_{\text{old}}\text{(t)}/\text{R}_{\text{o}})}{\text{tl}}
$$

We can then return to the general solution of the differential equation described
in the 'standard' analysis strategy section:

$$
\begin{align}
\text{R(t)} &= \frac{ k_{\text{syn}}}{k_{\text{deg}}} +  (\text{R}_\text{o} - \frac{ k_{\text{syn}}}{k_{\text{deg}}})* e^{-\text{tl}*k_{\text{deg}}} \\
\text{R}_{\text{new}}\text{(t)} &= \frac{ k_{\text{syn}}}{k_{\text{deg}}}*(1 - e^{-\text{tl}*k_{\text{deg}}})
\end{align}
$$
where in the second line we set $\text{R}_\text{o} = 0$ as there was no new
RNA at the start of labeling, and derive an estimator for $k_{\text{syn}}$:

$$
k_{\text{syn}} = \frac{\theta*(\text{normalized read count)}*k_{\text{deg}}}{1 - e^{-k_{\text{deg}}*\text{tl}}}
$$
This yields the same estimator as in the steady-state case when 
$\theta = 1 - e^{-k_{\text{deg}}*\text{tl}}$. To make this strategy work, we 
just need to estimate the starting levels of old RNA. Thus, you need some
sort of RNA-seq data from a timepoint equivalent to the start of labeling.

One thing you may be wondering is just how much of non-steady-state analysis 
this really is if we are still assuming constant rate constants? What if the
rates of synthesis and degradation are being actively modulated by the cell
during the label time? I won't present a detailed proof here, but it turns out
that the strategy laid out above is the best that you can do. In the case where
rate constants are changing during the label time, you can think of the estimates
from this strategy as time averages of the true rate constants, averaged over
the labeling period.

## Pulse-chase analyses

The infrastructure is in place for EZbakR to support pulse-chase analyses of
kinetic parmaeters (e.g., metadf can accept specification of a pulse and chase
time). Thus, in the future, EZbakR will provide such support. That being said,
I would like to take some time to argue against *almost* ever doing a pulse-chase
NR-seq experiment. Pulse-chase experimental designs suffer from a number of 
shortcomings and are almost never necessary in NR-seq experiments. 

There are many classic examples of pulse-chase experimental designs being used 
to assess the kinetics of RNA metabolism. In these cases though, the point of 
the pulse was to create a species of RNA whose dynamics are completely driven 
by the degradation kinetics of the RNA. For example, the pioneering studies 
from the Parker lab elucidating the mechanisms of RNA degradation involved 
pulsing with a particular nutrient to stimulate transcription of a construct. 
The chase then washed out that nutrient to shut off transcription from the 
construct, meaning that RNA produced from it will exclusively degrade post-chase. 
When you pulse with a metabolic label like s^4^U though, you have already 
created a species of RNA whose dynamics are completely degradation driven: the 
unlabeled RNA. Following the metabolic label pulse with a nucleotide chase is 
thus redundant, switching the degradation driven population from the unlabeled 
RNA to the labeled RNA. 

In addition, pulse-chase experiments often necessitate extensive exposure of 
your cells to the metabolic label, which increases the 
possibility of cytotoxic effects of metabolic label. Finally, analysis of 
pulse-chase experiments is more complicated than a pulse-label design. To 
estimate the kinetics of degradation, you need to know what fraction of the RNA 
for each feature of interest was labeled after the pulse. This means that 
estimating the synthesis and degradation kinetics requires two separate samples, 
RNA extracted after the pulse and RNA extracted after the chase. With the pulse 
alone, you can get all of the information that you can get with both though!! 
You want multiple time points? Just do multiple pulses. Pulse-labels are easier 
to perform and easier to analyze than pulse-chases.
