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