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.
