Normalize the number of cells per patient to achieve equal contribution of patients for Wilcoxon rank sum test I'm analyzing a single-cell RNAseq dataset where each patient has a few hundreds of sequenced cells before and after a medical intervention (blue and orange colors respectively in the graph below):

I want to test whethere the expression of a specific gene changes after the intervention.
The expression of this gene before and after the intervention is displayed in the next boxplot, where the Y axis is the expression score for all of the cells for every patient:

You can clearly see that in several patients, the cells reduce the median expression after the intervention while the median expression is higher after the intervention in other patients.
But I suspect that patient_8 (for example), which demonstrates a decrease of the expression, affects the total trend as it has lots of cells compared to the other patients, causing the entire data to demonstrate a decrease of this gene's expression when we compare all of the cells from all of the patients together before and after the intervention:

So overall, when I consider all of the cells from all of the patients together, it seems that this slight decrease that we see (which is in fact significant in the level of the single-cell) was actually seen because patient_8 showed a decrease of the expression and has lots of cells.
How can I correctly normalize the amounts of cells per patient so all patients will be considered equally while conducting wilcoxon rank sum test on all of the cells from all of the patients together?
I considered to randomally choose cells before and after the intervention from each patient, such that the number of cells will be equal to the minimal number of cells of the patient that had the lowest count of cells - But something in this thought feels wrong...
 A: A couple of warnings to start.
First, single-cell RNA sequencing (scRNAseq) typically finds its strength in distinguishing subclasses of cells within a mixture. It seems that you are assuming a single cell type in your analysis. Make sure that assumption is correct, or your results will at best be difficult to interpret cleanly.
Second, the Bioconductor project has well developed tools for such analysis. Those tools take many sources of variation into account, beyond what I suspect that you have included in your "Gene Expression Score" (which seems to be some normalized ratio of counts mapped onto the gene of interest divided by the total reads per cell) and probably beyond those that I'm familiar with as a non-expert. This sc-RNAseq course from the Wellcome Sanger Institute is a useful introduction.
The following outlines how to handle some major normalization issues with basic R, but you should consult the above references for current best practice.
The major issue is that you want to take the reliability of observations into account while not letting an individual case overwhelm less reliable cases inappropriately. You note this problem for the patients. You don't want to over-weight patient 8 by just looking at averages among all cells, as that might hide true differences among patients. But you also don't want to treat all patients equally, as patients with more cells do provide more information. You certainly don't want to throw away most of the data, as you would have to if you were to reduce the numbers to those of patients 4 and 5.
Those same considerations apply to the scRNAseq per-cell data. Cells with more total sequencing reads will provide more reliable estimates of the expression of your gene of interest. It's not clear whether that's taken into account by your "Gene Expression Score."
Both the among-cell and among-patient issues might be handled with a count-based generalized linear mixed model, for example with a Poisson model using the glmer() function in the R lme4 package. That would keep separate information for cells and for patients while weighting appropriately according to the reliability of the associated estimates.
For data, each cell would provide a separate row with an observation of counts mapped to the gene of interest, further labeled according to patient, pre/post intervention status, and total read counts for the cell. The modeling would thus take into account the reliability of each cell's gene-expression value, via the total read counts.
The pre/post-intervention status would be a fixed effect and patients would be treated as random effects. The mixed model provides an among-patient average pre/post-intervention difference, but with patients appropriately weighted according to their numbers of cells. A random intercept for patients would allow for differences among patients in baseline expression. An additional random "slope" would allow the intervention effect to differ among patients around an overall estimate of the effect.
The model with a random intercept might look something like this:
glmer(mappedCounts ~ interventionStatus + offset(log(totalCounts)) + (1|patient), family = poisson(link="log")) 

That models the log (the link function) of the mean number of the gene's mappedCounts as a function of interventionStatus. The offset term effectively places the result into a "per total read count" ratio scale, while the Poisson model recognizes that cells with more counts are more reliable in that ratio. The random intercept models Gaussian differences among patients in the baseline gene expression.
In practice, you might need to use a negative binomial model (often used in RNAseq) instead of a Poisson model, and to account for other sources of variability (like different sequencing runs among patients). But the principles are the same: weight cases, whether cells or patients, appropriately based on the reliability of the data they provide.
