Calculate change in frequency for nodes in tree #595
Conversation
|
858d2d0 updates the logic of these annotations to exclude tips by default and also omit annotations for zero-valued change in frequencies. The tree view now looks like this:
Even though the tips are gray, they still obscure the pattern in the internal nodes a bit. The scatterplot view is much more helpful for identifying recently growing/declining clades:
This most recent commit also adds an annotation for the current frequency of each internal node with a non-zero change in frequency. This annotation allows us to make plots like this one where we can identify rapidly growing clades that are still at lower frequencies, for example.
|
|
The build from the example screenshots above is available in the staging bucket and through the Nextstrain development server. |
|
Awesome! This PR looks like it does a straight-forward additive change in frequency? Ie going from 5% to 10% in 4 weeks would get labeled as dfreq = 0.65? This is not the best way to capture logistic growth. I would recommend a Fisher-Pry transform. More information can be found here:
But all you need to do is run:
on the input frequencies and then a simple regression will give you logistic growth rate. Most people are reporting these in terms of daily logistic growth rate, where for example B.1.1.7 is ~0.04 per day in the US and ~0.08 per day in Switzerland.
The above comes from working with lineage frequencies for a little while now and trying to assess their growth rates. Linear regression on Fisher-Pry transformed data has been the most robust method for me. I would:
Or we could give it a different name to avoid confusion with usage of Curious for @rneher's thoughts here as well. I do see these estimates breaking down when calculating the growth from 99% to 99.9%, etc... Perhaps we'd want to remove estimates for nodes with current frequency >99%, but I'd wait do to this until we can look at behavior. |
|
And separate from the issue of exactly how to calculate growth rate from vector of node (clade) frequencies, there is an issue of how to treat tips. I see that you've made the decision to not give attach a growth rate to tips. This is well justified in a lot of ways, however, I think it would be better to cast growth rates to tips. Take a look at this zoomed portion of B.1.1.7 for example:
Every internal branch gets an estimate of growth rate, but the final branch leading to each tip does not. These 2 tip internal branches will be about as wonky as single tip branches. This causes it to look like the internal bits of B.1.1.7 aren't growing, but they really are. Here, I'd set a tip count threshold (perhaps something like 10) in which only internal nodes with 10 or more subtending tips would have the logistic growth directly calculated. Other branches (with 9 or fewer subtending nodes) would inherit the closest calculated "upstream" growth rate estimate. This will emphasize single-tip branches coming from basal expansions, but I think this is okay. I'd also calculate current frequency for tips in addition to 2-plus-subtending-tip internal nodes. This will highlight the contribution of more recent tips. Both these will result in (I believe) more useful scatterplots where you can actually mouseover tips. |
Yes, this is just the original dfreq formulation from flu, so 5 to 10% in 4 weeks (0.0833 years) would be
Ok, I’ve added this logistic transform and regression to the delta frequency script as an option along with a user-configurable option to limit these calculations to nodes with at least 10 tips. The script propagates the values of parents to nodes that are too small to have values calculated. I reannotated the original tree that I’ve shown in the PR with a
One issue here is that the color scale is linear, so growth >3 looks the similar across the whole tree and it is difficult to identify clades at the higher end of the distribution. The coloring of tips also obscures the signal of internal clades even if those tips have the same color as their parents. Would it be difficult to implement a log color scale in Auspice or should we define our own scale in a different way through a config file? Next, we can look at the scatterplot of logistic growth by sample date.
This view is complicated by older tips that have inherited the logistic growth values of their parent clades. We also see there is a wide range of growth values with the most extreme values representing small clades (those with <10 tips). This pattern seems consistent with the logit transformation’s greater emphasis on small changes at the tails of the input distribution. Just for comparison, I repeated these growth annotations but without propagating values from parents to children and got the following view that more clearly highlights rapidly growing smaller clades like the one at the top. I understand the need to include tip annotations in the tree view, though. For this kind of analysis, I wish that the scatterplot view had a way to toggle tips on an off the way it can toggle branches.
Back to the original tree with tip annotations, if we plot growth by current frequency, as shown below, we still see the wide range of growth values on the left, but we can also more clearly identify medium-sized clades (by current frequency) that have larger growth values (like 501Y.V1).
We can use this view to zoom into a recent 501Y.V1 subclade with higher growth. From the tree and map views, we can see that this growing clade has a declining subclade with a ORF1b:P218L mutation in West Africa and a growing subclade with a N:R204P mutation in Australia and central Europe.
One issue that you hinted at before is that propagating growth values to smaller clades can give the appearance that these small clades are doing as well as their parents when we really don’t know. This pattern is easier to see with the linear dfreq measurement and linear color scale for the following clade where its very small subclades appear to have the highest change in frequency, but they’ve really inherited their parent’s value.
Finally, here is the correlation between the logistic growth values and linear change in frequency values. For early detection of small clades that are rapidly growing, the logistic growth values provide a much better signal.
If this looks like a better direction, I'll update all the Auspice config files to reference these new logistic growth and current frequency fields. |
|
Thanks a lot for pushing this, @huddlej ! The large negative values after log-transform could be driven by the KDE smoothing kernel and will depend rather sensitively on the value of The behavior of frequency trajectories near one are likely never a perfect logistic. Geographic and demographic penetration will slow a variant down before reaching 100%. But logit should still be a better transform than the simple log that we have used in the past. log transform has the advantage though that variants will receive smaller growth values once they become common (growth from 90 to 95% is arguably not terribly important). Overall, frequencies are normalized such that they sum to one at each pivot, right? I have been thinking of using the explicit derivative of your smoothing kernel, but this is complicated by unequal weighing pivot by pivot. |
That makes sense. For the figures above, I used the default value from the diffusion frequencies code (
With
That's right. We could always switch to a more flexible implementation though. |
|
Thanks for all the detailed work here @huddlej. Starting with a replication check using different methodology to check that final numbers make sense. Here, I looked at the base of B.1.1.7 which has a logistic growth rate of 4.0 in your analysis (https://dev.nextstrain.org/staging/ncov/global-with-dfreq?c=logistic_growth&f_clade_membership=20I/501Y.V1&l=scatter). It looks like this is measured in If I take the same 3984 strains present in your dataset and run my logistic rate estimate I get:
This rate of 0.04 per day corresponds to a rate of 13.0 per year, so about 3X what you have. (This is a little funny of comparison in that the data goes until Feb 15 in this example, but the final pivot is March 22, so the final 5 weeks in the PR example is almost all projection from the KDE.) I followed up here by directly looking at the final 4 pivots of clade 501Y.V1 via Auspice which are 27%, 29%, 30%, 32%. Logit transform followed by linear regression gives a slope and logistic growth rate of 1.1 per year. I think I need to run things myself to dig in deeper here to figure out what's going on with the difference between values of 1.1 and 4.0, but do you have any thoughts on what may be causing this John? |
|
Also, I think we're going to want 6 pivots (weeks) back. If you look at the latest live dataset, the final 2 weeks have ~30 genomes and the 2 weeks before this have ~330. I'm comfortable making estimates using 4 weeks of data, but the most recent two weeks shouldn't count for much here. We see consistently that the nowcast of B.1.1.7 frequency based on KDEs is less than "real" frequency once the data comes in. For example: https://nextstrain.org/ncov/global/2021-03-16 has pivot from 2021-02-28 at 36%, while https://nextstrain.org/ncov/global/2021-03-01 has this same pivot at 25%. Either in addition to or instead of this, we could also pull back the latest pivot we calculate in general. This may be more principled actually. As we don't trust the current estimate. Maybe lop off the final pivot in current tip frequencies and move to 5 pivots back for growth estimates? |
I think this is a bug with how I ran the build on March 22 but did not download the latest GISAID data. Hopefully, this isn't going to be an issue in production, but I'll doublecheck.
It seems like the issue might be differences in how we calculate our date units. One small point is that the delta frequency script asks for "4 pivots back from the last pivot" to get the frequency 4 weeks back, so we end up with 5 data points to fit. I attached a Jupyter notebook showing the frequencies and pivots for 501Y.V1 and how the model fitting works for those data: In your example of 32% latest frequency and earliest of 27% frequency, the delta time between those pivots is
That is a great idea. With this change, then the |
|
Thanks for the well thought out notebook @huddlej. This really makes things clear. Seeing how close to perfect fit the linear regression is on the logit scale is really comforting. There's perhaps a bit of evidence for the final pivot to be a slightly below expectation, but it's really not as big of an effect as I was expecting:
I was proposing an "upstream" change so that when That said, the current B.1.1.7 example looks pretty spot-on without the final pivot being obviously off. I'm not sure I'd worry about this at the moment. Given notebook, I'm quite convinced about the logistic growth calculations and think the current parameters (4 pivots back, so 5 pivots total) are probably good.
Ah! Yes. I messed up with the delta time and was off by 1 week vs 4 weeks to give the factor of 4 error. Thanks for catching this. |
|
Also, thanks for all the work implementing the tip count cutoff for calculation of growth rate. I was a bit surprised by this behavior here:
where the highest growth rate clade does not have anything special about it. It's just a clade of 10 tips with no nucleotide mutations and nothing particularly special about it.
It looks like it just happens to have 10 recent tips associated with it that are sampled and so there is apparent clade growth:
Maybe the only recourse here is to have a larger cutoff value? Like 15 instead of 10? Because even if we didn't propagate internal node growth rates to the tips, this would still register as a internal node with 10 tips with high growth rate. Or perhaps dropping the final pivot would help in this scenario? I can imagine these small clades being influenced by a particular very recent sample pulling up the final pivot. |
Ah, that seems much nicer!
I bet dropping the final pivot would help. Initially, I tried setting a minimum current frequency of 1% instead of a fixed number of tips (after talking with @kistlerk about her similar work). But it turned out that viewing these data with current frequency on the x-axis helped focus on clades that are large enough to consider their high logistic growth as more relevant without applying any filters during the analysis. I will re-run the analysis with the latest available data so the last pivot and available samples are more closely matched. We can see from those results whether we need to modify any filter thresholds or alter the pivots. |
|
Thanks John. I had run this trail build https://github.com/nextstrain/ncov/actions/runs/717668977 from commit 53053e2. You should be able to work from this: I haven't looked at the results here yet. |
|
@trvrb I ran a fresh build with the latest available GISAID data, to see if logistic growth annotations were less likely to pick up spurious signals from KDE distributions. As a sanity check, I inspected the frequency distribution of strains by their submission date. Although there were strains in the tree from early April (a couple of days ago), 48% of frequencies at the last timepoint came from strains submitted >7 days ago and 21% submitted more than 1 month ago.
This distribution suggests that calculating frequencies from one week back as the maximum date would be a better implementation. Next, I looked at the logistic growth values by sampled date and current frequency. The same issue that we noticed in the earlier analysis appeared here where the clade with the highest logistic growth was a small one (N=10 strains) with two more recent samples whose KDE distributions most likely drove the apparent recent growth.
To address this issue, I added a configuration parameter for frequencies called Without defining this parameter, its value defaults to 0 and the max date is (As an aside, while I was inspecting the frequencies command, I noticed that the KDE bandwidth is 0.05. At 0.00274 years per day, this bandwidth corresponds to 18 days. That seems kind of wide now that we have such densely sampled trees.) Other options to tune parameters and avoid these spurious signals would be to:
I reverted the frequencies censoring back to original behavior and then I tried increasing the numerical precision threshold to 1e-2, as a way to avoid censoring data while limiting extreme logistic growth values. This change set an upper threshold of ~12 for the upper end and ~-7 for the lower end. The resulting tree view more clearly highlights recently growing clades.
The highest growth clades with current frequencies >1% are actually reasonable clades with many amino acid mutations.
Some larger clades with high growth rates did not have amino acid mutations.
Their growth appeared to be driven by smaller subclades with a couple of recent samples.
This result suggests that censoring the last couple of days in frequency calculations could still be helpful. Adding censoring of the last 7 days back to the frequency calculations, the same larger clade still has a relatively high growth rate, but the subclade that appeared to be driving that growth is no longer solely responsible now that its two latest strains do not contribute to frequencies.
For what it is worth, this approach of censoring the last 7 days for frequency calculations and increasing the numerical precision threshold seems to produce the clearest signal from the tree view, too.
I’ve pushed this latest tree to the staging server, if you want to explore it more. |
|
Nice progress here John. However, I do think that the shrinkage from using
The base of the clade is estimated to be growing as expected from just looking at frequency plot. However, subclades within it quickly dive to 0. My desired behavior would be a tree that's much more blocky. Ie we expect subclades within 20H/501Y.V2 to be growing at the same rate as the parent clade unless we have good evidence otherwise. I'd suggest to go in the direction of something like |
|
Ok, I see what you mean, @trvrb. Here is the same view for 501Y.V2 using
I added the frequency threshold ( The "logistic growth" and "current frequency" fields now have proper human-readable names now, too. |
|
Thanks John! This is behaving rather nicely I think. I'm quite happy with the "blockiness" of the tree where moderately sized clades are being estimated to have similar logistic growth rates. A couple observations: The clade with the highest logistic growth rate is a subclade of 20I/501Y.V1. And generally as expected subclades of 20I/501Y.V1 rank above much of other circulating viruses. You can see the growth of 501Y viruses here: https://nextstrain.org/staging/ncov/global-with-dfreq?c=gt-S_501&d=tree,frequencies&l=scatter&p=full&scatterY=logistic_growth, which is really quite striking.
I'm mostly okay keeping "Change in frequency" in, but it does have some odd behavior and I prefer how logistic growth is working. For example take a look at this parent clade of 20H/501Y.V2.
You can see that this parent clade (part of 20C in green) is estimated at higher "Change in frequency" than the 20H subclade (in orange). However, I think this just because this 20C parent clade is at higher frequency and getting carried along by a rising 20H subclade. If we look at logistic growth it nests in as expected and is mostly lower than 20H.
This parent clade to 20H is actually the thing with the highest "Change in frequency" and rather stands out in the scatterplot above all the 501Y viruses. (This is the main thing that gives me pause about including this, people will see clades like this one and think something is going on when it isn't) My only remaining complaint is color ramp, where for logistic growth the vast majority of nodes have values greater than -10, but because the minimum value is around -30 the color ramp is compressed to basically be oranges and reds. I think we should supply a min and max in the Auspice config JSON to ramp from -10 to 10. The other way to tackle this would be to increase This deserves a bit of testing, which I'll try to do. |
|
As specified in the above commit I thought that linear delta frequency was just too confusing. For example, parent clades of B.1.1.7 were consistently getting assigned higher values of delta frequency. Generally, I think is too biased to assigning high values to high frequency clades, which is exactly what we don't want to target. I spent a fair bit playing around with different parameters for logistic growth calculation. There were a few broad considerations:
I'm pretty happy with the results. By my eye the trees generally have the appropriate behavior with known clades and also with clicking around to look directly at frequency behavior. Results are available at:
|
|
I'm sure this could use some further finessing, but I'm happy enough with things as they stand that I'd now be okay merging. Thanks for all the detailed work on this PR John. |
|
I think this generally works very nicely. But it does generate some odd dead branches: I am wondering whether a more robust approach would be to consider only clades with tips in the last 3 months or so. |
|
Thanks Richard. Your example here is the edge case I was having the most trouble dealing with. It's been difficult to get the small basal clades that don't hit the 50 tip cutoff to register appropriately. If we want the basal 20I/501Y.V1 clades to be appropriately red, it's really difficult not also have the basal clades to the parents of 20I also be inappropriately red. The actual parent clade to 20I/501Y.V1 should have high logistic growth, but small offshoots to this parent lineage should not. That said, a cutoff for only calculating logistic growth rates for clades with tips sampled in the last ~3 months seems like it could be a generally helpful thing. Worth experimenting with. Also, I see that the merged PR is failing the CI tests. Should look into this as well. |
Adds a script, Snakemake rule, and configs to calculate the change in frequency (i.e., "delta frequency") for nodes in a given tree over a fixed period of time. The default time period is 4 weeks which corresponds to 4 timepoints or "pivots" in the tip frequencies data. The workflow exports this new `delta_frequency` annotation in auspice JSONs for visualization as a color-by or scatterplot attribute. Fixes #594
Do not annotate change in frequency for tips by default (but allow users to override this with a flag) and do not emit entries for nodes with zero-valued delta frequencies. The resulting new annotations only color internal nodes and branches with non-zero delta frequency values which allows us to more rapidly identify clades of interest.
Adds a `--method` argument to the delta frequency calculation script and logic to calculate a regression from the logistic transform of frequencies. The resulting annotation estimates logitic growth for clades. This new annotation is included as "logistic_growth" in the auspice JSON along with the original linear annotation, for comparison.
Limits delta frequency annotations to nodes in the tree that are currently at >=1% frequency and propagates the delta frequency (or logistic growth) values of these internal nodes to their lower-frequency children. Also, separates the logic for annotating the current frequency of each node such that tips have their own frequency annotations. These changes should help us identify clades of interest from both tree and scatterplot views.
…calculations Adds a configuration parameter, `recent_days_to_censor`, to the top-level `frequencies` configuration group. When this value is defined, the max date used for frequencies is today's date minus the number of days requested. The value defaults to 0, for backward compatibility.
Increase the threshold for values in the logit calculations used for logistic
growth. This change has the effect of restricting the range of logistic growth
value such that small clades with a few recent samples do not appear to dominate
the signal of rapidly growing clades.
args = parser.parse_args()
@@ -116,7 +115,7 @@ if __name__ == "__main__":
# transform (as when frequencies equal 0 or 1).
y_frequencies = logit_transform(
node.frequencies[first_pivot_index:],
- pc=1e-4
+ pc=1e-2
)
# Fit linear regression to pivots and frequencies and use the
Adds config and Snakemake param for the frequency threshold parameter of logistic growth calculations. Also, adds a separate logistic growth section to the config, so these calculations can be parameterized differently from linear delta frequency calculations. Sets minimum tips to 25, to make sure annotations only occur for clades that are large enough to warrant inspection.
Parent clades of B.1.1.7 were consistently getting assigned higher values of delta frequency. Generally, I think is too biased to assigning high values to high frequency clades, which is exactly what we don't want to target.
1. Increase min tips requirement 2. Increase delta pivots 3. Include max frequency cutoff 4. Include min frequency cutoff
…te logistic growth The latest CI builds have fewer samples than previous CI builds, so the new requirement that nodes have >=50 tips to calculate logistic growth is never true for the CI data. This leads to a bug in the logistic growth script where the root node tries to inherit its logistic growth value from a parent it doesn't have. This commit adds a specific check for the existence of a parent node before trying to pull from that node and assigns a default missing value in the edge case where the parent does not exist.
Instead of hardcoding the minimum frequency for clades to get logistic growth calculated, expose this value as a parameter on the command line and in the workflow configuration. This change should allow us to tune parameters more easily.
Reduces the stringency of default logistic growth settings such that values can be calculated for the smaller tree in the CI build and we make sure the growth script is tested properly.
huddlej
force-pushed
the
add-delta-frequency
branch
from
cb41320
to
41e60ea
Compare
@trvrb, I rebased this branch on to master so we're up-to-date with the new CI dataset on master. 2f875d6 fixes the CI bug and 41e60ea enables logistic growth for the CI build by making the internal node requirements less strict.
Do you want to continue experimenting in this branch or merge and then continue experimenting separately? It seems like a higher minimum current frequency threshold (instead of a high minimum number of tips) would provide a way to get at the "clades with tips in the last 3 months" approach. 2130500 adds the minimum frequency as config parameter that we could tune for this purpose. |
Only calculate logistic growth for clades with frequency greater than a threshold (set to e-06). This effectively only calculates logistic growth for clades sampled in the previous 3 months.
|
This did reveal one idiosyncrasy to how
where very recent tips are gray. These are generally viruses sampled in the past week. I looked into things and these viruses (like "Australia/NSW4424/2021" in global above) have a collection date after the with all zeros. This is not what I would have expected. I'd expect I don't think this makes a big difference however. But was surprising to me and slightly shifts clade frequencies. |
|
(I just realized I was supposed to update the change log as part of PR with this as a feature. I'll edit |
|
Love this new feature, @trvrb and @huddlej . For use cases with time periods where the absolute number of samples is low and/or the variant is at very high or very low frequency, a weighted-least squares regression in logit-space may be more appropriate. (The strategies you have taken in terms of thresholding may also address these issues, but it may be worth trying for your offline analyses @trvrb.) The loss function would be Sum (f(x_i) - y_i)^2/(var(y_i|x_i)^2 where x = week, y = logit(frequency of variant in week), and var(y|x) is due to finite sample effects, rescaled by the squared derivative of logit. This is implemented in statsmodels and sklearn This can help with time windows lacking good data and also with low/high probabilities, while still being as fast as linear regression. Back-of-the-envelope: Finite sample correction: if you observe k out of n, then the variance of k/n is
(If n is very low, you may want to use a beta(1/2,1/2) prior, which amounts to adding 0.5 to k and 1 to n.) Jacobian of logit the logit will rescale a frequency x to log(x) - log(1-x), with derivative 1/x + 1/(1-x). Plugging in the observation k/n to this formula gives n/k + n/(n-k). Put this all together to get a variance in the logit observation of approximately var(y|x) = k(n-k)/n * (1/k + 1/(n-k))^2 When k is small and n is large, this is about 1/k; when k = n/2, this is about 4/n. The appropriate weights for the regression would be the reciprocals of those. So an error at small k would be weighted by about k, where an error for large k would be weighted by about n/4. |
Description of proposed changes
Adds a script, Snakemake rule, and configs to calculate the change in
frequency (i.e., "delta frequency") for nodes in a given tree over a
fixed period of time. The default time period is 4 weeks which
corresponds to 4 timepoints or "pivots" in the tip frequencies data. The
workflow exports this new
delta_frequencyannotation in auspice JSONsfor visualization as a color-by or scatterplot attribute.
Related issue(s)
Fixes #594
Testing
Tested by CI and with the Nextstrain "global" build. Examples include scatterplot view:
And the tree view:
Based on the tree view, we may want to limit the delta frequency calculations to internal nodes instead of calculating them for all tips.