Skip to contents

Compute the posterior log2 fold change distributions of RNA synthesis and degradation

Usage

EstimateRegulation(
  data,
  name.prefix = "Regulation",
  contrasts,
  reference.columns = NULL,
  slot = DefaultSlot(data),
  time.labeling = Design$dur.4sU,
  time.experiment = NULL,
  ROPE.max.log2FC = 0.25,
  sample.f0.in.ss = TRUE,
  N = 10000,
  N.max = N * 10,
  CI.size = 0.95,
  seed = 1337,
  dispersion = NULL,
  sample.level = 2,
  correct.labeling = FALSE,
  verbose = FALSE
)

Arguments

data

the grandR object

name.prefix

the prefix for the new analysis name; a dot and the column names of the contrast matrix are appended; can be NULL (then only the contrast matrix names are used)

contrasts

contrast matrix that defines all pairwise comparisons, generated using GetContrasts

reference.columns

a reference matrix usually generated by FindReferences to define reference samples for each sample; can be NULL if all conditions are at steady state (see details)

slot

the data slot to take f0 and totals from

time.labeling

the column in the Coldata table denoting the labeling duration, or the numeric labeling duration itself

time.experiment

the column in the Coldata table denoting the experimental time point (can be NULL, see details)

ROPE.max.log2FC

the region of practical equivalence is [-ROPE.max.log2FC,ROPE.max.log2FC] in log2 fold change space

sample.f0.in.ss

whether or not to sample f0 under steady state conditions

N

the sample size

N.max

the maximal number of samples (necessary if old RNA > f0); if more are necessary, a warning is generated

CI.size

A number between 0 and 1 representing the size of the credible interval

seed

Seed for the random number generator

dispersion

overdispersion parameter for each gene; if NULL this is estimated from data

sample.level

Define how the NTR is sampled from the hierarchical Bayesian model (must be 0,1, or 2; see details)

correct.labeling

Labeling times have to be unique; usually execution is aborted, if this is not the case; if this is set to true, the median labeling time is assumed

verbose

Print status messages

Value

a new grandR object including a new analysis table. The columns of the new analysis table are

"s.A"

the posterior mean synthesis rate for sample A from the comparison

"s.B"

the posterior mean synthesis rate for sample B from the comparison

"HL.A"

the posterior mean RNA half-life for sample A from the comparison

"HL.B"

the posterior mean RNA half-life for sample B from the comparison

"s.log2FC"

the posterior mean synthesis rate log2 fold change

"s.cred.lower"

the lower CI boundary of the synthesis rate log2 fold change

"s.cred.upper"

the upper CI boundary of the synthesis rate log2 fold change

"s.ROPE"

the signed ROPE probability (negative means downregulation) for the synthesis rate fold change

"HL.log2FC"

the posterior mean half-life log2 fold change

"HL.cred.lower"

the lower CI boundary of the half-life log2 fold change

"HL.cred.upper"

the upper CI boundary of the half-life log2 fold change

"HL.ROPE"

the signed ROPE probability (negative means downregulation) for the half-life fold change

Details

The kinetic parameters s and d are computed using TransformSnapshot. For that, the sample either must be in steady state (this is the case if defined in the reference.columns matrix), or if the levels at an earlier time point are known from separate samples, so called temporal reference samples. Thus, if s and d are estimated for a set of samples x_1,...,x_k (that must be from the same time point t), we need to find (i) the corresponding temporal reference samples from time t0, and (ii) the time difference between t and t0.

The temporal reference samples are identified by the reference.columns matrix. This is a square matrix of logicals, rows and columns correspond to all samples and TRUE indicates that the row sample is a temporal reference of the columns sample. This time point is defined by time.experiment. If time.experiment is NULL, then the labeling time of the A or B samples is used (e.g. useful if labeling was started concomitantly with the perturbation, and the steady state samples are unperturbed samples).

By default, the hierarchical Bayesian model is estimated. If sample.level = 0, the NTRs are sampled from a beta distribution that approximates the mixture of betas from the replicate samples. If sample.level = 1, only the first level from the hierarchical model is sampled (corresponding to the uncertainty of estimating the biological variability). If sample.level = 2, the first and second levels are estimated (corresponding to the full hierarchical model).

if N is set to 0, then no sampling from the posterior is performed, but the transformed MAP estimates are returned

Examples

banp <- ReadGRAND(system.file("extdata", "BANP.tsv.gz", package = "grandR"),
          design=c("Cell","Experimental.time","Genotype",
                       Design$dur.4sU,Design$has.4sU,Design$Replicate))
contrasts <- GetContrasts(banp,contrast=c("Experimental.time.original","0h"),name.format="$A")
reference.columns <- FindReferences(banp,reference= Experimental.time==0)
banp <- EstimateRegulation(banp,"Regulation",
                             contrasts=contrasts,
                             reference.columns=reference.columns,
                             verbose=TRUE,
                             time.experiment = "Experimental.time",
                             N=0,               # don't sample in the example
                             dispersion=0.1)    # don't estimate dispersion in the example
#> Computing Regulation for 2h...
#> Sampling from non-steady state for mESC.2h.dTag.90min.4sU.A,mESC.2h.dTag.90min.4sU.B,mESC.2h.dTag.90min.4sU.C (reference: mESC.0h.dTag.2h.4sU.A,mESC.0h.dTag.2h.4sU.B)...
#> Sampling from steady state for mESC.0h.dTag.2h.4sU.A,mESC.0h.dTag.2h.4sU.B...
#> Beta prior for mESC.2h.dTag.90min.4sU.A,mESC.2h.dTag.90min.4sU.B,mESC.2h.dTag.90min.4sU.C: a=2.516, b=587.881
#> Beta prior for mESC.0h.dTag.2h.4sU.A,mESC.0h.dTag.2h.4sU.B: a=2.460, b=609.598
head(GetAnalysisTable(banp))
#>                       Gene  Symbol Length     Type Regulation.2h.s.A
#> Acsl3   ENSMUSG00000032883   Acsl3   3950 Cellular          583.4872
#> Pnkd    ENSMUSG00000026179    Pnkd   3599 Cellular          250.0747
#> Ankrd16 ENSMUSG00000047909 Ankrd16   3688 Cellular          863.1847
#> Dst     ENSMUSG00000026131     Dst  23279 Cellular         8054.4326
#> Vim     ENSMUSG00000026728     Vim   2926 Cellular          637.1811
#> Hdlbp   ENSMUSG00000034088   Hdlbp   6231 Cellular         3370.3347
#>         Regulation.2h.s.B Regulation.2h.HL.A Regulation.2h.HL.B
#> Acsl3            431.1494           3.768871           3.768871
#> Pnkd             181.5204           3.258101           3.258101
#> Ankrd16          722.3386           1.079148           1.079148
#> Dst             8175.6293           3.328200           3.328200
#> Vim              466.9153           4.845723           4.845723
#> Hdlbp           2615.4669           5.992583           5.992583
#>         Regulation.2h.s.log2FC Regulation.2h.s.cred.lower
#> Acsl3               0.43651305                       -Inf
#> Pnkd                0.46222685                       -Inf
#> Ankrd16             0.25699402                       -Inf
#> Dst                -0.02154682                       -Inf
#> Vim                 0.44854250                       -Inf
#> Hdlbp               0.36582333                       -Inf
#>         Regulation.2h.s.cred.upper Regulation.2h.s.ROPE Regulation.2h.HL.log2FC
#> Acsl3                         -Inf                   NA              -0.3138421
#> Pnkd                          -Inf                   NA              -0.3808574
#> Ankrd16                       -Inf                   NA              -0.3082710
#> Dst                           -Inf                   NA              -0.3483107
#> Vim                           -Inf                   NA              -0.9835942
#> Hdlbp                         -Inf                   NA              -0.5041620
#>         Regulation.2h.HL.cred.lower Regulation.2h.HL.cred.upper
#> Acsl3                          -Inf                         Inf
#> Pnkd                           -Inf                         Inf
#> Ankrd16                        -Inf                         Inf
#> Dst                            -Inf                         Inf
#> Vim                            -Inf                         Inf
#> Hdlbp                          -Inf                         Inf
#>         Regulation.2h.HL.ROPE
#> Acsl3                      NA
#> Pnkd                       NA
#> Ankrd16                    NA
#> Dst                        NA
#> Vim                        NA
#> Hdlbp                      NA