Perform composite test on results from mashr

Perform composite test evaluating the specificity of an effect. Evalute the posterior probability that an a non-zero effect present in _all_ or _at least one_ condition in the inclusion set, but _no conditions_ in the exclusion set.

Usage

compositePosteriorTest(
  x,
  include,
  exclude = NULL,
  test = c("at least 1", "all")
)

Arguments

x

"dreamlet_mash_result" from run_mash()

include

array of conditions in the inclusion set

exclude

array of conditions in the exclusion set. Defaults to NULL for no exclusion

test

evaluate the posterior probability of a non-zero effect in "at least 1" or "all" conditions

Details

The posterior probabilities for all genes and conditions is obtained as 1-lFSR. Let prob be an array storing results for one gene. The probability that _no_ conditions in the exclusion set are non-zero is prod(1 - prob[exclude]). The probability that _all_ conditions in the inclusion set are non-zero is prod(prob[include]). The probability that _at least one_ condition in the inclusion set is non-zero is 1 - prod(1 - prob[include]). The composite test is the product of the probabilties computed from the inclusion and exclusion sets.

See description in section Identifying shared and cell type specific genetic regulatory effects of Zeng, et al. (https://doi.org/10.1101/2024.11.02.24316590).

See also

run_mash()

Examples

library(muscat)
library(mashr)
library(SingleCellExperiment)

data(example_sce)

# create pseudobulk for each sample and cell cluster
pb <- aggregateToPseudoBulk(example_sce[1:100, ],
  assay = "counts",
  cluster_id = "cluster_id",
  sample_id = "sample_id",
  verbose = FALSE
)

# voom-style normalization
res.proc <- processAssays(pb, ~group_id)
#>   B cells...
#> 0.016 secs
#>   CD14+ Monocytes...
#> 0.029 secs
#>   CD4 T cells...
#> 0.017 secs
#>   CD8 T cells...
#> 0.016 secs
#>   FCGR3A+ Monocytes...
#> 0.016 secs

# Differential expression analysis within each assay,
# evaluated on the voom normalized data
res.dl <- dreamlet(res.proc, ~group_id)
#>   B cells...
#> 0.01 secs
#>   CD14+ Monocytes...
#> 0.009 secs
#>   CD4 T cells...
#> 0.0083 secs
#>   CD8 T cells...
#> 0.0063 secs
#>   FCGR3A+ Monocytes...
#> 0.0087 secs

# run MASH model
# This can take 10s of minutes on real data
# This small datasets should take ~30s
res_mash <- run_mash(res.dl, "group_idstim")

# Composite test based on posterior probabilities
# to identify effect present in *at least 1* monocyte type
# and *NO* T-cell type.
include <- c("CD14+ Monocytes", "FCGR3A+ Monocytes")
exclude <- c("CD4 T cells", "CD8 T cells")

# Perform composite test
prob <- compositePosteriorTest(res_mash, include, exclude)

# examine the lFSR for top gene
get_lfsr(res_mash$model)[which.max(prob), , drop = FALSE]
#>       B cells CD14+ Monocytes CD4 T cells CD8 T cells FCGR3A+ Monocytes
#> FBXO6      NA              NA          NA          NA      1.981561e-17

# Test if *all* cell types have non-zero effect
prob <- compositePosteriorTest(res_mash, assayNames(res.dl))