Skip to contents

Extract a table differential expression results from models fit with lucida()

Usage

results(fit, coef = NULL, method = c("classic", "Bayesian"), expand = FALSE)

# S4 method for class 'lucida'
results(fit, coef = NULL, method = c("classic", "Bayesian"), expand = FALSE)

Arguments

fit

object from lucida()

coef

specify which coefficient, or combination of coefficients to test.

method

"classic": multiple testing correction using Benjamini–Hochberg method, or "Bayesian": shrink coefficient estimates and report local FDR from ashr::ash() using empirical Bayes approach

expand

include sigSq_g and theta for each gene

Value

tibble storing hypothesis test for each gene and cell type

Details

The coef argument can be 1) a single coef name to test one coefficient: "x1", 2) an array coef names to test multiple coefficients: c("x1", "x2"), 3) string of linear functions of coefficients: "x1 - x2".

Multiple tesing correction and shrinkage are applied across all cell types as the same time. This gives a 'study-wide' rather than a `cell type level` multiple testing correction.

"classic": logFC is maximum likelihood estimate, se is frequentist standard error, and FDR is computed using Benjamini–Hochberg method.

"Bayesian": uses empircal Bayes shrinkage and multiple testing method of Stephens (2017) implemented in ashr::ash(). logFC is the posterior mean shrinkage estimate, se is the posterior standard deviation, lFDR is the local false discovery rate, and lFSR is the local false sign rate.

References

Stephens, M. (2017). False discovery rates: a new deal. Biostatistics, 18(2), 275-294. doi:10.1093/biostatistics/kxw041

Examples

library(SingleCellExperiment)

# Load example data
data(example_sce, package="muscat")
sce <- example_sce

# Compute library size for each cell
sce$libSize <- colSums(counts(sce))

# Specify regression formula and cell annotation 
form <- ~ group_id 
fit <- lucida(sce, form, "cluster_id", verbose=FALSE)
#> B cells 
#> CD14+ Monocytes 
#> CD4 T cells 
#> CD8 T cells 
#> FCGR3A+ Monocytes 

# Differential expression results
# classic
results(fit, "group_idstim")
#> # A tibble: 4,841 × 7
#>    cluster_id        ID                AveExpr  logFC     se   P.Value       FDR
#>    <chr>             <chr>               <dbl>  <dbl>  <dbl>     <dbl>     <dbl>
#>  1 FCGR3A+ Monocytes ISG15               40.8   5.02  0.182  2.69e-166 1.30e-162
#>  2 CD14+ Monocytes   ISG15               61.2   5.39  0.197  3.17e-165 7.67e-162
#>  3 CD14+ Monocytes   ISG20                9.82  3.96  0.160  5.84e-135 9.43e-132
#>  4 CD14+ Monocytes   APOBEC3A_ENSG000…   14.1   4.44  0.184  1.70e-128 2.06e-125
#>  5 FCGR3A+ Monocytes ISG20               10.1   3.44  0.154  9.96e-110 9.64e-107
#>  6 CD14+ Monocytes   TMSB10               8.79  1.29  0.0581 1.25e-109 1.01e-106
#>  7 FCGR3A+ Monocytes IFITM3               9.23  2.88  0.147  3.04e- 85 2.10e- 82
#>  8 B cells           ISG20                5.79  2.47  0.127  5.55e- 84 3.36e- 81
#>  9 CD14+ Monocytes   DYNLT1               3.34  1.57  0.0816 1.14e- 82 6.11e- 80
#> 10 CD14+ Monocytes   PFN1                10.9  -0.957 0.0506 1.26e- 79 6.10e- 77
#> # ℹ 4,831 more rows

# Bayesian
results(fit, "group_idstim", method="Bayesian")
#> # A tibble: 4,841 × 8
#>    cluster_id        ID      AveExpr  logFC     se   p.value      lFDR      lFSR
#>    <chr>             <chr>     <dbl>  <dbl>  <dbl>     <dbl>     <dbl>     <dbl>
#>  1 FCGR3A+ Monocytes ISG15     40.8   5.01  0.194  2.69e-166 2.37e-162 2.37e-162
#>  2 CD14+ Monocytes   ISG15     61.2   5.39  0.197  3.17e-165 2.62e-161 2.62e-161
#>  3 CD14+ Monocytes   ISG20      9.82  3.96  0.160  5.84e-135 8.20e-132 8.20e-132
#>  4 CD14+ Monocytes   APOBEC…   14.1   4.36  0.146  1.70e-128 2.85e-125 2.85e-125
#>  5 CD14+ Monocytes   TMSB10     8.79  1.29  0.0591 1.25e-109 5.28e-107 5.28e-107
#>  6 FCGR3A+ Monocytes ISG20     10.1   3.42  0.167  9.96e-110 1.23e-106 1.23e-106
#>  7 FCGR3A+ Monocytes IFITM3     9.23  2.87  0.144  3.04e- 85 1.79e- 82 1.79e- 82
#>  8 B cells           ISG20      5.79  2.46  0.136  5.55e- 84 3.55e- 81 3.55e- 81
#>  9 CD14+ Monocytes   DYNLT1     3.34  1.55  0.0762 1.14e- 82 3.69e- 80 3.69e- 80
#> 10 CD14+ Monocytes   PFN1      10.9  -0.957 0.0508 1.26e- 79 1.91e- 77 0        
#> # ℹ 4,831 more rows