Hi everyone,

I’m working on a scRNA-seq dataset with multiple biological samples. After standard normalization/integration, PCA, and clustering, I identified several cell populations.

Now, I’m interested in one gene, G. Within some cell populations, some cells express G (G+) while others do not (G-). I would like to test, within each cell population, whether there are DEGs between G+ and G- cells.

My plan is to create pseudobulk profiles by aggregating raw counts per:

sample_id × cell_population × G_status

Then, for each cell population, run DESeq2 comparing G+ vs G-.

My question is about the design formula. Many samples contribute both G+ and G- cells, so a blocked/paired design seems appropriate:

design = ~ sample_id + G_status

However, the design is not perfectly paired: for some cell populations, a few samples only have G- cells because they have no G+ cells.

I tried both:

design = ~ G_status

and

design = ~ sample_id + G_status

and for some populations the results are completely different. The unblocked model gives thousands of DEGs, while the sample-blocked model gives almost none, sometimes only G itself, even when most samples contribute to both groups.

So my main question is:

For pseudobulk DE comparing G+ vs G- cells within the same cell population, should I include sample_id as a covariate even when the design is only partially paired? Or can this overcorrect and remove meaningful signal?

Any advice on the best design, or on whether this comparison is conceptually appropriate, would be very appreciated.