Construct heatmap for non-steady state (NSS) analysis with improved mechanism score

This uses the output of bakR and a differential expression analysis software to construct a dataframe that can be passed to pheatmap::pheatmap(). This heatmap will display the result of a steady-state quasi-independent analysis of NR-seq data.

Usage

DissectMechanism(
  bakRFit,
  DE_df,
  bakRModel = c("MLE", "Hybrid", "MCMC"),
  DE_cutoff = 0.05,
  bakR_cutoff = 0.3,
  Exp_ID = 2,
  sims = 1e+07
)

Arguments

bakRFit: bakRFit object
DE_df: dataframe of required format with differential expression analysis results. See Further-Analyses vignette for details on what this dataframe should look like
bakRModel: Model fit from which bakR implementation should be used? Options are MLE, Hybrid, or MCMC
DE_cutoff: padj cutoff for calling a gene differentially expressed
bakR_cutoff: padj cutoff for calling a fraction new significantly changed. As discussed in the mechanistic dissection vignette, it is best to keep this more conservative (higher padj) than is typical. Thus, default is 0.3 rather than the more standard (though admittedly arbitrary) 0.05.
Exp_ID: Exp_ID of experimental sample whose comparison to the reference sample you want to use. Only one reference vs. experimental sample comparison can be used at a time
sims: Number of simulation draws from null distribution for mechanism p value calculation

Value

returns list of data frames: heatmap_df and NSS_stats. The heatmap_dfdata frame can be passed to pheatmap::pheatmap(). The NSS_stats data frame contains all of the information passed to NSS_stats as well as the raw mechanism scores. It also has p values calculated from the mechanism z scores.

Details

Unlike NSSHeat, DissectMechanism uses a mechanism scoring function that is not discontinuous as the "degradation driven" vs. "synthesis driven" boundary. Instead, the score approaches 0 as the function approaches the boundary from either side.

In addition, DissectMechanism now defines a null model for the purpose of p value calculation using the mechanism score. Under the null hypothesis, the mechanism score is the product of two normal distributions with unit variance, one which has a non-zero mean. Simulation is used to estimate the integral of this distribution, and the number of draws (which determines the precision of the p value estimate) is determined by the sims parameter.

DissectMechanism also provides "meta-analysis p values", which can be interpreted as the p-value that a particular RNA feature is observing differential expression or differential kinetics (or both). This meta_pval is estimated using Fisher's method for meta analysis.

Examples

# \donttest{
# Simulate small dataset
sim <- Simulate_bakRData(100, nreps = 2)

# Analyze data with bakRFit
Fit <- bakRFit(sim$bakRData)
#> Finding reliable Features
#> Filtering out unwanted or unreliable features
#> Processing data...
#> Estimating pnew with likelihood maximization
#> Estimating unlabeled mutation rate with -s4U data
#> Estimated pnews and polds for each sample are:
#> # A tibble: 4 × 4
#> # Groups:   mut [2]
#>     mut  reps   pnew    pold
#>   <int> <dbl>  <dbl>   <dbl>
#> 1     1     1 0.0502 0.00101
#> 2     1     2 0.0501 0.00101
#> 3     2     1 0.0502 0.00101
#> 4     2     2 0.0501 0.00101
#> Estimating fraction labeled
#> Estimating per replicate uncertainties
#> Estimating read count-variance relationship
#> Averaging replicate data and regularizing estimates
#> Assessing statistical significance
#> All done! Run QC_checks() on your bakRFit object to assess the 
#>             quality of your data and get recommendations for next steps.

# Number of features that made it past filtering
NF <- nrow(Fit$Fast_Fit$Effects_df)

# Simulate mock differential expression data frame
DE_df <- data.frame(XF = as.character(1:NF),
                       L2FC_RNA = stats::rnorm(NF, 0, 2))

DE_df$DE_score <- DE_df$L2FC_RNA/0.5
DE_df$DE_se <- 0.5

DE_df$DE_pval <- 2*stats::dnorm(-abs(DE_df$DE_score))
DE_df$DE_padj <- 2*stats::p.adjust(DE_df$DE_pval, method = "BH")

# perform NSS analysis
NSS_analysis <- DissectMechanism(bakRFit = Fit,
               DE_df = DE_df,
               bakRModel = "MLE")
#> Combining bakR and DE analyses
#> Calculating mechanism p-value. This could take several minutes...
#> Assessing differential synthesis
#> Constructing Heatmap_df

# }