Merging mashr analyses

merge_mashr.R

library(dreamlet)
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)

# Differential expression analysis within each assay,
# evaluated on the voom normalized data
res.dl <- dreamlet(res.proc, ~group_id)

# 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)

# same results from feeding in probabilities directly
prob.v2 <- compositePosteriorTest(1 - get_lfsr(res_mash$model) , include, exclude)

identical(prob, prob.v2)

# posteriors from first dataset
post1 = 1 - get_lfsr(res_mash$model)
colnames(post1) = paste0('test1_', colnames(post1))

# posteriors from second dataset
# (pretend this is a second dataset)
post2 = 1 - get_lfsr(res_mash$model)
colnames(post2) = paste0('test2_', colnames(post2))

post_merge = cbind(post1, post2)

# define new include and exclude
prob_merge <- compositePosteriorTest( post_merge, include, exclude)