Our recent post elicited this comment from Bob that was so detailed I thought it deserved its own post.
This one comment from Bob has enough material for 4 standalone posts:
1. Fitting hierarchical models in genetics: why the full Bayesian approach gives the right answer, whereas an existing shortcut method using approximations and likelihood ratio tests does not.
2. Stan’s HMC scales well in high dimension, but not so well in bad geometry, so for this problem it was better to use a Poisson/gamma model with 400,000 latent parameters than a much lower dimensional but computationally awkward negative binomial model. (Recall the dictum that mixture models all have computational problems and that all computational problems are essentially mixture models.)
3. “Super-scalable penalized maximum likelihood inference for biome problems with some collaborators in optimization. With 150K biomes averaging 3 megabases each and arbitrary amounts of data, it fits in about 20m on a couple decent servers (doesn’t need GPU).”
4. The struggles of interdisciplinary, or interlab, collaboration within certain subfields of biology.
Bob’s summary of that last point: “It’s trivial to rewrite a statistically sound, fully Bayesian version of DESeq2, but biologists are too conservative to use new tools, so I haven’t even tried to fix this. In the end, I basically gave up working on biology because of the politics.”
And here’s Bob’s comment in full:
There are already hierarchical models in wide use in genomics. For example, rMATS is a system for inferring splice variants. The statistical computation in that package is a mess, which led them to the ridiculous conclusion that fewer variants from a population with deeper sequencing is better for estimating the population values than more variants with less deep sequencing (basic stats will show you the latter is better if your goal is estimating the population). I rewrote it all as a Bayesian model. It gets the right answer statistically, namely that it’s better to have more samples at lower sequence depth, unlike rMATS. The problem with rMATS is their Laplace approximation for max marginal likelihood and then their likelihood ratio tests, which together give them the wrong results statistically, even for their model.
Negative binomials induce really bad geometry, as does just about anything that gives you multiple ways to explain the same data (i.e., either higher mean or higher dispersion can explain a large value). So I reparameterized in two ways. First, I coded the differential expression with a gamma-Poisson rather than negative binomial, which with 20 sequencing runs over 20K splice variants led to 400K latent parameters. Second, I use an overall mean plus difference parameterization on Gelman’s advice. Stan’s HMC scales well in high dimension, but not so well in bad geometry, so fitting is OK.
For political reasons, I haven’t been able to publish. As an example, my original collaborator’s Ph.D. supervisor literally won’t let her work on anything other than his project and the person whose lab she rotated through is only interested if it produces “breakthrough biology” we can publish in Nature, which is not something I can deliver. I haven’t been able to find anyone who can help me evaluate this and write it up in a way that the biologists would be OK with. I’d be happy to share the code. I even hired an intern specifically to work on this, but the intern showed up and refused to work on the problem!
More recently, I worked on super-scalable penalized maximum likelihood inference for biome problems with some collaborators in optimization. With 150K biomes averaging 3 megabases each and arbitrary amounts of data, it fits in about 20m on a couple decent servers (doesn’t need GPU). It’s like Rob Patro et al.’s salmon, but I figured out how to do all the bias adjustments for things like hexamer and positional bias (just push the uncertainty through Bayes style—not hard), which is why Rob told me that he abandoned the approach. I’ve had the same problem getting a biologist to help evaluate and write up. Even my collaborators on the autism project where I have a ton of gut biome data from Simons Foundation’s Autism Research Initiative (SFARI) and we did publish a paper, but they say they have zero time to do this. Everyone says they know someone, but when I contact them, they’re too busy. If I could come up with half a year funding for a postdoc, they’d loan me one (they really do assign their postdocs work—it’s not at all like CS and stats). It’s amazing to me how constrained biologists are in working with other people. The code we produced was really great, because Robert Gower (an optimization specialist) helped add a mirror descent optimization algorithm, which is about an order of magnitude faster than the L-BFGS default in Stan. And Robert Blackwell (it’s the 3 Roberts project), a high-performace compute specialist, added some high performance sparse matrix computation in C++ and figured out how to scale it out on a (CPU) cluster. What’s cool is that time doesn’t depend on the data size other than the really fast conversion to k-mer counts.
The autism gut biome paper for which I helped with the Stan modeling is Multi-level analysis of the gut–brain axis shows autism spectrum disorder-associated molecular and microbial profiles. This project was led by Jamie Morton, who was an amazing postdoc at Flatiron. Jamie wrote a really nice overview of how to think about compositional analyses like differential expression in Establishing microbial composition measurement standards with reference frames. But he has zero time to help build a state-of-the-art differential expression tool.
I would suggest looking at DESeq2 to see what a mess current differential expression inference is. It’s pretty easy to see why none of their p-values are well calibrated. It’s trivial to rewrite a statistically sound, fully Bayesian version of DESeq2, but biologists are too conservative to use new tools, so I haven’t even tried to fix this.
In the end, I basically gave up working on biology because of the politics.
I am a biologist and I’ve got some cool RNASeq data to play with and other gene data as well. I’d totally love a Bayesian version of rMATS and DESeq2. Most of my PhD work was bringing Bayesian models and tools into ecology across different scales. These days it’s all about bringing Stan to solve problems for ecologists, bioinformatics folks, and managers.
Mostly an aside the point of the post, but what’s the intuition for why the gamma-Poisson has better geometry than the neg binomial? I know the fact that one can get a neg binomial from a gamma-Poisson mixture, but I’m curious about the computational differences.
I actually dont know but here’s some intuition in general.
The things that are hard, are staying within narrow channels, and wide range of scales. Imagine youre flying down the channel on the Death Star in Star Wars… The slightest lapse of control and you slam into the wall and die (divergence) but to get multiple pictures along the channel (samples) you have to fly in tight control for a long time.
Well, imagine that tight channel for Xbar is the result of the average of 100 different other variables… X_i Then you could just sample those X_I and calculate their average to get Xbar… Each of the X_I are easy to sample… Steps of size 1 maybe… So sampling the bigger model is less constrained and so easier to adapt to. Maybe like spaceships that individually dive down and take a picture of their little slice of the channel bit are otherwise very lightly constrained… The analogy is not perfect but maybe you get the idea.
There are at least two fields of Genetics (genetics is a vast field) that use Bayesian statistics. Mine is artificial selection (and evolution in ecology) in plant and animal breeding, which use hierarchical models (aka Mixed Models) since at least the 80’s. Essentially this theory makes predictions of an unobserved feature (how good is this animal/line as reproducer), and Bayesian theory is widespread, for prediction itself but also to derive new methods on top of old ones. There is a long tradition of teaching Bayesian in the community, e.g. “Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics” Sorensen & Gianola but also e.g. Appendixes in “Genetics and Analysis of Quantitative Traits”, Lynch and Walsh.
The other one is “Physiological Genomics” in which people try to understand the machinery from DNA to phenotype. Very often they need to resort to some sort of rank of genes, gene expression, levels, etc. and hypothesis testing is important. And it’s also important that co-workers understand what you do. That’s why (I think) biologists are conservative.