Skip to contents

Calculate logFCs and associated p-values for a given comparison, using either limma or the Negative Binomial quasi-likelihood framework of edgeR. The observed counts for the WT variants can be used as offsets in the model.

Usage

calculateRelativeFC(
  se,
  design,
  coef = NULL,
  contrast = NULL,
  WTrows = NULL,
  selAssay = "counts",
  pseudocount = 1,
  method = "edgeR",
  normMethod = ifelse(is.null(WTrows), "TMM", "sum")
)

Arguments

se

SummarizedExperiment object.

design

Design matrix. The rows of the design matrix must be in the same order as the columns in se.

coef

Coefficient(s) to test with edgeR or limma.

contrast

Numeric contrast to test with edgeR or limma.

WTrows

Vector of row names that will be used as the reference when calculating logFCs and statistics. If more than one value is provided, the sum of the corresponding counts is used to generate offsets. If NULL, offsets will be defined as the effective library sizes (using TMM normalization factors).

selAssay

Assay to select from se for the analysis.

pseudocount

Pseudocount to add when calculating log-fold changes.

method

Either 'edgeR' or 'limma'. If set to 'limma', voom is used to transform the counts and estimate observation weights before applying limma. In this case, the results also contain the standard errors of the logFCs.

normMethod

Character scalar indicating which normalization method should be used to calculate size factors. Should be either "TMM" or "csaw" when WTrows is NULL, and "geomean" or "sum" when WTrows is provided.

Value

A data.frame with output from the statistical testing framework (edgeR or limma).

Author

Charlotte Soneson, Michael Stadler

Examples

se <- readRDS(system.file("extdata", "GSE102901_cis_se.rds",
                          package = "mutscan"))[1:200, ]
design <- stats::model.matrix(~ Replicate + Condition,
                              data = SummarizedExperiment::colData(se))
                              
## Calculate "absolute" log-fold changes with edgeR
res <- calculateRelativeFC(se, design, coef = "Conditioncis_output", 
                           method = "edgeR")
head(res)
#>             logFC    logCPM         F       PValue          FDR logFC_shrunk
#> f.0.WT  1.9721334 20.083627 56.768150 2.160897e-10 3.601495e-09    1.9721331
#> f.1.AAC 1.3136230  4.688782 17.640288 8.439695e-05 2.910240e-04    1.2984544
#> f.1.AAG 1.0392549  5.677955 13.791731 4.313291e-04 1.327167e-03    1.0342409
#> f.1.ACC 2.0874024 12.862431 62.628899 4.638644e-11 1.159661e-09    2.0873503
#> f.1.ACG 1.8698694 13.513672 50.764347 1.134839e-09 1.194568e-08    1.8698442
#> f.1.AGC 0.9309018  4.754181  9.170734 3.545487e-03 8.975917e-03    0.9228638
#>         df.total df.prior df.test
#> f.0.WT  63.85553 60.85531       1
#> f.1.AAC 63.85462 60.85531       1
#> f.1.AAG 63.85524 60.85531       1
#> f.1.ACC 63.85553 60.85531       1
#> f.1.ACG 63.85553 60.85531       1
#> f.1.AGC 63.85487 60.85531       1
## Calculate log-fold changes relative to the WT sequence with edgeR
stopifnot("f.0.WT" %in% rownames(se))
res <- calculateRelativeFC(se, design, coef = "Conditioncis_output", 
                           method = "edgeR", WTrows = "f.0.WT")
head(res)
#>                 logFC    logCPM         F       PValue          FDR
#> f.0.WT   4.209817e-16 19.757846  0.000000 1.000000e+00 1.000000e+00
#> f.1.AAC -6.078796e-01  4.553595  9.224434 1.337186e-02 1.601421e-02
#> f.1.AAG -8.740695e-01  5.655807 24.839338 6.574158e-04 1.125187e-03
#> f.1.ACC  1.207258e-01 12.506900 28.957264 3.790068e-04 7.018644e-04
#> f.1.ACG -1.016593e-01 13.217924 67.771782 1.324404e-05 4.905202e-05
#> f.1.AGC -9.855150e-01  4.738963 18.559705 1.767705e-03 2.658202e-03
#>          logFC_shrunk df.total df.prior df.test
#> f.0.WT   4.802398e-16 9.426352 6.426305       1
#> f.1.AAC -6.661972e-01 9.423708 6.426305       1
#> f.1.AAG -8.823465e-01 9.425925 6.426305       1
#> f.1.ACC  1.207031e-01 9.426352 6.426305       1
#> f.1.ACG -1.016627e-01 9.426352 6.426305       1
#> f.1.AGC -1.030935e+00 9.424817 6.426305       1

## Calculate "absolute" log-fold changes with limma
res <- calculateRelativeFC(se, design, coef = "Conditioncis_output", 
                           method = "limma")
head(res)
#>             logFC      CI.L     CI.R   AveExpr        t      P.Value
#> f.0.WT  1.9474761 1.5292108 2.365741 19.757630 9.144191 9.237987e-19
#> f.1.AAC 1.2874854 0.6745388 1.900432  4.417300 4.125197 4.227413e-05
#> f.1.AAG 1.0494394 0.5382396 1.560639  5.581250 4.031733 6.250330e-05
#> f.1.ACC 2.0664575 1.6160089 2.516906 12.504065 9.009616 2.738680e-18
#> f.1.ACG 1.8444895 1.3981560 2.290823 13.217517 8.115994 2.737555e-15
#> f.1.AGC 0.9173399 0.3223443 1.512336  4.627474 3.027900 2.568326e-03
#>            adj.P.Val         B  se.logFC df.total df.prior
#> f.0.WT  3.079329e-17 34.491225 0.2129741      600      Inf
#> f.1.AAC 1.537241e-04  1.802578 0.3121028      600      Inf
#> f.1.AAG 2.232261e-04  1.302156 0.2602949      600      Inf
#> f.1.ACC 7.824800e-17 33.291805 0.2293613      600      Inf
#> f.1.ACG 3.421944e-14 25.725045 0.2272660      600      Inf
#> f.1.AGC 6.941422e-03 -2.063388 0.3029624      600      Inf
## Calculate log-fold changes relative to the WT sequence with limma
stopifnot("f.0.WT" %in% rownames(se))
res <- calculateRelativeFC(se, design, coef = "Conditioncis_output", 
                           method = "limma", WTrows = "f.0.WT")
head(res)
#>                 logFC        CI.L        CI.R   AveExpr             t
#> f.0.WT   9.552321e-08 -0.02134725  0.02134744 19.757630  1.061086e-05
#> f.1.AAC -6.050453e-01 -1.19863707 -0.01145359  4.417300 -2.417050e+00
#> f.1.AAG -8.921801e-01 -1.41471440 -0.36964582  5.581250 -4.048771e+00
#> f.1.ACC  1.233369e-01  0.05349146  0.19318235 12.504065  4.187362e+00
#> f.1.ACG -1.002978e-01 -0.12552201 -0.07507357 13.217517 -9.428860e+00
#> f.1.AGC -9.852809e-01 -1.66066239 -0.30989941  4.627474 -3.459365e+00
#>              P.Value    adj.P.Val           B    se.logFC df.total df.prior
#> f.0.WT  9.999918e-01 0.9999918338 -10.2853001 0.009002401 6.904162 3.904162
#> f.1.AAC 4.676895e-02 0.0534502279  -4.8181840 0.250323924 6.904162 3.904162
#> f.1.AAG 5.021241e-03 0.0079702233  -2.7275623 0.220358245 6.904162 3.904162
#> f.1.ACC 4.226773e-03 0.0071038198  -4.8399208 0.029454563 6.904162 3.904162
#> f.1.ACG 3.435219e-05 0.0002862682   0.1043135 0.010637319 6.904162 3.904162
#> f.1.AGC 1.078800e-02 0.0151943682  -3.3263517 0.284815526 6.904162 3.904162