Scientist Stories: David Balding, How the causes of complex disease are distributed in the genome

Hello, I'm David Balding, and speaking to you from Melbourne, Australia, where we're in the middle of winter. Well, towards the end, fortunately. And I'd like to speak to you today about sort of one of the big picture questions in human genetics is where is all the heritability? Where are the causal effects in the human genome for complex traits? So we know the basics, that GWAS have really kind of revolutionized the field in the last decade and a half, and they've been getting better and better, and now we routinely have the order of a million individuals and 10 million SNPs. And statisticians will recognize immediately the problem that we have typically fewer observations than parameters requiring estimation. And so without further assumptions, we can't estimate the effect of every SNP. Now you might think, oh, we don't have a real problem there because we just test one SNP at a time. And of course, that was the original approach. And it does work to some extent, but the limitation soon became clear. You really need— because complex traits turned out to be very complex— you really need to look at the effect of each SNP in the context of the effects of all the other SNPs. So there's a lot of information in my talk, but I hope to just, you know, stick to the big picture, don't worry too much about details. I put some details in so you can look at slides later or our paper that I'll cite at the end. And the technology involved, the statistical technology, is not that great. It's linear regression, Gaussian distribution, not too much more. And so here it is, the workhorse of genome-wide statistical analysis is linear regression. There is an implicit assumption of additivity here, so that's a simplification. Seems to be OK. I'm not going to question that here. I mean, obviously there's quite a few things I have to skip over, and that's one of them, but most people assume additivity. And also most people assume, as I've indicated here, that the phenotypes and the genotypes are standardized. And by the way, for the purpose of this talk, I'm going to be talking about mainly about continuous phenotypes. Similar conclusions hold for binary phenotypes case control, but it's— there's a few extra complications that I won't go into today. So we've got standardized phenotype, standardized genotype, and under those conditions the regression coefficients when you square them are just the heritability of a SNP. And so heritability is a good currency for causal effects because it to some extent combines the effect size with MAF, and both are important, and we'll be discussing that more as we go along. So as I've said, we can't estimate all these betas. Even nice simple linear regression doesn't work because there's too many coefficients. Classic statistical problem, overfitting, so on. Now there's lots of ways to deal with this, and I just want to emphasize that there's no correct solution. They all— all the approaches require additional assumptions. None of the assumptions is true, But there is better and worse, and you can check that by fitting to data. So what I'm presenting today, I hope, is better than previous models, although— and I'll review the reasons. Now, because of this no perfect solution, the area has become quite a bit mired in controversy and confusion, and some people might think, oh, that's hopeless, we'll never get anywhere. I mean, I think we can, we are getting somewhere, and we can get a lot further. I hope for this audience I don't have to bother too much with why do we care. I mean, where causal effects are in the genome is crucial to really everything. The whole enterprise of the genomic era is based on this idea of understanding the genetic architecture of many traits, and things like effects that are common across traits, and also better prediction of traits. And so what I'm talking about here today is key to to all of those things, I think. So what are some of the solutions? Well, the early approaches were based on this sparsity idea that there were not too many SNPs involved in any trait. And unfortunately, we came to realize that nature isn't sparse. And so trying to select the right set of SNPs turned out to be not so fruitful, because the reality is it's not just a few loci controlling every trait. Everything's interconnected. There's many weak causal effects more or less everywhere in the genome for many traits. And crucially, we can't distinguish zero effect from weak effect. I mean, I've often been dismissive of any attempts to count the number of non-null SNPs or things like that, because you can't tell the difference between a null SNP with zero effect and a weak effect often. And weak effects on their own are not very interesting, but collectively they turn out to be important for many complex traits. So in recognition of that, nowadays we focus more on shrinkage, penalized, or Bayesian regression methods. So they're not necessarily fully Bayesian approaches, but they always involve a lot of Bayesian ideas, and I'm going to use the Bayesian language. The penalty function doesn't necessarily have to be a prior probability distribution, but really When you start thinking about the problem, there's no other way than, you know, in trying to decide what's the right penalty function to penalize large values of beta, the only way to do that is the classic way of basing a prior on our previous knowledge and understanding of the problem. The priors have to be symmetric because of the kind of arbitrary labeling of alleles. We don't know which one is going to be positive and which is the negative effect typically.. And in that case, the variance becomes the key parameter. And because I've said beta squared is the heritability, the variance of the effect size distribution is the expected heritability. And I'm going to be really focusing on that as an important parameter. So here's the idea. We deal— we regularize the problem by— with all of these, you know, million or so beta coefficients, we put a prior on them to reduce the flexibility in them, to reduce the overfitting problem. And we can show that in this way we can get good estimates providing— and we can get better estimates if we have better prior distributions. So that's a big focus of the talk. Now, we do have a lot of prior distributions. And this is a busy slide. And I don't want you at all to study all the details. This is just really to make the point that lots of people have come up with a ton of different priors. I want to emphasize this key one, which underlies the very popular GCTA software, but it really— when you start looking for it, it's everywhere in human genetics, a similar kind of assumption, which is that the effect sizes have a Gaussian distribution with some variance, but the key thing is there's no j subscript there because the variance is assumed to be the same for all SNPs. So this is basically ridge regression, but the sigma squared is not a tuning parameter. It is the key parameter of interest. Some priors use the t distribution instead of the Gaussian. We can have discrete mixtures of Gaussians or t's. We can have continuous mixtures. We can have a double exponential, as in lasso regression. I mentioned this hyperlasso, one of the early contributions from me and my colleagues, which didn't fly, but I think it had some good qualities. But yes, if you want to learn more, There's a nice review by Zhu and Zhao in 2020. This figure over here is just to sort of make the point that there's lots of different— this is to illustrate the idea of mixtures of Gaussians. So we have a whole lot of prior distributions, but the key thing I would sort of argue that people have been focusing on details of these priors and missing the wood for the trees, because the key thing is How does the prior differ from SNP to SNP? And the thinking so far has generally been, well, you know, we don't know. There might be different classes of SNPs with small effects and big effects, and we can impose a mixture distribution to describe the different classes of SNPs, but we don't know in advance which SNPs are in which class. So we have to compose these different components into a single prior distribution, which is the same for every SNP.. But we would argue that we do know a lot about SNPs prior to observing the GWAS phenotype data. We can get estimates of minor allele fraction or linkage disequilibrium parameters either from the genetic data, ignoring the phenotypes of the study at hand, or from reference databases. And we've got now a ton of information on functional annotations available from public databases. And our argument is that these should inform the prior in a way that has not generally been the case in the past. So I've already hinted at this focus on the variance of the prior as being a key parameter, and in fact, if you look at the comment on the bottom, we'd say that nothing else matters. And in fact, in our analysis, we don't really pay attention to much else other than the variance, which is the expected heritability of each SNP. And we've introduced the term a heritability model for something that assigns variances to the prior for every SNP. Now, if you're paying attention, you might notice that if you've got one variance for every SNP, you haven't solved this too many predictors problem. You've still got all these variance parameters that we don't know. And the answer is that although that's true, we can model the variances in terms of a lower degree of freedom model. So we can have a number of parameters describing the the variances, which is much less than the number of SNPs. So every SNP gets its own variance, but it's not a free parameter. It's going to be a model depending on properties of the SNP, like MAF and linkage disequilibrium and functional annotations. So as I've said, a huge amount of analyses in human genetics, just without any comment, implicitly assume a uniform heritability model, that this prior variance is the same for each SNP. So you just take the total SNP heritability and divide by the number of SNPs, and then we're going to say that we can do better. So how do we do that? Well, we set up a— sorry, I've got— I skipped a slide, something went wrong. So before I shift onto that, I just want to revise LD score regression, which was another landmark development in 2015. So you will appreciate that this was a key paper in the transition away from analyzing individual genetic towards summary statistic data, and there's a whole lot of reasons for that. And the fundamental idea of LD score was this idea that the more LD there is of a SNP, the more heritability it will tag, and that's represented by this linear regression here, and you can estimate parameters of that and get out meaningful heritability estimates and other properties. So that was a big advance, but we can develop that idea. But another point to emphasize is in these models there's no effect size anymore. We just— when we look at summary statistics, we don't even need an effect size parameter, but the variance parameter that I've emphasized is there, and you'll see now. So here's the LD score model. I might have slightly changed the notation from what you've seen before, but it's basically It's a generalized LD score regression model. And so there's the intercept term, which can be interpreted in terms of population structure, although we've shown in another paper that it doesn't actually do that very well. But I'm not going to focus on that here. Here's the LD score parameter. Now, it's a little bit different from what you've seen before because I put this variance parameter in for each SNP. The original LD score paper just had the uniform heritability model there. And so then that's terms— the LD score term simplifies quite a bit. But we're going to say that we can do better by having better heritability models there, and therefore an improved LD score. And again, just to reiterate that point, there's no effect size in sight. So whether the effect size distribution is t or Gaussian or mixtures is just— isn't even there in this analysis, but the variance of the effect size is there. And that's the expected heritability. So we're obviously facing a key question of how to choose a heritability model, and we started off thinking about this nearly 10 years ago and thinking about the effect of LD on sigma squared. Basically, higher LD— if you've got a whole bunch of SNPs in high LD, they behave a bit like a single SNP, and so we would expect them to have less heritability per SNP than if they were in low LD. And we proposed a model for that and showed that it worked better and so on. In retrospect, I think we would have been better to focus on MAF. So MAF and LD are obviously tightly connected. And arguably, at the time, we said after adjusting for LD, we had a look and said, oh, if you also adjust for MAF, you get a gain, but it's not so great. But actually, we didn't do the other way around, which is to do MAF first, which would have been better., and we would have seen that it is actually more important, but it's also easier to model, and also MAF is trait-specific in a way that sort of has a natural interpretation that I'll come to now. So to see that, remember that our effect size parameter was standardized. It's not strictly speaking statistical standardization, but it's very close according to a formula like this, and Now standardization has, together with a uniform heritability model, has a really strong implicit assumption about the relationship between MAF and effect size, and it's very strong and it's wrong, as I'm going to argue in a moment. So, you know, but because standardization is so prevalent, we're going to leave the standardization there, but adjust instead this variance parameter in terms of the MAF via this free parameter alpha. And alpha can be different for different traits, or you can use a kind of default value as we'll discuss in a moment. So alpha equals -1 is back to the uniform heritability model, and we can choose other alphas, and we discussed that first in our 2017 paper. So we've got a free range of ways to scale the effect size distribution,, and that has different implicit assumptions about the relationship between MAF and effect size. So 2 special cases: -1 gives you the uniform heritability model, which is very widely assumed in human genetics, often adopted uncritically because we're so used to standardizing variables we don't think about it. But because it implies constant expected heritability over SNPs, that turns out to say the unscaled effect sizes without any scaling at all increased rapidly as MAF decreases, and unrealistically, and this corresponds to very strong negative selection. Now, all the time that human geneticists have been adopting this bad assumption, a lot of plant and animal breeders have been adopting a different assumption, and there seemed to be little comment on the fact that these two different schools were doing very different assumptions about the MAF effect size relationship. So They assumed effect size constant over SNPs, which means the heritability decreases as the MAF decreases, and constant effect size, that's unscaled effect size, corresponds to assumption of no selection because it's selection that has differential effects on different MAF. So obviously then the question is now, well, which of these alpha values is better, or is there some different value? And we answered that most recently in our 2020 paper. I'll give you the references at the end. And here are some estimates for 14 UK Biobank traits. The red ones are the ones significantly different from zero according to some criterion. But remember, over here is the -1 that's assumed in nearly all human genetic analyses, and we've never seen any data that come anywhere near supporting that value. So it's a terrible assumption. The plant and animal breeders were on stronger ground if this human data can also be relied on in plants and animals. But we found, you know, surprisingly strong consistent value of -1/4 fits data pretty well for many traits. And because there's a fair bit of computation involved in estimating alpha separately, we've suggested -1/4 as a kind of default value, you'll notice that college education shows stronger, more negative alpha, which corresponds to stronger negative selection, which is a bit worrying for us who are positive about college education. But anyway, you know, there are possible explanations for this, and I won't go into that discussion here. Also, height selected against, which seems a little bit surprising, but I have to say, so, you know, the selection story is complicated and And I'm not really making any strong claims about the selection that's being imposed here. The striking thing is that overall there is evidence in the MAF effect size distribution for weak negative selection affecting most traits. But that's for the whole trait, a genome-wide average. What about particular bits of the genome? And this is an estimate of the value of alpha within certain regions of functional importance indicated on the side. We looked at more, and I'm only looking at the significant ones that are significantly different inside from outside the region here. And you'll see most of these regions, because they are functionally important, they show stronger effects of negative selection. So that kind of seems plausible. There's one exception here, this repressed category, but all the others show there is more negative selection within these functional regions than outside them, which is not too surprising. And also coding and coding-related things like promoter and 3' UTR have the strongest negative selection of the categories that we considered here. So having focused a little bit on genome regions like this, there was a bit of a breakthrough in 2015 when Finnucane et al., the Boston group, proposed the first— their first deviation away from the uniform heritability model, but in terms of these functional annotations. So they basically have some indicators A for whether or not a SNP is in a particular functional category, and then they have a parameter, the effect size of that category tau. And LD score is the— the original LD score is the special case of just one category,. But now they allow many functional categories. It doesn't really matter how many in a sense, because we've got on the order of millions of SNPs. So we ought to be able to estimate these tau precisely because we have a lot of observations. But in practice, the marginal precision is not that great because of strong correlation among the tau, and I'll mention that in a moment. And I also mentioned another step forward in 2017 from, again, from the the Boston group, that rather than just have binary indicators of whether you're in or out of a category, you can have a continuous variable describing, for example, LD properties of the SNP. And so they developed their model in that direction. So now we have lots of heritability models. There's the Finucane et al. one I just mentioned, has 53 parameters, a lot of functional categories. They also— each functional category has a buffer region of 500 base pairs either side, and a couple of other categories. Then the baseline LD model came along with these extra continuous variables in the model, MAF and LD, and so they bumped it up to 75 parameters. Some statisticians might be horrified that 75 is a lot of parameters to estimate, but as I said, in principle there's no problem. Now our contribution after that was to come along with an approximate likelihood and a likelihood based on summary statistics. I'll say a bit more about that on the next slide because I think it's important. And using our approximate likelihood, we showed that the baseline LD model was the best model available at the time. Not surprising because it had more parameters, and there was no problem really with overparameterization, so it's not surprising that those extra parameters did buy something worthwhile, and it was the best model. But with this insight, we could go ahead and improve on it straight away by getting rid of the 10 MAF strata. It's not a good way to model MAF. The continuous Sink-1 parameter model that I showed you earlier is much better. So we put that in, also changed a little the measures of LD, and reduced the number of degrees of freedom, and improved the performance of the model. And we proposed two versions of the model, whether you want to use this default value of -1/4 for the alpha parameter which specifies this MAF-effect size relationship, or you can have alpha as a free parameter to be estimated, which of course increases the degrees of freedom by 1. And if you're looking at many phenotypes, that increases the computational burden quite a bit. So likelihoods for GWAS summary statistics. I just want to mention that Zou and Stevens did develop a likelihood, and I've just given some brief details of it here, but they were focused on on doing Bayesian inference for heritability. And so they used MCMC. Their approach is coded in MATLAB. It's fairly computationally demanding, but it's also a little bit different in focus from what we had in mind. So they assumed you had parameters— sorry, the summary statistics were effect size estimates and their variances or standard deviations. We've assumed even less in terms of the data available. Heritability, just a single chi-square statistic. And then basically we model the distribution as a scaled chi-square. So now that the null distribution is the case of this expectation being 1, and then non-null SNPs have this expectation different from 1, and we can model that with our heritability model, which leads to our more general version of LDSC model. Now, so the marginal density assumed here is pretty good, and that's easy to check. The question is what to do about the correlations among the summary statistics, and we've done something pretty simple, dividing the marginal log likelihoods by a local sum of R-squared values to mimic the effects of correlation. We did some checks that that seemed to perform quite well, but I think that's an interesting area to see if we can improve improve on that and whether we can somehow merge our approach with Izu and Stephens' likelihood. But anyway, I won't say any more about that. Let's get back to the story. And what can you do with these improved heritability models? Well, here's a recent paper from Zhang and Speed in Nature Communications, just come out. And this line in the middle represents the predictive performance of a typical kind of PRS in predicting a number of phenotypes. They looked at both whole genome data, individual genome data prediction, and summary statistic prediction. This is the results for the individual genome, but you can look in their paper for summary statistic results as well. And then there are various improvements possible by regularization, because the classical PRS doesn't have any regularization. So just the simple GCTA model with the uniform heritability model already gives you a big improvement. And the LDAC-Thin model is a very simple simplified version of our model, but our best model is the BLD-LDAC model, and it gives quite a bit better prediction. And what I want to emphasize now is these different bars here correspond to different prediction models. So you can take any existing prediction model and improve on it by taking out its implicit assumption of uniform expected heritability and putting in a better heritability model. So we can improve everyone's prediction in this approach. Sorry, we— that's— I'm getting excited. It's Zhang and Speed. I'm not an author on this paper. And the They found that on average their improvement in prediction performance by just by putting in a better heritability model corresponded to about a 25% increase in sample size. So I think that's pretty impressive, and as I just emphasize again, we can improve pretty much anyone's prediction model. And so of these ones here, the BoltPredict and the BayesRpredict seem to have been the best,, but they get improved even further by our heritability model. So what else can we conclude from our improved heritability model? Here are some further updated estimates of heritability and this alpha parameter, which is interpreted as a measure of selection, and we're improved. This is including more data here and including more SNPs down to half a percent frequency, and so now that -1/4 estimate of alpha has changed a little bit. To -0.18 is our best effort with the low— our best estimate with the low MAF SNPs. But it's kind of interesting that still the alpha parameters are pretty much all consistent with this -0.18. Again, only college education is significantly different from that overall mean value, and these starred ones are significantly different from zero. The heritability estimates are for range from high for height down to some estimates that are quite small, but they're all significantly different from zero, the heritability estimates. And they're all, as shown here, approximately equally precise, so I didn't put in the SDs for the heritability estimates. And this kind of big question of where is heritability in the genome, we've estimated two parameters on the— you know, this is a fairly substantial analysis, with a lot of details that we can— that I'll have to skip over here. But the two measures are the relative concentration or enrichment of heritability, which just says how much extra heritability is there in this region than in a genome-wide baseline without caring about overlapping categories. But this tau is the parameter in our heritability models, and it's estimated in a multiple linear regression framework So each tau estimate does take into account the other tau values. So that's part of the reason why these lambdas are estimated much more precisely. The taus are estimated less precisely, and because they are correlated with each other because of overlapping categories. But now, and so there is room for model improvement here, and we've been working on that, but it's not completely straightforward because of these correlations I'm talking about. But what can we see? So conserved is not quite a functional category, it's conservation across mammalian species, but it shows that, you know, by both measures of the parameters, it has the biggest impact. But in a sense, that's not that interesting. It's not surprising that conserved regions have a high concentration of heritability. We have to drill down into other functional annotations to understand why, and of course, a lot of these conserved regions turned out to be coding, so there's some overlap there, but there's other things here. But after conserved coding has the highest concentration of heritability by both measures, and then 5' UTR, transcription sars-cov-2, 3' UTR. So none of those is too surprising. What maybe is a bit surprising is that many of these functional categories that we think are important are not that highly enriched for heritability. Maybe, you know, many of them are 2- to 3-fold and many even less promoters. So introns are slightly enhanced for heritability. Many of these tau values, as I said, are not marginally different from zero, but one interesting one is DHS. That's been much disputed in the past with estimates that are all over the place, but our latest estimate is there's the 70% enrichment of heritability in DNA's hypersensitivity sites, but actually tau is negative. So So being in a DNA's hypersensitive region actually decreases your heritability. The 1.7 being bigger than 1 is because there is higher heritability in those regions, but they are explained by other functional annotations within those regions, not just being in the DHS region. So we think there's some fascinating results here. We think there's a lot we can do— further we can do with this. It's going to be really important for understanding genetic architecture of traits. And I just need to get up my last slide, which for some reason I'm having trouble with. So I— there we are. There, I'll leave you with our recent papers, and a lot of the ideas in this talk are in a paper that's submitted. I'll put it up on an archive in a moment. We're just hoping to get referee reports and we'll do revisions and then we'll put it up, or you can email me for that, for my slides, or for the, for a preprint of the paper. So thank you for your attention.