Scanpy: kNN graph and leiden clustering differences between scanpy and Seurat

[mxposed]

In our analyses we wanted to try SCTransform normalization instead of default log-norm. I have done it quite crudely, but it works: I run SCT in Seurat and dump the counts on disk to load in scanpy.
While verifying that this approach worked, we encountered slight inconsistencies between clustering using (1) vanilla log-norm scanpy (2) SCT imported scanpy and (3) SCT in Seurat.
After investigation, it appears that vanilla scanpy sometimes better picks up some clusters than SCT+scanpy, despite the latter having more relevant genes in its HVG list.

Here is the investigation: https://github.com/mxposed/notebooks/blob/master/sct-scanpy.ipynb

And here are the main questions that remain:

1. Why Vanilla scanpy could resolve those populations, despite operating on less marker HVGs?
2. What is the difference between kNN graph construction and clustering between Seurat and scanpy?
3. How to be sure we did not undercluster and miss some smaller cell populations?

I would be glad for any feedback or input, and of course if someone knows the answers, that's great!

Best wishes,
Nick

PS. Thank you for scanpy!

[fidelram]

Thanks for the nice investigation. However, is very difficult to say what the differences could be. Have you checked the sc.tl.leiden parameters. I think that by default weights=True. Also sc.pp.neighbors has as default n_neighbors=15 ( I think). This may account for some differences.

[LuckyMD]

Hey @mxposed!

Nice work on the scArches covid model!

Love these types of comparisons :). There are a lot of potential differences. Just a few that come to mind:

1. Louvain and leiden are stochastic... random seeds will differ. This can generate quite different clusterings. I would check against different choices of random seed in Scanpy and/or Seurat to see what the background distribution is for clustering quality.
2. Seurat uses what it calls a "shared nearest neighbour graph". That is weighted by Jaccard Index by default. I'm not sure that's the case in scanpy.
3. Numerical library differences means you will never get the same result :/.

[dawe]

My two cents.
In additions to everything mentioned in this issue (parameters, rng, libraries) and the fact sNN and kNN are slightly different, neither Leiden nor Louvain method will save you from missing small populations (see here), it is a well-known issue when working with modularity.

[LuckyMD]

True, the resolution limit could play a role. But I would assume not yet at this scale, as we would expect a more even degree distribution in a kNN graph, and we're not looking at modules that make up a very small proportion of the total nodes. It would be interesting to check when this starts being important for larger atlases.

Have you found this to be a problem in some of your datasets?

[dawe]

@LuckyMD Mine was a general answer, I agree that at this scale it may not be an issue but it may indeed be for larger atlases.

@mxposed about your first question: cell distances depend on the HVG in absolute terms, but the overall structure of your data is more "relative". If the kNN graph topology is overall conserved you'll end up with similar populations.

[mxposed]

@fidelram Thank you for pointing this out. I did miss the n_neighbors parameter for sc.pp.neighbors function. It has default 15, while Seurat's have default 20 (and yes, they then prune the kNN). Adding this parameter did solve the discrepancy in the first dataset, but not the second, in the investigation (I have not updated the notebook).

@LuckyMD Thank you! I really liked applying scArches, and it's also a very natural approach: having a reference, mapping to it. I hope we're moving towards that direction generally.

Thank you for pointing out that Leiden is stochastic, I didn't realize that, and the fixed default random seed obscures it a little. I'll try to look at different seeds and assess the distribution of clustering. Can't estimate to which degree different runs would disagree.

Indeed, it appears that scanpy does kNN and doesn't do any pruning (judging from my brief glance at the code). I honestly expected that some kind of pruning of the kNN graph would be there. I remember two talks, one from Dana Pe'er and one from Dominic Grün, that mentioned kNN pruning as a strategy to improve analysis.

@dawe Thank you for linking to the resolution limit. However, I don't think it's the case here, because 2 of the 3 strategies that I tried did resolve those populations.

If we focus on dataset 2 (SC167) in the investigation, obviously, there's some small kNN topology difference between the strategies tried, that leads to SCT+scanpy strategy being slow to separate DC1 cells from B cells. I am mostly surprised that vanilla (log-norm) strategy does separate those cells.

I wonder how to go about investigating what drives that behaviour?