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.
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
sefor 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"whenWTrowsisNULL, and"geomean"or"sum"whenWTrowsis provided.
Examples
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: ‘MatrixGenerics’
#> The following objects are masked from ‘package:matrixStats’:
#>
#> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#> colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#> colWeightedMeans, colWeightedMedians, colWeightedSds,
#> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#> rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#> rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: ‘generics’
#> The following objects are masked from ‘package:base’:
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: ‘BiocGenerics’
#> The following objects are masked from ‘package:stats’:
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from ‘package:base’:
#>
#> Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#> mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
#> unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: ‘S4Vectors’
#> The following object is masked from ‘package:utils’:
#>
#> findMatches
#> The following objects are masked from ‘package:base’:
#>
#> I, expand.grid, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#>
#> Attaching package: ‘Biobase’
#> The following object is masked from ‘package:MatrixGenerics’:
#>
#> rowMedians
#> The following objects are masked from ‘package:matrixStats’:
#>
#> anyMissing, rowMedians
se <- readRDS(system.file("extdata", "GSE102901_cis_se.rds",
package = "mutscan"))[1:200, ]
design <- model.matrix(~ Replicate + Condition,
data = 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 3.461728e-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 5.181525e-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