Best method to compare differentially expressed genes between 2 single cell clusters for GSEA

Question

I want to get a ranked gene list for GSEA analysis from single cell data, which needs a group differential expression statistic for ranking.

Which is the best method to compare the differently expressed genes between 2 clusters scanpy.tl.rank_genes_groups (https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.rank_genes_groups.html) or MAST.zlm (https://rdrr.io/bioc/MAST/man/zlm.html) or FindMarkers function of Seurat (https://satijalab.org/seurat/reference/findmarkers)?

Thanks

Answer 1

I wish to focus my contribution singly on ScanPy vs. Seurat FindMarkers.

ScanPy's claim is it is essentially a speeded up version of Seurat FindMarkers with better performance (discussed below) written in Python.

What they are doing are essentially datamining the expression signal using multivariate statistics (PCA) focused through tSNE. There is no controversy or dispute with this approach at all (please see Caveat I shouldn't use the term "datamining" in context - bad habit).

What ScanPy is doing is using Graph theory Louvian groups to extract the differential signal, presumably from the tSNE focused PCA. Cool.

If you've already got to grips with ScanPy then leveraging it as a data mining approach - to me - looks sensible. The statistics will not be hugely powerful, that doesn't mean they are bad, it approximately means they are solid.

What is ScanPy doing? The speeding up is via parallelisation, which is very easy in Python. The implication is Seurat is not parallelised, which would amaze me. The authors talk about some form of 'chunking' (they say 'partially loading the object') via HDF5. This is officially what they say

ANNDATA is similar to R’s EXPRESSIONSET [26], but supports sparse data and allows HDF5-based backing of ANNDATA objects on disk, a format independent of platform, framework, and language. This allows operating on an ANNDATA object without fully loading it into memory

I am confused by this, because HDF5 is a rapid I/O format - its not a datachunking format which his what the authors imply. Yes it facilitates datachunking but simply via rapid I/O. HDF5 is of course independent Python, but why discuss that at length (we all know that) when the critical information is how the partially loaded object is achieved. There's a lot of strategies to achieve the last bit - none of which are discussed.

The overall levels of high performance Python (outside ScanPy) - because the reality is Python's core code is slow - are three steps:

1. Parallelisation and code optimisation.
2. Data chunking - using algorithm strategies to partial load data e.g. Hadoop strategies.
3. Compiling via JIT, numba and static typing via Cython.

The authors have done point 1, however again HDF5 doesn't facilitate point 2 above, all it does is enable faster read/write. HDF5 is not actually the fastest I/O either both Pacquet (sp?) and feather are faster but normally for pandas dataframes. HDF5 is actually just as fast as pickle under many tests (!).

So the paper is a bit weird in that it is sparsely described - but this does often happen admittedly. The partial loading stuff via HDF5 is weird, because they are not describing how that is done.

The other bit that is weird is the ML and ANNDATA class. They state,

ANNDATA is similar to R’s EXPRESSIONSET [26]

I rarely use R so to me that isn't clear. The other weird thing is the authors talk about ML, they use an ML library, but give no further details nor why it is used. What I would assume is either: that they use it to assist with pseudotemporal ordering analysis OR they use it for unsupervised learning (PCA, tSNE). The latter is more likely. Louvian groups are from the igraph package they use

So in my personal view, despite the paper being a bit (to very) obscure I see no reason why ScanPy cannot replace Seurat purely in terms of the datamining being performed. Python is just as good at PCA and tSNE etc ... The use of Louvain groups is cool to trap that information. The key benefit of ScanPy is:

• your already familiar with it.
• faster - it definitely could be accelerated further.
• good for Python

The drawbacks I can see are:

• RNA analysis is dominated by R and you'll likely pick up a R focused reviewer (which could be fun).
• As per your previous post there appear shortcomings with ScanPy - but like the development team can't think of everything first time round.
• ScanPy2 is not released.
• its not generating a critical (test) probability to distinguish between expression levels just like Seurat - there are real advantages to this however.

Caveat I keep saying "data-mining" and I shouldn't. If the reviewers are perfectly happy with the above analysis alone, that's cool. This is because what I repeatedly call "data-mining" are other investigators final results, particularly in RNA expression. If your aim is to publish a big RNA-seq data set and identify the core features of e.g. a cellular immune response - good approach. I don't generate data, mine is purely analytics and algorithms, so my perspective (and reviewing standards) are different.

---

Just to state my interest in DEGs (or DGE). I use unsupervised learning, PCA, tSNE to generate a hypothesis and get a feel for the data. Then however, I do customised full-scale ML (supervised learning) to interrogate the data. Finally, I use DEG to assess which way around the ML result is pointing (under- or over-represented). I will not let the DEG override the ML (supervised learning) analysis - which I see as much stronger analysis. My context is different however and I'm consolidating the result is via extensive custom databases, i.e. very different.