Scanpy: Introduce GLM-PCA to scanpy.external

[flying-sheep]

I agree with @adrianveres here: GLM-PCA by @willtownes seems like a rigorous way to deal with zero inflation, since we use PCA as base for all subsequent processing. It’s apparently much faster than e.g. ZINB-WaVE

• twitter thread explaining the problem and solution
• code

[fidelram]

@flying-sheep I think that the code is here: https://github.com/willtownes/glmpca

[willtownes]

Hi to all, thanks for your interest in glmpca. I have been thinking of doing a python package now that the R package is finished and it would be an honor to have it included in scanpy. Can you give me a sense of how urgently you would need the package (ie what is the typical release cycle)? Also let me note a few caveats about the method:

• It does not handle zero inflation (which ZINB-WAVE does). However, we argue in our paper that despite large numbers of zeros, UMI data are not zero-inflated. We do not make any claim about the appropriateness of the glmpca model for non-UMI data (eg Smart-Seq read counts), which may actually be zero-inflated, although you could certainly run it with eg the negative binomial likelihood.
• glmpca is an alternative to PCA but not necessarily a replacement to PCA. For example, it is at least 10x slower than PCA and we are still working on the big data implementation for sparse matrices (in other words, we assume you can load the data matrix in dense form, which can be limiting).
• We describe a fast approximation to GLM-PCA in the paper which involves transforming raw counts to either Pearson or deviance residuals from a null model then applying standard PCA to that. This approach is just as fast as PCA as long as the null model can be computed in closed-form, which is what we have implemented here: https://github.com/willtownes/scrna2019/blob/master/util/functions.R#L164 . The idea is similar to the sctransform approach used by seurat, but the computation is simpler and faster.
• We also provide a deviance-based gene filtering method which is an alternative to using highly variable genes. This and the residuals functions will be available as an R package on bioconductor.

I look forward to collaborating with you all to help make these methods available to a wider community!

[willtownes]

Oops sorry here is the correct line for the null residuals function, please ignore the link in previous comment. https://github.com/willtownes/scrna2019/blob/master/util/functions.R#L206

[flying-sheep]

OK, seems like I misunderstood the point about zero inflation here. You just meant “large number of zeroes” as in “pretty sparse” then?

A factor of 10 isn’t that bad for something that’s more complex, and I doubt PCA speed is the bottleneck for most datasets.

So not a replacement, but an enhancement. As such, it would probably live in scanpy.external except if you want to develop it within scanpy instead of as a separate package (which is possible, but would tie you to our – currently slow but we’ll get better) release cycle.

[willtownes]

A rough implementation of glmpca in python is now available here: https://github.com/willtownes/glmpca-py . I will try to get it organized as an installable package tomorrow and add unit tests. Issues/ pull requests welcome.

[willtownes]

The package is now available on pypi and there is an automated test suite.

[flying-sheep]

Great! I’m sure this will soon find its way into scanpy.external!

[flying-sheep]

@falexwolf @willtownes @LuckyMD Valentine Svensson suggests that zero inflation does not exist in droplet protocols, but that log-transforming data could be responsible for the apparent zero inflation. Further, the high number of zeros can be accurately modeled with a non-zero-inflated model: https://www.nature.com/articles/s41587-019-0379-5

Since GLM-PCA doesn’t model zero inflation, it’s probably a really good base for distance calculations in scanpy in cases where its performance is sufficient. From the paper:

The multinomial model adequately describes negative control data, and there is no need to model zero inflation.

[davidsebfischer]

I agree that there is definitely no need to model zero inflation.

[flying-sheep]

My point is: The paper suggests that the reason for the zero inflation idea might be the log-transform, so we should offer a better path from counts to distances.

Our current path is to go via normally distributed “expression” values which can be used to calculate distances. Something like fold changes, i.e. what we currently do by log-transforming (because log(count) = foldchange):

counts (→ normalization) → expressions (→ normalization) (→ embedding) → distances → analyses

We use PCA as a latent space embedding here for efficiency purposes but it’s not required, we could calculate distances directly from expressions.
One alternative is to stay with (possibly bias-normalized) counts, and create our latent space from those directly using a suitable model (like GLM-PCA):

counts (→ normalization) → embedding → distances → analyses

The other alternative is to offer something like SCTransform and stay with our original path, bu better.

---

All of this is of course only super important if the log transform turns out to be maximally problematic (which the amount of successful data analyses using it doesn’t suggest), but I think offering alternatives will definitely be very beneficial!

[willtownes]

Hi, thanks again for your interest in GLM-PCA. We welcome its inclusion in scanpy, but some caveats are that it is about 10x slower than PCA and we are still working to improve its numerical stability and ability to handle sparse data matrices.

With that in mind, we have put together an implementation of Pearson and deviance residuals as an approximation to GLM-PCA via the scry R package. These residuals, based on binomial and poisson approximation to multinomial, can be computed in closed form so they are computationally as fast as log-transforming. The sctransform method uses a negative binomial likelihood which doesn't have a closed form solution and is more complicated to implement (although we do recomment it from a statistical validity standpoint).

In addition to the null residuals, the scry package has an implementation of feature selection via deviance, which may also be of interest as an alternative to highly variable genes. This is also a closed form computation. Both the feature selection and null residuals functions allow adjusting for categorical batch labels.

I do hope to implement both of these in python eventually but it's pretty far down my to-do list. Given the functions are fairly simple, I welcome anyone to go ahead and copy them into python if they find it potentially useful.

[flying-sheep]

Sounds good! Maybe there’s takers in @theislab

[davidsebfischer]

We can compute any kind of residual in principle (NB right now, Poisson + log normal soon) in batchglm, this is also not restricted to categorical covariates and automatically selects whether closed form is available. The implmentation is numerically very stable. I think it could be a potential backend for this, @willtownes?

[willtownes]

@davidsebfischer that sounds very appealing, I will probably use it whenever I do the python implementation. Thanks!