Steady-state quasi-independent mechanistic investigations

Violating the steady-state assumption

Before we talk about what to do when the steady-state assumption is questionable, it’s worth briefly mentioning when the steady-state assumption is worth throwing out. Assuming steady-state RNA dynamics means that you are assuming the rate at which RNA synthesized is the same as the rate at which it is degraded. “Wait, does that mean we are assuming \(ksyn = kdeg\)?” No. The degradation rate constant \(kdeg\) is not the rate of degradation. The rate of degradation (i.e., the number of RNA molecules degraded per unit time) depends on \(kdeg\) AND the amount of RNA present. So steady-state actually means that \(ksyn = kdeg*[RNA]\). In the context of NR-seq experiments, assuming steady-state means that we are assuming this relationship holds throughout the entire metabolic label feed time. Therefore, another way to phrase this assumption is to say that we are assuming that our population of cells aren’t actively responding to some stimulus and regulating \(ksyn\) and/or \(kdeg\). RNA levels remain constant during the experiment.

Here’s a fun question: are cells ever at steady-state? Yes and no. An individual cell is almost never at steady-state. Gene expression is constantly regulated as the cell progresses through the cell cycle and performs all of the biochemical functions necessary to stay alive. BUT, bakR is designed to analyze bulk nucleotide recoding RNA-seq data. This means the data comes from a population of cells, not a single cell. Thus, when I say that I am assuming steady-state, what I really mean is that I am assuming the population average is at steady-state. If the population of cells are asynchronous (so they aren’t all going through the exact same cell cycle stages at the same time) and haven’t been perturbed recently (e.g., treated with a drug), then this is probably a pretty solid assumption.

Lots of cool experiments are done with cells far from steady-state though. When using nucleotide recoding RNA-seq, it is increasingly popular to perform rapid degradation of a protein of interest, treatment with a specific drug, or some other acute perturbation leading up to and during metabolic labeling. Therefore, we need a way to analyze data without the helpful steady-state assumption. This is challenging though, because if we can’t assume steady-state, how do we relate what we measure (the fraction new) to the RNA metabolic kinetics we care about? The answer: keep it non-parametric.

Analyzing NR-seq data far from steady-state

I am using the term non-parametric here to suggest getting away from using explicit mathematical functions to model RNA metabolism. Rather, we can think very generally about what effects changes in RNA synthesis and degradation kinetics can have on the fraction new in an NR-seq experiment. Let’s say we compare an unperturbed, steady-state population of cells to a cells that have been treated with a drug of interest. If the fraction new for a particular transcript is higher in the drug treated cells than the untreated cells, how could that have come about? Well, generally speaking, there are two possibilities:

1. The RNA may have been destabilized by the drug treatment. This is the same explanation you would give for the change if everything were steady-state.
2. The synthesis rate of the RNA may have been ramped up following drug treatment. This is an intuitive explanation that is surprisingly only possible away from steady-state. At steady-state, changes in synthesis kinetics affect the amount of new and old RNA equally, leading to no change in the fraction new.

Here is the key realization: these two possibilities have different impacts on the LEVELS of RNA. If synthesis is increasing upon drug treatment, then the amount of that transcript will also increase. If the transcript is destabilized though, then there will be less of it around at the end of the experiment. Therefore, combining differential expression with differential fraction new analyses reveals the kinetic differences. More generally, here is a diagram that covers the range of possible kinetic mechanisms and conclusions (\(\theta\) is the fraction new):

Diagram of NR-seq NSS analysis strategy.

The statistics of DissectMechanism

DissectMechanism attempts to formalize the analysis discussed in the last section. It combines a differential fraction new test statistic with a differential expression test statistic. The differential fraction new test statistic is just the L2FC(kdeg) z-score output by bakR. Isn’t this invalidated by violations of the steady-state assumption? Yes and no. While it is true that bakR’s estimated L2FC(kdeg) may no longer truly represent the change in the degradation rate constant, it still does represent an accurate assessment of the change in the fraction new. The kdeg is estimated as:

\[ k_{\text{deg}} = -\text{ln}(1 - fn)/\text{tl} \] where tl is the label time. While the transformation on the right of this expression may no longer represent the degradation rate constant away from steady-state, it is simply a monotonic transformation of the fraction new. Thus, differences in bakR’s estimated “kdeg” can still be interpreted as differences in the fraction new, assuming the label times in the conditions being compared are the same. If the label times differ, this also represents a naive (though somewhat inaccurate in the non-steady-state case) normalization for these differences.

The mechanism test statistic used by DissectMechanism is quite simple. It’s just:

\[ \text{mech_stat} = [\text{L2FC(kdeg) z-score}]*[\text{L2FC(RNA) z-score}] \]

Why this is a reasonable test-statistic choice is best demonstrated via visualization of this function:

## Grid to plot on

# Grid size in each axis
n <- 300

# DE z-scores
DE_score <- rep(seq(from = -10, to = 10, length.out = n), times = n)

# bakR z-scores
bakR_score <- rep(seq(from = -10, to = 10, length.out = n), each = n)

grid_df <- tibble(DE = DE_score,
                  bakR = bakR_score)


# Mechanism score
grid_df <- grid_df %>%
  mutate(Mech = ifelse(DE > 0, 
                       (bakR + 2)*DE,
                       (bakR - 2)*DE))

# Calculate p-value
grid_df <- grid_df %>%
  mutate(Mech_pval = stats::df(abs(Mech), df1 = 2, df2 = 2, ncp = 2))

# Visualize test statistic
ggplot(grid_df, aes(x = bakR, y = DE, z = 0.03*Mech)) + 
  geom_contour_filled() + 
  theme_classic() + 
  xlab("bakR score") + 
  ylab("DE score") 

In this plot, dark blue represents high confidence post-transcriptional regulation and bright yellow represents high confidence transcriptional regulation. There are four somewhat asymetric quadrants in this plot:

Top-left: Fraction new went down and expression levels went up. This has to be stabilization.

Top-right: Fraction new went up and expression levels went up. This has to be transcriptional upregulation.

Bottom-left: Fraction new went down and expression went down. Transcriptional downregulation must be at play

Bottom-right: Fraction new went up and expression went down. The RNA was destabilized.

This particular function is admittedly a bit arbitrary, but it has a lot of the trends that you would want in a function attempting to formalize to the strategy discussed in the last section:

1. The quadrants described above match the mechanistic conclusions in the schematic shown in the last section.
2. Low confidence in differential fraction new or differential expression can be balanced out by high confidence in the other. Moderate confidence in both is as good as high confidence in one and lower confidence in another.
3. In keeping with the possibility of steady-state dynamics, differential expression can be transcriptionally driven without a strong fraction new signal. More on this in the next section when I discuss the vertical dividing lines that determine just how strong the fraction new signal has to be to qualify as “real” signal.

The assymetry of the mechanism score

I briefly alluded in the last section to the asymmetry of the mechanism score function. The x-axis (0 differential expression) represents one dividing line, but what about the two separate vertical dividing lines? The reason for this is that transcriptional regulation can yield no fraction new signal. Thus, low confidence changes in the fraction new should be interpreted as no change in the fraction new. When no change in the fraction new is accompanied by differential expression, this is most likely transcriptional regulation.

The locations of the vertical lines are set by the bakR_cutoff parameter in DissectMechanism. bakR_cutoff represents the bakR_padj cutoff for calling a change in the fraction new “significant”. DissectMechanism determines what z-score is necessary to yield that padj (since multiple-test adjustment makes the relationship between padj and z-score a bit difficult to predict a priori) and uses that as the vertical dividing lines.

The reliance on statistical confidence is an important shortcoming of DissectMechanism’s analysis strategy. A low confidence L2FC(kdeg) could be the result of no large changes in the fraction new, but it could also be the result of low statistical power. In other words, there may be strong post-transcriptional regulation lurking in the region between these vertical dividing lines, but the coverage or number of replicates may have been too low to assess this confidently.

The consequence of all of this is that bakR_cutoff tunes the kind of false positives that are more likely to crop up:

1. If bakR_cutoff is small, then only large fraction new z-scores are going to be called significant. Any other low confidence fraction new changes that are accompanied by differential expression will be called transcriptional regulation. Thus, small bakR_cutoff’s yield more transcriptional false positives.
2. On the other end of the spectrum, a bakR_cutoff close to 1 will cause any small change in the fraction new to be interpreted as a real change. Thus, if the direction of fraction new differences is opposite the direction of differential expression, DissectMechanism will call this post-transcriptional regulation, no matter how inconsequential the fraction new shift. Thus, large bakR_cutoff’s yield more post-transcriptional false positives.

The conclusion: be wary. As discussed in the summary at the top of this vignette, bakR_cutoff should be chosen in a way to limit the false positives that would be more likely to lead you astray.

The mechanism p-value

The last wrinkle in DissectMechanism is the calculation of the mechanism score’s p-value. Multiple test adjusted p-values are nice because they give you a clear sense of how confident you can be that a given set of mechanistic assignments are accurate. Without these, it can be difficult to interpret the mechanism scores. It turns out though that nature of the mechanism score null model makes it difficult to calculate these p-values with great precision.

Work through a bit of mathematical statistics and you’ll find that the null hypothesis distribution is the esoteric “doubly non-central scaled F-distribution”. There is the aptly named sadist R package which implements functions for working with little known distributions, but even it only accomodates the unscaled version of this distribution. While one R package exists to calculate integrals of the doubly non-central scaled F-distribution, the precision of its estimates is about as good as I was able to shakily code up myself. Therefore, I went forward with a home brew Monte Carlo approach to p-value calculation. In short, I simulate from the null model a bunch of times (the number of which is set by the sims parameter in DissectMechanism), and calculate the p-value as the fraction of the time the simulated mechanism score was larger in magnitude than a given mechanism score.

This means that the precision of the p-value is set by the sims parameter. By default, it is set to 10 million, which means that the minimum meaningful p-value it can estimate is 1/10,000,000, or 1e-7. In this vignette, I decreased this by an order of base-10 magnitude to speed up the runtime, as it can take 10 or more seconds to collect this many draws from the null. You can get more precision by increasing sims, but it’s going to cost you heavily on the runtime.

NOTE: if 0 of the draws from the null are more extreme than a particular mechanism score, that mechanism score is assigned a p-value of 0.5/sims. Because of this, relative ordering of mechanism scores is not perfectly preserved in the relative ordering of p-values. Thus, if you are attempting to generate ordered lists of high confidence mechanistic assignments, order by mech_score, not mech_pval!!