3
$\begingroup$

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.: enter image description here

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): enter image description here

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!

Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – fra
    Feb 3, 2020 at 18:07
  • $\begingroup$ 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. $\endgroup$
    – Charles
    Feb 3, 2020 at 18:09
  • $\begingroup$ 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 $\endgroup$
    – alan ocallaghan
    Feb 6, 2020 at 10:32
  • 1
    $\begingroup$ This paper might be interesting to you biorxiv.org/content/10.1101/2020.01.15.906248v1 $\endgroup$
    – alan ocallaghan
    Feb 6, 2020 at 10:33
  • 1
    $\begingroup$ Thanks for the paper @alanocallaghan! $\endgroup$
    – Charles
    Feb 9, 2020 at 17:42

2 Answers 2

Reset to default
3
$\begingroup$

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/

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Charles
    Feb 9, 2020 at 17:52
  • $\begingroup$ 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. $\endgroup$
    – Andrew Judell
    Mar 2, 2020 at 15:07
  • $\begingroup$ @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. $\endgroup$
    – Andrew Judell
    Mar 2, 2020 at 15:22
2
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.