Hi, wondering if anyone has thought about the best way to do single cell DE testing including subject as a random effect in a multilevel design including many subjects contrasting timpoint and disease state effects within a single "cluster"

For example, consider data similar to the design in the Limma guide section 9.7, except observations are cells, there are n cells per timepoint for each subject. I'm only looking at cells from a single cell type ("cluster" if you like), within that cluster I have many subjects with 2 timepoints and 2 states, say, diseased and healthy.

If I convert my SingleCellExperiment object into a DGElist, using y = convertTo(sce, type = "edgeR") I can test the contrast of e.g. the delta between timepoint 1 and 2 in diseased vs the delta between timepoint 1 and 2 in healthy subjects without including the intra-subject correlation in the model as a random effect e.g.

subject = c("sub1", "sub1", "sub2", "sub2", "sub3", "sub3",
             "sub4", "sub4", "sub5", "sub5", "sub6", "sub6") 
condition = c(rep("diseased", 6), rep("normal", 6)) 
timepoint = c(rep(c("1","2"),6)) 

model = cbind(subject,condition,timepoint) %>% 
  as.data.frame() %>% 
  mutate_if(is.character, as.factor) %>% 
  unite("condition_time", c("condition", "timepoint"))

ct.matrix = model.matrix(~0 + model$condition_time + model$subject)

y = convertTo(sce, type = "edgeR")
y$samples = cbind(y$samples, model)

v <- voom(y, ct.matrix)
fit <- lmFit(v, ct.matrix)

contrast_matrix  = 
    makeContrasts(t1_disease_vs_normal = 
                 (condition_timediseased_2 - condition_timediseased_1) -condition_timenormal_2 - condition_timenormal_1))

c_fit = contrasts.fit(fit = fit, contrasts = contrast_matrix)
eb_c_fit = eBayes(c_fit)

In this case though, the sce object has many replicates for each subject-timepoint combination where the replicates are cells of a single cell type.

Question 1 It seems that it's not possible to add subject as a random effect in edgeR https://support.bioconductor.org/p/87292/
So I would have to embed single cell data into an elist object...and not use voom / edgeR in order to accomodate subject as a random effect? For example I would have to create an list "y" with single cells as observations and use :

corfit = duplicateCorrelation(y, response_time_matrix, block=design$subject)
corfit$consensus

I'm not sure this is advisable as I'd like to include the size factors calculated by scran as normalization factors for use in voom. One alternate idea is to summarize the scran normalized library by averaging each subject's cells at a given timepoint. Thus you would get a "meta-cell" from cluster k, subject j, timepoint t, of size n after normalizing out cell-specific differences like cell size, capture efficiency differences between cells with scran... this is a slightly different philosophy than using a "pseudobulk" library pool (as in section 6 of the DE testing vignette in the SimpleSingleCell workflow) of all cells as it avoids overrepresentation biases from e.g. large cells in the sample pool.

Question 2 Even if I did this, is there a way I could use such averaged summarized data into a linear model with random effects in limma by embedding the data into a elist object for limma?

This level of model complexity is uncommon in single cell analysis but I thought it would be interesting to look at.

thanks for reading.
-matt

Question 1: With respect to duplicateCorrelation, we explored this approach (among others) some time ago. I would say that it wasn't worth it. Running duplicateCorrelation() across hundreds of samples takes a long time without providing clear benefits with respect to type I error control.

With respect to making your suggested "meta-cells": it is debatable whether or not you want to normalize before summing. On one hand, it ensures that each cell contributes equally to the resulting meta-cell; but on the other hand, large cells already contribute more RNA in real bulk libraries, and that doesn't cause any real dramas in interpretation. I can't say one is better than the other, only that they are different. What I will say, though, is that summing normalized counts has less predictable statistical properties than summing the raw counts. The latter sum still is a count, has a lower bound of Poisson variance, and can be reasonably modeled with a NB distribution; while anything goes with the former.

For the record, there is a third option of weighting the pseudo-bulk samples by the number of cells used in its summation. These weights can, in theory, be used in edgeR's GLM fitting to account for differences in precision for each pseudo-bulk sample. However, weighting assumes that all cells were i.i.d. NB distributed within each group, which is obviously untrue. As the number of cells increase, the inter-sample variation dominates over the imprecision of each sample, and the violations of the i.i.d. assumption become more obvious.

Question 2: Well, you could just pretend it's a matrix from bulk RNA-seq, and get limma to treat it as such. Don't take this as a recommendation, though.