Also estimate an innovation rate to deal with empty families at the root

[tomwenseleers]

I was wondering if instead of fitting a birth-death model it would be possible to allow for fitting a birth-death-innovation model to also explicitly allow for gene families with zero members at the root. I keep on running into problems where all the gene families that are of interest of me get kicked out because of being empty at the root. I know about the option --zero_root but that doesn't seem to fundamentally solve the problem and in the documentation is explicitly advised against. Or if an extra innovation rate parameter cannot be introduced, would you have suggestions for alternatives that can fit such models or deal better with empty families at the root? Maybe https://github.com/csurosm/count or https://github.com/ssolo/recount-public?

[benfulton]

CAFE works with explicit transition probabilities for its numbers, and I don't think the model you suggest would give us those.

[hahnlab-user]

I would say that either of the methods you suggest would work well, though they are giving slightly different kinds of answers. You could also try just analyzing sub-trees of your species tree in which these families do exist at the root, if you want family-specific inferences (like ancestral states).

[tomwenseleers]

Thanks both for the helpful responses. Analysing subtrees is what I was doing already, but even when working with such subtrees I still have a significant number of zero-root families, many of which with interesting biological functions. So that problem remains unfortunately... The specific problem I was analysing was looking at gene expansions and contractions across several independent origins of sociality in insects (the final analysis will have 250+ species). And so the smallest set I can prune my tree to would be a single origin of sociality plus some closely related solitary species. But even then many zero-root families remain, including many that are highly relevant to the evolution or elaboration of sociality (odorant receptor families etc).

I should clarify that by innovation I meant a copy-independent gain/immigration term, with transition rates

$$q_{n,n+1}=\lambda n+\nu,\qquad q_{n,n-1}=\mu n,$$

rather than an additional per-copy duplication rate.

My understanding is that this standard linear birth–death–immigration process does have explicit transition probabilities. For example, even with equal duplication and loss rates, starting from zero gives a negative-binomial count distribution, with

$$P_{00}(t)=(1+\lambda t)^{-\nu/\lambda}.$$

Transition probabilities for other starting counts can also be computed explicitly using generating functions and positive recurrences. This is an established model class, including in the work underlying [Count](https://doi.org/10.1093/bioinformatics/btq315) and the [gain–loss–duplication likelihood framework](https://arxiv.org/abs/0902.0970). Am I overlooking a restriction specific to CAFE’s implementation, rather than a lack of transition probabilities for this process?

I now have a [research fork with an optional BDI implementation](https://github.com/tomwenseleers/CAFE5/tree/innovation-bdi), with the assumptions and calculations described in the [mathematical documentation](https://github.com/tomwenseleers/CAFE5/blob/innovation-bdi/docs/innovation_mathematics.md). Aside from allowing for a copy-independent innovation (gain) rate, it also allows for asymmetrical birth & death rates, as well as for gamma rate variation, observation error, and separate refitted-bootstrap significance calculations. I recognise that these additional components require their own assessment, particularly regarding root assumptions, ascertainment, and calibration. This was coded with the help of Codex/GPT 6-Astra. I also apply it to my example data, for which it is shown that a model with asymmetric birth & death rates and an innovation rate fits much better than the currently allowed CAFE models with symmetric birth & death rates & no innovation rate parameter.

A Bayesian BDI model coded in cmdstanr would probably also be feasible. And aside from CAFE-type analyses I was thinking I could also use negative binomial phylogenetic mixed models to test if social species have more copies of a particular gene family than solitary ones (e.g. using brms), but not sure of the statistical power one would obtain.

Would you be interested in reviewing a smaller, opt-in contribution focused initially on the BDI transition kernel and root-state support? If so, what scope and validation would you want for an upstream contribution? Alternatively, if this is outside CAFE’s intended scope, I could also develop it as an entirely separate package. What do you reckon could be best?

[hahnlab-user]

Hi Tom,

Thanks for further explaining what you were interested in. I think we'd be very interested in extending CAFE to include this model.

There is no inherent limitation, other than a lack of signal with which to estimate the innovation parameter. My attempts more than a decade ago at this suggested that these parameter estimates were very unstable, likely due to a lack of signal. But of course genomes have proliferated since then, so it may be possible to do much more.

Perhaps we can discuss further over email? You can reach me at mwh@iu.edu.

cheers,
Matt