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)?


  • $\begingroup$ Just to add I personally think this is a cool question, others may disagree. There is quite a bit of discussion on DEGs per se ... ScanPy hits that R to Python threshold, which is no trivial issue IMO. $\endgroup$
    – M__
    Jan 2 at 4:17

1 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.

  • $\begingroup$ In the ScanPy paper Wolf et al 2018, when the authors say "HDF5" I wonder if they mean "HDFS"? This would be the Hadoop Distributed File System and would actually make sense. I can't imagine it - it's a HUGE error - but the rationale with partial object loading would be completely understandable. It does demonstrate how sparsely this paper is written some algorithms in Genome Biology papers are so detailed it's incredible. This is the opposite spectrum. $\endgroup$
    – M__
    Jan 2 at 13:06
  • $\begingroup$ Additional weird point The authors calling a class ANNdata for machine learning is also weird. ANN stands for Artificial Neural Network and is a basic type of deep learning. The authors definitely do not use deep learning in ScanPy ... their analysis appear to be singly unsupervised learning. It's a weirdly written paper. I guess the authors must be carrying kudos in algorithm development - that does happen. The reviewers can wave through a paper if the authorship is carries kudos and editor doesn't pick it up. $\endgroup$
    – M__
    Jan 2 at 14:06

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.