# Differential gene expression bias due to effect of an individual sample

I am analysing a human single cell RNA seq experiment, where we have 4 groups, four samples each. Data has been analysed using Seurat, with the canonical workflow. I have tried DE using various methods (Wilcox, MAST, LR), and for each of them the top genes (biggest fold change, lowest adjusted p) are genes over expressed by specific samples, e.g.:

Looking at the average expression for these genes per each subject, and plotting on a Boxplot, it does not seem to be a genuine difference (e.g. for the gene above):

Besides the subject, there is no other parameter (e.g. age, sex) that is biased across groups and that shows up as differentially represented in the PCA plot.

Is there a DE method that takes into account this behaviour and corrects for it? How can these results be presented? Besides increasing number of subjects (not really feasible), is there a solution to analyse these data?

Thanks a lot!

• How did you scale the data (ScaleData or SCTransform)? You can usually regress out any source of variation in your data using vars.to.regress = "Group". In your case, "Group" would be the column of object@metadata containing the Group information. – fra Feb 3 at 18:07
• Thanks for the reply, I did not do SCTransform, and with ScaleData I did not regress any variable. I do not think that scaling by "Group" will be useful: 1) I do not want to mask genuine changes between groups, 2) the DE methods use normalised data, not scaled data. – Charles Feb 3 at 18:09
• Regressing out group would eliminate most if not all differential expression between groups. Regressing out sample would be a slightly more reasonable approach but would likely end up being pretty similar – alan ocallaghan Feb 6 at 10:32
• This paper might be interesting to you biorxiv.org/content/10.1101/2020.01.15.906248v1 – alan ocallaghan Feb 6 at 10:33
• Thanks for the paper @alanocallaghan! – Charles Feb 9 at 17:42

Leave-one-out cross validation is an iterative method that leaves out one sample until each sample has been left out once. You can then identify the intersection of the differentially expressed genes from the sample space of repeated analyses to produce a set of genes that are consistently differentially expressed independent of individual sample biases. more info here: http://efavdb.com/leave-one-out-cross-validation/

• Thanks a lot Drew J., this seems a valid approach! However, I do wonder why there is not a DE method addressing this issue directly. I am not a statistician, so I may have missed some theory behind this, though. – Charles Feb 9 at 17:52
• The statistical models generally used for DE involve population level assumptions that are applied to incredibly small sample sizes. The main problem is that sequencing is still prohibitively expensive for small labs to produce data from experiments that include sufficient numbers to address effect size and power. – Drew J. Mar 2 at 15:07
• @Charles if you found the answer useful, please up-vote it. I am trying to build up my reputation points on this forum specifically to enable voting. Thank you. – Drew J. Mar 2 at 15:22

One potential problem here is that, as far as I can see, you're not taking the hierarchical structure of your scRNA-seq data into account. I.e. the fact that cells from the same patient / sample will be more similar than cells from different samples, even within the same group. This is also the problem addressed in the paper mentioned by @alan. Also, it usually makes more sense to do this kind of group vs. group comparison for each subpopulation (cluster/cell type) separately (assuming you have multiple subpopulations in your data and did some kind of clustering).

One way of doing this is by using a pseudobulk-level approach, where you basically aggregate (by either averaging or summing) your data over each sample, for each cluster separately. So you end up with one counts matrix per subpopulation with one row for each gene and one column for each sample. The sample columns are then the sums (or averages) of all cells belonging to that cluster-sample combination.

You can then use bulk-RNAseq methods such as edgeR or DESeq2 on the aggregated data to test differences between groups of samples for each cluster. The muscat package implements all these steps and has quite extensive documentation. So you could check that and test it out on your data.

The paper describing the method of muscat can be found here: https://www.biorxiv.org/content/10.1101/713412v2