2.1 Measures of effect sizes

At the heart of a PRS is the application of genetic effect sizes estimated from one study, called here the “training” study, and then used to predict phenotype in another study, called here the “test” study. In order for this procedure to make sense at all, the investigator must be explicitly assuming that the estimated effect in the training study is similar or the same as the effect in individuals in the test study. Thus, to understand PRS we must first examine what sort of genetic effects might reasonably be assumed to be the same. To be specific, begin by considering a single locus in Hardy-Weinberg equilibrium with two alleles A₀ and A₁ [23]. Let p be the frequency of the A₀ allele and q = 1 − p be the frequency of the A₁ allele, and assume that we have labeled the alleles such that p ≥ q.

2.1.2 Binary traits

Human disease studies are often most interested in binary phenotypes, diseased or not diseased. For binary traits generally the only reported effect size is an odds ratio, OR (defined below). There are very practical reasons this is the case [23]. However, to use our quantitative genetics tools, we generally [25] model a binary phenotype as resulting from the existence of a threshold t on an unobserved normally distributed phenotype, L, that we usually call “liability” to the disease in question. Here we make a stronger assumption than needed for most continuous trait analysis: we assume liability is normally distributed, parameterized to have mean 0 and variance 1 for computational convenience. Thus, liability is assumed to follow a standard normal density ϕ(x), with standard normal cumulative distribution Φ(x). The overall population prevalence of the disease, ψ, is uniquely determined by the threshold t (Figure 1), and vice versa.

Figure 1 Normally distributed liability with disease-determining threshold at liability greater than 2.

(18)

ψ=∫t∞ϕ(x)dx.

(19)

t=Φ−1(1−ψ)

The penetrances, ζ₀ and ζ₁, of alleles A₀ and A₁ with additive effects α₀ and α₁ are defined as the probability an individual is diseased given they have the allele in question. Using the approximation [23] that 1 − V_a ≈ 1, as it is for virtually all known alleles contributing to complex disease phenotypes [26],

(20)

ψ=pζ0+qζ1.

(21)

ζ0≈∫t−α0∞ϕ(x)dx,

(22)

=1−Φ(t−α0).

(23)

α0≈t−Φ−1(1−ζ0).

(24)

ζ1=ψ−pζ0q.

(25)

α1≈t−Φ−1(1−ζ1).

In this we see that the additive effect of an allele uniquely determines its penetrance, and vice versa. The odds ratio (OR) of allele A₁ to A₀ is defined as

(26)

OR=ζ1(1−ζ0)(1−ζ1)ζ0

(27)

=ζ1(1−ψ−qζ11−q)(1−ζ1)ψ−qζ11−q

(28)

=ζ1(1−q−ψ+qζ1)(1−ζ1)(ψ−qζ1)

(29)

=ζ11−ζ1⋅p−(ψ−qζ1)ψ−qζ1

Presented in this format, it seems far more intuitive that an OR is fundamentally a function of three things: the frequency of the alleles, p, q, the penetrance of alleles, ζ₁, and the overall prevalence of the disease, ψ. That it is possible to estimate the OR without explicit knowledge of the prevalence or allele frequency does not change the fact the underlying quantity is per se a function of frequency and prevalence. If two populations have the same odds ratio, OR = OR*, then

(30)

ζ11−ζ1⋅p−(ψ−qζ1)ψ−qζ1=ζ1*1−ζ1*⋅p*−(ψ*−q*ζ1*)ψ*−q*ζ1*.

(31)

ζ1(1−ζ1*)ζ1*(1−ζ1)=(p*−(ψ*−q*ζ1*))(ψ−qζ1)(p−(ψ−qζ1))(ψ*−q*ζ1*).

Thus, if the odds ratio of allele A₁ to A₀ is the same in two populations, then the odds of allele A₁, ζ11−ζ1, in one population divided by the odds of A₁ in the other population, ζ1*1−ζ1*, is a ratio of the allele frequencies in the two population multiplied by a ratio involving the prevalence, frequency and penetrance of the allele. If we further believe the “reason” the two populations have the same OR is because the penetrance of A₁, the chance of developing disease given you have an A₁, is same in the two populations then

(32)

ζ1=ζ1*=ζ.

(33)

(p*−(ψ*−q*ζ))(ψ−qζ)=(p−(ψ−qζ))(ψ*−q*ζ)

(34)

ψ−qζp−(ψ−qζ)=ψ*−q*ζp*−(ψ*−q*ζ).

Two populations with same OR can have the same penetrance for an allele if and only if they have the same allele frequency in both populations, p = p*, and the same prevalence in both populations, ψ = ψ*, or the ratio of allele frequencies between the populations is somehow “forced” to be determined in a rather complex way by the penetrance, prevalence and frequency in each population. While we can imagine it is possible for this complex relationship to exist at some particular moment in time, if the allele frequency in either population were to fluctuate slightly, the prevalence and penetrance would have to simultaneously alter in a precise manner to track the allele frequency perturbation. Effectively this can only be true for a disease model where allele frequency, allele penetrance, and population prevalence are all tightly linked together, effectively a Mendelian disorder where the frequency of a variant determines its penetrance and population prevalence. For any complex human disease with many genetic and environmental factors each contributing only a tiny fraction of the total variance, the precise relationship required between the two populations’ allele frequencies, disease frequencies and allelic penetrances will never hold. Put more intuitively, the model where the OR can be the same between two populations naturally arises out of a model where the frequency of a disease is determined by the frequency of the allele. For a complex human disorder where no allele has any strong effect on disease, the frequency of disease is more realistically thought of as nearly independent of the frequency of any one allele, and as such it is impossible for the OR of a complex trait to be the same between two populations with differing allele frequency.

2.2 PRS implementations

As a general rule, PRS analysis begins with the estimation of the effect sizes, generally reported as a β for a continuous trait or as an OR for a binary trait, at a large number of SNPs in a large collection of individuals with some phenotype of interest in the training study. After determination in the training study, effect sizes are generally taken to a different study population, the test study.

Assuming the training study estimated effects at a very large collection of genomewide SNPs, generally the first step in forming a PRS is to “thin” the markers. Thinning markers is a very practical approach to the challenges induced by LD. As discussed in some detail in the first in this series, LD, correlation between allelic states at neighboring SNPs, causes departures from additivity between those neighbors. The combined genetic effect of two markers in LD is, as a rule, very different than the sum of their individual effects, and in many circumstances smaller [23]. In the Kempthorne framework, we see that LD often leads to a “negative” interaction between neighbors, and the joint variance explained by two markers is smaller than the sum of their individual variances. From a Falconer viewpoint, one would say the estimated effect at one SNP is inflated by genetic effects of its neighbors. With this world view, one might intuitively think about a given region where there is only one SNP with a “real” effect, but all of its neighbors who themselves have no “real” effect will have inflated estimated effects caused by LD with the one “real” SNP.

As we previously suggested, there are straightforward approaches to model LD and estimate what effect sizes would be in the absence of LD (i.e., the “true” effects with a Falconer view), but most PRS applications take a slightly different tactic. In the most straightforward approach, often imagining that within any small genomic region only a single SNP will have a large LD-independent effect on phenotype, the SNP with the largest effect in the training study is selected first. Next all other SNPs with sufficiently large LD with the picked SNP are eliminated. This process then repeats, picking the remaining SNP with the largest effect size, and then eliminating all others in high LD with it, until all SNPs have either been picked or eliminated. The collection of picked SNPs, the SNPs after thinning, and their estimated effects will be the set of SNPs used to construct the PRS. There can be innumerable subtle variations on the SNP thinning algorithm [27, 28, 29, 30, 31, 32], including the application of some quite sophisticated machine learning algorithms, but the essential logic for all of them will likely be similar: find SNPs thought to have large effects individually or in combination, remove neighbors whose apparent effects appear to be explained by previously chosen variants, repeat. These processes might use explicit estimates of LD in the process or any sort of surrogate involving physical distances or other measures of correlation in allelic state.

2.2.3 Numerical toy example

To give the reader an intuitive understanding of what the above mathematics is describing, here we show a very simple-minded “toy” example demonstrating how changing allele frequency alone can dramatically affect a PRS. The situation described is not meant to be realistic, but instead to illustrate what could happen if there were systematic correlation between minor allele frequency and the direction of effect, combined with systematic differences in allele frequency between training and test studies.

Imagine what might be a “typical” successful PRS, one with 100 different SNPs, all unlinked to one another which in total explain five percent, V_PRS/V_P = 0.05, of the total phenotypic variance in the training study. For numerical simplicity assume that V_P = 1 and that the frequency of the common allele is the same at all loci with p = 0.8, all loci contribute equally to the trait, and so that the difference in allelic effects is the same at all loci with β ≈ 0.03953. Notice that in this toy example the minor allele is associated with increasing trait (liabilty) value at all sites, i.e., the sign of β is positive at all sites. If we were to plot the distribution of PRS across individuals in the training study we would see an approximately normal distribution with standard deviation ≈ 0.22 (Figure 2).

Figure 2 PRS distribution in the training study (Blue) and test study (Green). In this numeric example the only difference between the test and training studies is minor allele frequency which is systematically lower in the test population.

Now suppose we were to apply this PRS in a test study where the difference in allelic effects, β are the same at all 100 SNPs as they are in the training study. However, imagine that by some extremely unfortunate chance all the allele frequencies are different in a systematic way. Imagine that in the second population p = 0.9 at all loci. Thus, in these two populations the genetic effects are exactly the same, and the only difference is in the frequency of the alleles. As is likely to be intuitive to careful readers, the variance due to the PRS, V_PRS, in the test population will be less than the training study, because 2pq is smaller at each locus by a factor of 916, and therefore in the training population the overall variance explained by the PRS will be smaller by a factor of 916. The standard deviation of the PRS distribution in the test population will be 75% as large as the training population. The means, however, will be vastly different. The mean PRS in the test study will be ≈ 7 standard deviations below the training study (Figure 2). For most practical purposes these distributions do not overlap. In this numeric example, systematic differences in allele frequency lead to potentially detectable differences in PRS variance, but also lead to profound differences in mean. The intuition for why this happens derives from individual SNPs each having only a moderate effect on the PRS (β < 1), resulting in β > β². Since the contribution to variance for each SNP is 2pqβ², but the contribution to mean is 2qβ, the effect on mean is vastly larger than the effect on variance. This is not meant to provide a realistic picture of any real PRS data set but instead to establish an intuition for what could happen if there are two specific types of correlation not generally anticipated by most PRS investigators. If the sign of β is correlated with minor allele frequency, and allele frequencies vary in a systematic way between training and test studies, then PRS will vary between training and test studies in predictable ways. That there is likely to be correlation between the sign of β and minor allele frequency within studies is the topic of Section 2.6. That there is likely to be correlation in minor allele frequency between studies is discussed in Section 2.7 and beyond.

2.3 Using a PRS

PRSs are seen primarily in two related contexts. In the first context, we might wish to ask to what extent does a PRS developed in one collection of individuals (a training set) accurately predict the known phenotype of individuals in a second study (a test set), where we compare the predicted phenotype to the known phenotype in the test set. In this context we are fundamentally assessing what fraction of the phenotypic variance is explained by the PRS in this new collection of individuals. We could be doing this to compare and contrast differing methodologies used to construct PRSs [28, 34, 35, 36], or we might be doing so to establish the extent to which the phenotype of individuals in the second study appear to have a similar genetic basis as the first [37, 38, 39, 40]. In a different context, we might wish to apply PRSs developed in the training set to predict the unknown phenotype of individuals in the test set, given their known genotypes. It is in this context that a PRS might have significant medical utility [41, 42, 43], if, for instance, individuals with a substantial chance of developing a disease later in life can be identified long before such a condition develops, allowing interventions to be attempted earlier to reduce that risk.

2.3.1 Continuous traits

Developers of PRS methodologies frequently wish to show how well PRS predicts phenotype in a test study with known phenotype. For a quantitative phenotype the usual method of evaluation is to calculate the squared correlation between the PRS and the known phenotype, with the notion being increasing squared correlation implies improving PRS. For a perfectly trained PRS, the squared correlation between the PRS and phenotype should converge to the heritability explained by the markers used to construct the PRS whenever there is no dominance within a locus, there are no interactions between genes, nor between genes and the environment, and the state of all genotypes is independent of one another (i.e., no LD between sites contributing to the phenotype). To see this, imagine a PRS where the estimated β_v was equal to its true value at all included sites v, and imagine a test study drawn from a population with the same β_v at all sites. Letting PRS_j and P_j be the PRS and true phenotype of individual j in the test study, and letting SNPs n + 1, …, N be sites contributing to phenotype not included in the PRS, with M environmental factors also contributing, we find

(54)

Pj=(∑v=1nβvSvj)+(∑v=1N2αv0)

(55)

+(∑k=n+1NβkSkj)+(∑m=1Mem).

(56)

PRSj=∑v=1nβvSvj.

(57)

E[Pj]=0.

(58)

E[PRSj]=∑v=1n2βvqv=−2∑v=1nαv0

(59)

E[Pj−PRSj]=E[Pj]−E[PRSj]

(60)

=−E[PRSj]

(61)

Var[Pj]=VP.

(62)

Var[PRSj]=E[PRSj2]−E[PRSj]2

(63)

=E[ (∑v=1nβvSvj)2 ]−(∑v=1n2βvqv)2

(64)

=∑v=1nβv2E[Svj2]+2∑v=1n∑w=v+1nβvβwE[SvjSwj]

(65)

−4∑v=1nβv2qv2−8∑v=1n∑w=v+1nβvβwqvqw

(66)

=∑v=1nβv22qv(1+qv)+8∑v=1n∑w=v+1nβvβwqvqw

(67)

−4∑v=1nβv2qv2−8∑v=1n∑w=v+1nβvβwqvqw

(68)

=∑v=1nβv2(2qv−2qv2)

(69)

=∑v=1n2pvqvβv2

(70)

=∑v=1nVav,

where we use the lack of correlation between genotypes (no LD) extensively, i.e., E[SvjSwj]=E[Svj]E[Swj]. With similar analysis we find

(71)

Cov[PRSj,Pj]=E[PRSjPj]−E[PRSj]E[Pj]

(72)

=E[PRSjPj]

(73)

=E [ PRSj( PRSj+(∑v=1N2αv0)

(74)

+ (∑k=n+1NβkSkj)+(∑m=1Mem) ) ]

(75)

=E[PRSj2]+E [ PRSj( (∑v=1N2αv0)

(76)

+ (∑k=n+1NβkSkj)+(∑m=1Mem) ) ]

(77)

=E[PRSj2]+E[PRSj]E [ (∑v=1N2αv0)

(78)

+(∑k=n+1NβkSkj)+(∑m=1Mem) ]

(79)

=E[PRSj2]+E[PRSj](−E[PRSj])

(80)

=E[PRSj2]−(E[PRSj])2

(81)

=Var[PRSj],

where this adds the assumption of no interaction between any genetic and environmental factors. It should be further pointed out that the allele frequencies and additive effect of the A₀ alleles, pv,qv,αv0, above are the values from the test population. Nothing about this changes if the allele frequencies, and therefore additive effects of alleles, differ between the training and test populations. Calling VPRS=∑v=1nVav=∑v=1n2pqβv2 the additive variance due to the loci contributing to the PRS in the test population, and hPRS2 the heritability due to the loci contributing to the PRS in the test population, we find the squared correlation coefficient, r², between phenotype and PRS is

(82)

r2=Cov2[PRSj,Pj]Var[PRSj]Var[Pj]

(83)

=(VPRS)2VPRSVP

(84)

=VPRSVP

(85)

=hPRS2.

Thus, if a PRS in a training study correctly estimates all β’s, and those β’s are shared with the test population, and there are no interactions, the squared correlation between the PRS in the test study and the phenotype of those individuals is the heritability explained by the markers included in the PRS in the test population. The careful reader will have noticed that this proof was made substantially harder looking than it might have been because the PRS does not have 0 mean. Nevertheless, nothing about this changes if the training and test populations have differing allele frequencies, differing additive effects, and therefore differing PRS means, so long as they have the same β’s. Recall though, that we began this section by asserting that the reason one was calculating the squared correlation between PRS and phenotype was as a measure of how “useful” or even “good” a PRS was. We have just shown, though, that this measure is the heritability (additive variance over total variance) of the markers used in the PRS in the test population. We have already seen above that if one population has lower minor allele frequency, the additive variance due to the locus will be less, and therefore heritability due to the marker will be less. So, if a PRS is applied to two different populations, and one of them happens to have lower minor allele frequency at most loci, then the PRS will appear to be “doing worse” in the population with the lower minor allele frequencies, even if all β’s are identical and perfectly estimated. Using the numbers from our toy numerical example, we would find that the correlation between PRS and phenotype was about half (916) as “good” in the test study as it was in the training population. More generally, since the squared correlation between PRS and phenotype is a linear function of the variance due to the markers used, which is a function of the allelic frequencies, it becomes intuitive that the “utility” of a PRS is likely to be a linear function of the variance in allele frequency between training and test populations (e.g., Figure 3 from [44]).

Figure 3 The PRS distribution is divided into percentiles. The odds ratio of average members of a percentile group to average members of the population is plotted for a disease with threshold t = 2, ψ ≈ 0.02275 and V_PRS of 1%, 5% and 10% of the total variance. Notice what might appear to be a slightly counterintuitive result. The OR of an individual at the 50^th percentile of PRS is necessarily less than 1. This derives from the fact that the residual variance, 1 − V_PRS is necessarily less than the total, and an individual at the PRS median has a contribution to phenotype from the PRS loci of exactly 0, the same mean as any randomly chosen individual with unknown PRS. However such a random individual has residual variance 1.

2.4 Unification of continuous and binary traits with correction for allele frequency

To simplify the presentation moving forward, we will generally assume that for binary traits OR’s have been converted to effects on the liability scale. After estimating an OR in a training study, we convert the OR to a β value, measured on the liability scale, using prevalence, ψ, and allele frequencies, p and q, from the training study. For any common disorder, a natural approximation would be

For a more rare disorder, where the OR might potentially be very large but ψ small, the simpler approximation of

will likely be more practical. Using one of these estimates for the penetrance of A₁, the conversion to effects on the liability scale finishes with

We should remind ourselves that this used the fact that the variance explained by this locus was relatively small. In this fashion, regardless of whether the trait is continuous or binary, we assume the effect sizes have been estimated as a β which represents the difference in allele effects on the phenotype scale for a continuous trait, or on the liability scale for a binary trait. For all that follows whenever we refer to β we mean this value for either a continuous or binary trait, as appropriate. Furthermore, to overcome the problem of PRS means being a function of allele frequency, we will normalize PRS to have mean 0, by using the appropriate allele frequencies, i.e., the allele frequency in the test population when applied to a test population.

2.5 Medical genetics application

One of the most important practical applications of a PRS is to identify individuals at risk of developing a disease who have not yet been diagnosed with such a disease. If a prediction can be made sufficiently early, there is a greater potential for early interventions that may reduce their risk. In this medical genetics application, a PRS developed in a training study is applied to individuals with unknown phenotype in the test population. This is done to identify individuals with unusually large PRSs with the understanding that individuals with unusually large PRSs have unusually high probability of developing disease. Here we find the potential of this approach is determined by the variance explained by the PRS, V_PRS, in the test population.

To model this, assume a binary disease with threshold t and total prevalence ψ. Continue to assume that the PRS has been normalized to have mean zero and calculated on the liability scale as described above. Also assume that there are a sufficiently large number of nearly independent loci contributing to the PRS so that the score itself is approximately normally distributed with variance explained by the PRS of V_PRS, and therefore with variance due to all other factors 1 − V_PRS. Imagine dividing the test population into PRS percentiles. The first percentile includes anyone with PRS less than, tp₁, and in general tp_i, 1 ≤ i ≤ 99 is given by

where Φ^–1 is the inverse of a standard normal distribution, and tp_i is in units of VPRS. Individuals in the i^th percentile are those with whose PRS is greater than or equal to tp_i−1 and less than tp_i. Let PRS_i be the PRS of a randomly chosen individual with PRS in the i^th percentile.

where ϕ is a standard normal density. These results are simple applications of the mean of truncated normals [45]. Calling the penetrance of a randomly chosen individual in the i^th percentile ζ_PRSi,

where OR_PRSi is the odds ratio of a randomly chosen individual in the i^th percentile relative to a randomly chosen individual in the whole population. Figure 3 plots this odds ratio for varying levels of V_PRS for a disease with threshold t = 2, and prevalence ψ ≈ 0.02275, i.e., a relatively common disease. The utility of this approach is very much an increasing function of V_PRS in the test population. It should also be clear why PRSs are thought to be so promising for medical genetic applications. For a PRS that explains as little as one percent of the total liability for disease, individuals in the highest percentile have an odds ratio over 1.75. For a PRS explaining 10% of the total variance, a PRS in the top 10% gives an odds ratio above two, and a PRS in the top percentile has an odds ratio above 5, which is similar to some of the highest ever estimated single locus contributors to any complex disorder, e.g. the odds ratio of APOE4 for Alzheimer Disease has been estimated to be ≈ 4.6 [46]. Odds ratios this large could justify early life interventions intended to reduce risk. These V_PRS values are from the test population, and are therefore a function of allele frequencies in the test population.

When PRSs are not calculated with a step that normalize the mean to 0, it has been obvious to many that PRSs for binary traits measured as OR’s cannot be meaningfully compared between studies. To find something that could be compared between studies, some investigators [47] measure the difference in PRS means between cases and controls within a study. The difference between case and control means within one test study might be then compared to the difference in case-control means in a second test study. A larger difference between case and control means might be taken as the PRS performing “better.” By now it is perhaps intuitive that the mean difference between case and control PRS is also fundamentally a measure of the variance due to the PRS, and this case-control difference will vary between studies with differing allele frequencies. To understand why this is, continue to assume that PRS has been converted from the OR scale to the liability scale and normalized to have zero mean, as described above. Furthermore, continue to assume a sufficiently large number of nearly independent loci contribute so that the PRS distribution in the test population is well approximated by a normal distribution with variance V_PRS. To find the mean of PRSs in cases (E[PRS_j|D]) and controls (E[PRSj|¬D]), we think of L_j, the liability of individual j, as being the sum of two approximately normally distributed factors, their PRS, PRS_j, and all their remaining (residual) liability, R_j.

where PRS_j is approximately normally distributed with mean 0 and variance V_PRS, and the residual liability R_j is also approximately normally distributed with mean 0 and variance V_R = 1 − V_PRS. Therefore, the mean PRS of diseased individual j with residual liability R_j = x is the expectation of a truncated normal (the PRS distribution) with mean x and variance V_PRS. Thus for a disease with prevalence ψ and corresponding threshold t,

To find the difference in average PRS between cases and controls note that

Figure 4 plots the difference between case/control PRS means as a function of V_PRS for a disease with threshold t = 2. The difference in mean PRSs between cases and controls is a function of the variance due to the PRS, which is a function of allele frequency. Increasing minor allele frequency increases the variance due to the PRS, and therefore increases the mean difference in the PRS between cases and controls. Populations with differing allele frequencies will therefore have varying differences in PRS mean.

Figure 4 For a disease with threshold 2, ψ ≈ 0.02275, the difference in PRS mean between cases and controls is given as a function of the proportion of the variance due to the PRS, V_PRS.

Another approach to evaluating PRS predictive power builds its framework in a machine learning / classification tradition. Here the goal is to find a rule based on the PRS that is used to classify individuals - predict their case/control status from their PRS. In this context the simplest classifier C predicts individual j to be a case if PRS_j > C and otherwise predicts j to be a control. For any given value of C we can find the sensitivity, i.e., the “true positive rate” (TPR), the fraction of individuals who are classified as a case and are actually a case, as well as the specificity, i.e., “true negative rate" (TNR), the fraction of individuals j classified as a control who are actually controls.

An arbitrarily large number of classifiers are possible by setting differing values of C, where any given C implies some specific TPR/TNR combination. The overall utility of this classifier is traditionally presented as a receiver operator characteristic curve (ROC), where 1 − TNR, also known as the “false positive rate,” is given on the x axis and the sensitivity on the y. Figure 5 presents the ROC for various value of V_PRS found in typical (0.01 ≤ V_PRS ≤ 0.1) biobank scale studies of a common complex disease with prevalence a little over 2 percent. Also given is the best possible ROC for that same disease assuming 80% heritability. As is clear, with current study sizes, a classifier based on PRS can achieve a very low false positive rate, almost all individuals classified as cases can be likely cases, by setting a very high threshold on the PRS scale (the top percentile, say), but such a classifier will miss nearly all true cases (poor sensitivity), even for a PRS that explains as much as 10 percent of the total variance. With perfect knowledge of the genetic basis of a trait with 80% heritability, a genetic classifier is still unlikely to be any more useful than a very good “screening tool,” i.e., while it will be possible to identify 95% of individuals likely to develop disease, such a classifier will have a false positive rate of approximately 8%, a value similar to or slightly better than the best prenatal screening tools for Down Syndrome [48].

Figure 5 For a disease with true prevalence ψ ≈ 0.02275, the receiver operator characteristic (ROC) is created by classifying samples with a threshold on the PRS. Samples above the PRS threshold are predicted to be cases. Any such classifier’s utility is determined by the V_PRS.

2.6 The infinitesimal model and its critiques

It has long been suggested [49] that the field of theoretical population genetics developed largely in the absence of any data sufficient to settle many of the key questions the field wished to understand. Perhaps the most basic question considered by the field is “What is the nature of genetic variation contributing to complex traits?” Theory and analysis surrounding this question reached a highly developed state long before there was anything like sufficient data to test most of the key assumptions introduced by that theory. Now that genetic data and estimates of genetic effects on phenotype are widely available, there is, perhaps less direct appeal to theory developed in the absence of this data than there might be. In this particular case, though, somewhat obscure theoretical results can help build our intuition.

In the human genetics community, the start of the GWAS era was accompanied by what was frequently called the “Common Disease, Common Variant Hypothesis” [50, 51]. Advocates and critiques of this hypothesis were, perhaps often unknowingly, echoing debates that had occurred more than a decade earlier in the theoretical quantitative genetics community. Theoretical quantitative geneticists advocating a view most similar to the Common Disease, Common Variant position were sometimes called “neo-Darwinians” [52]. Those arguing an opposing view were sometimes called supporters of the “infinitesimal” model. While the infinitesimal model has its origins with Fisher [53], the first presentation detailing the connections between individual locus population genetics forces (mutation, selection, and drift) with phenotypic level quantitative genetics measures (variance components, heritability etc.) is due to Lande [54]. In Lande’s development of the infinitesimal model, complex traits result from mutation/selection balance under stabilizing selection. The variants that contribute to traits are a combination of alleles of potentially large effect, but with very low frequency, combined with many common alleles with much smaller individual effects. The neo-Darwinians, often spearheaded by Turelli [55], argued that a substantial fraction of genetic variation was likely contributed by high-frequency alleles of large effect, whose frequency was maintained through balancing selection [56]. At the time these debates were most prominent, there was very little direct evidence to support either view. However, the neo-Darwinians developed several interrelated analyses of the infinitesimal model that they believed diminished the likelihood this view was correct. First, the neo-Darwinians used data from experimental breeding studies, in particular mutation accumulation studies, to argue that if the infinitesimal model was correct it requires that the distribution of allelic effects be highly leptokurtic, i.e., that a substantial fraction of the total genetic variance must be contributed by extremely rare alleles of very large effect. A corollary to this analysis showed most traits must be influenced by very many genes, hundreds or thousands generally [56]. The neo-Darwinian school argued that the only alternative to believing in this worldview was to suppose that a substantial fraction of the variation in complex traits was contributed to by common alleles of large effect, maintained by some form of balancing selection.

Over 20 years of GWAS has convinced nearly everyone in the human genetics community that common alleles of individually large effect generally explain little of the variance for most phenotypes studied, be they traits such as height [57], or common disease phenotypes [26]. This is not to say they never exist (e.g. APOE4 discussed above), but that such examples are very rare. The neo-Darwinian prediction that a significant fraction of the variance in most complex traits would be contributed by a few common alleles of very large effect is usually wrong. The natural conclusion, therefore, is that the infinitesimal model is likely largely correct. Nevertheless, if we are to accept this model, the critiques leveled by the neo-Darwinians are no less cogent. If we accept that there are almost no very common alleles of very large effect, then we must believe most complex disease is contributed to by very many loci and that a substantial fraction of the variation is contributed by rare alleles of large effect, likely scattered throughout many loci across the genome.

Placing these intuitions in the modeling context here, we believe that alleles that substantially increase liability for disease will more often than not be the less frequent allele because in this context, as disease per se can be a factor contributing to the stabilizing selection central to the infinitesimal model [54]. In fact, under a wide range of mutational models and effective population sizes, it may be fundamentally impossible to distinguish the effects of stabilizing selection from simple purifying selection on individual alleles [58].

2.7 A population created by a bottleneck

For over 100 years population geneticists have conceptualized a population as a self-contained entity, where one generation promulgates the next via binomial sampling of alleles, a process they generally call Fisher-Wright sampling. Population geneticists generally call alleles “neutral” when their probability of being sampled is independent of their identity and call the change in allele frequency from one generation to the next “genetic drift.” Importantly, the expected frequency of any neutral allele after sampling is the same as its frequency before sampling. Under drift alone, the frequency of a neutral allele is a martingale; a random variable whose expectation the next time it is observed is exactly equal to its current value.

Since the earliest days of fly genetics it was noticed that when a population of flies were kept in a large flask, and then a small number of flies were chosen at random, often by collecting whichever flies happened to move from one flask to the other via the flasks’ “bottleneck,” the frequency of alleles among flies that traveled through the bottleneck were often different. In general, alleles present in the first flask, the “source” bottle, are frequently absent from the second flask, the “destination” bottle. On the other hand, alleles that had been rare in the source bottle are, sometimes, seen in the destination bottle at frequencies much higher than in the source. This phenomenon has long been noted in human disease studies where the fact that many rare, high penetrance, disease alleles were discovered precisely because they were at unusually high frequency in a relatively isolated population [61, 62, 63].

A pair of extant populations may have historically had a source and destination relationship. When this relationship exists we will refer to the extant population which had been the source as the “historical-source” population, and other as the “historical-destination” population. Perhaps counterintuitively, the reason disease alleles are often at unusually high frequency in historical-destination populations is because allele frequency is a martingale unless acted on by strong natural selection. When a destination population was created, the expected frequency of all alleles in the destination population was same as their frequency in the source. However, if the destination population was created by sampling only a relatively small number of individuals from the source population, there is a substantial probability that any rare allele in the source population might not have been sampled and therefore been absent from the destination population, i.e., during a founding event, alleles are lost and are at frequency zero in the destination population. If we think N_f individuals are sampled from the source population to create the destination population then, for an allele A₁ with frequency q in the source population, the probability, ω, that A₁ is entirely absent from the destination population at the time of founding is the probability that its frequency, q′, in the destination population is zero.

The two in front of N_f derives from diploidy. Obviously, the smaller the number of individuals creating the destination population, N_f, the greater the chance an allele is lost. Similarly, the smaller the initial frequency, q, of an allele, the greater the chance it is lost. For any sufficiently rare allele and small bottleneck size, there will be a substantial probability the allele is lost because of the bottleneck.

Nevertheless, allele frequencies are a martingale. Therefore, E[q′] = q. In the destination population, either the A₁ allele will have been lost, in which case its frequency became 0, or it will have been present, in which case its frequency was greater than 0. From the law of total conditional expectation

Therefore alleles that are not lost during the population founding event will have higher average frequency in the destination population than they had in the source population. The rarer the allele and the smaller the number of founders, the greater their frequency will increase, on average, conditional on the allele making it through the bottleneck. Heuristically, we can think of this as a lottery for rare alleles. During the founding event, many rare alleles are losers and will be lost. However, if an allele happens to win, then its allele frequency will likely be higher than it was before.

Thus, at the moment of founding the expected frequency of rare alleles present in the destination population is higher than their frequency was in the source. After this moment of founding, the expected frequency of the alleles in both populations remains a martingale. Thus, at any time after the founding, the expected frequency of an allele present in both populations is higher in the historical-destination population than the historical-source. Of course, the actual frequency in both historical-source and historical-destination populations will drift over time. Nevertheless the average frequency remains higher in the historical-destination population.

That alleles are often at usually high frequency “by chance” in small isolated populations, historical-destinations, has been understood for nearly a century, and generally referred to as a “founder effect” in human genetics. While this phenomenon is widely known, it is, we believe, largely thought of as an isolated “random” effect. That belief starts with the true statement “whenever a bottleneck happens, allele frequencies change, and it is random whether or not the frequency increases or decreases," but fails to recognize that conditional on the allele being present in the historical-destination population, it is more likely to have increased in frequency during the bottleneck than it is to have decreased. Thus, conditional on an allele being observable in the historical-destination population, there is a systematic directionality to its frequency that originated during the founding; bottlenecks increase the frequency (on average) of alleles that survive the event, and this expected frequency in the current historical-destination is whatever it was during the founding. While this analysis contemplates a single founding event, the general logic and framework holds for any subsequent migration between the populations. A rare allele found in only one population before migration will increase in frequency in its new population when it migrates from the larger population to the smaller, and decrease in frequency migrating from the smaller to larger population.

This analysis is precisely correct for neutral alleles and approximates a transient analysis for selected alleles. If we imagine the frequency of the rarer disease allele is being tightly regulated by natural selection in the source population, perhaps because it is under very strong selection and therefore in mutation-selection balance (very large s [64]), then this difference in allele frequency caused by the founding will be present immediately after the bottleneck, but will presumably diminish over time from the effects of selection. In the presence of selection, allele frequency is not a martingale, and no matter how the frequency changed in the destination population, over time sufficiently strong selection will determine the allele frequency in the historical-destination population likely making it more similar to the historical-source if selection acts similarly in both. Under strong selection, the difference in allele frequency between historical-source and historical-destination populations will be on average larger the shorter the time since founding. On the other hand, if selection is weaker, allele frequency will fundamentally be a function of Ns (selection and population size) and the effect of selection on allele frequency will be “relaxed” in the smaller population (but see [65, 58] showing that under many circumstance the dependence on N is often minor). Thus, for “weakly” selected alleles, as long as the effective population size is smaller in the historical-destination population, allele frequency will likely be higher, on average.

Combining this result with the intuition we developed from the infinitesimal model, particularly when variants are ascertained in a case-control framework, we find that we expect the rarer allele, A₁, to increase liability for disease more often than the common allele, A₀, does, i.e., we expect β > 0 more often than β < 0. We expect the larger the effect on disease the lower the allele frequency, on average (β and q are inversely proportional). We think that when we compare two populations, allele frequencies may differ. If one of the populations is something akin to the historical-source population of the other, we expect many rare disease alleles present in the historical-source population will be absent from the historical-destination population. However, when the disease alleles are present in the historical-destination population, they will be on average at higher frequency than they are in the historical-source population. The careful reader will recognize that we have now described a mechanistic insight for the second major correlation assumed in our toy numerical example. If one population is closer to the historical-source and the other the historical-destination, we will expect there to be systematic differences in allele frequency between them. In our toy example, the training population is analogous to the historical-destination population and the test population to the historical-source. In real data, it will never be that all the variants are at higher frequency in the historical-destination than the historical-source, but on average they should be higher. Thus, we have reached the intuition behind the two fundamental correlations demonstrated in the toy example: positive correlation in the sign of β across sites, and positive correlation in allele frequency across loci. While the toy example was intended to be an extreme example of this, we expect the pattern to be general. Allele frequency differences per se ought to contribute to systematic and predictable differences in PRS.

2.8 Allele frequencies differ between studies

With these intuitions in mind, we can now predict what we expect to happen when we train a PRS in individuals from one population but test it in individuals from another. We expect allele frequencies to vary between populations. Thus, even if the allele effect sizes measured as β are the same between populations, the differing allele frequencies will cause each locus to contribute slightly different amounts of additive variance to the trait, have slightly different OR’s etc. Absent any other information, we might imagine that the various allele frequency differences across many loci would “cancel” out, and important averages might be approximately the same. However, if one population is the historical-source population, and other population its historical-destination, there will be consistent and predictable effects on PRSs.

Begin by imagining training a PRS in a historical-destination population and testing it in the historical-source population, and assume in both studies no correction was made for allele frequency. Assume that all alleles included in the PRS are truly disease influencing, have their effects accurately estimated as β’s, and all β’s are the same in both the historical-source and historical-destination population. Of course, if an allele is present in both the historical-source and historical-destination population we imagine the reason will often be because the allele was originally in the historical-source population, and was introduced into the historical-destination population at the time of founding, rather than being the result of multiple different mutational events to the same allele, or introduced from one population into the other by migration after the founding event. This assumption is fundamentally equivalent to assuming that the population founding is relatively recent relative to 4N_e generations, where N_e is the effective population size of the historical-source population. For humans with effective population sizes somewhere between 10, 000 and 100, 000 and a generation time of at least 15 years, this is equivalent to assuming the founding events happened less than, say, 150, 000 – 1, 500, 000 years ago, which is a near certainty. Therefore, when an allele is present in both historicalsource and historical-destination populations, we suspect those alleles were present in the historical-source population at the time of the founding event, survived the bottleneck, and will have, on average, higher frequency in the historical-destination population than in the historical-source population. This will not always be true, but we imagine it is true more often than not. In this scenario, on average we expect q in the training study (the historical-destination population) to be on average larger than in the test study. These are the conditions of the toy numerical example. If no correction has been made for allele frequency the PRS will have mean

Because we expect disease alleles to have positive β on average, see Section 2.6, increasing allele frequency q implies an increasing mean PRS. Thus, with no allele frequency correction, we expect the mean PRS to be larger in the training set (historical-destination population) than in the test set (historical-source population). If the disease is being contributed to by sufficiently rare alleles or the bottleneck included sufficiently few people, the means of the PRS could be so far apart that the two PRS distributions do not overlap (our numerical example). Not only will the means of the two PRS distributions differ, but the additive variance accounted for by the PRS will also differ between the two studies,

and we see the variance accounted for by the PRS will be higher in the training (historical-destination) population than in the test (founding, historical-source) population. Recalling that if we measure the effectiveness / utility / goodness of the PRS by the squared correlation between PRS and known phenotype, the training population will have higher squared correlation because it has a larger V_PRS. Similarly if we measure effectiveness as something like the difference between case / control PRS means, this too is an increasing function of V_PRS. In short, if we train a PRS in a historical-destination population and then test it in the historical-source population we will conclude the PRS “works better” in the training population than it does in the test population even if all genetic effects are identical. Moreover, the means of the PRS distribution will be lower in the test set, with the two PRS distributions potentially showing no overlap. Taken together these two observations could easily lead to the conclusion that the PRS “does not work well” in the test population.

Conversely, if we train the PRS in the historical-source population and test it in the historical-destination population, we might reasonably expect that some of the disease alleles present in the historical-source population will be absent from the test population. The effect of this will be discussed in detail below. However, as a thought experiment let us assume that all disease alleles used for training also happen to be in the test population, perhaps by chance, or perhaps because alleles absent from the test set were not included in the training set due to the fact that the reagent used for genotyping (a genotyping chip, say) in the training population only included SNPs known to be present in many populations. In this circumstance we would see the opposite result. The PRS would appear to have higher mean in the test population than in the training population, and the correlation between PRS and phenotype would be larger in the test population. We might appear to have the paradoxical result that a PRS trained in one population works better in a different population. If the sample sizes were vastly different in the training and test studies, one might be tempted to hypothesize that the PRS working “better” in the test study had something to do with the precision of the β estimates, but here we see the phenomena is entirely consistent with simply training in an historical-source population and testing in a historical-destination one.

2.9 An allele present in one study is absent from another

We expect many rare alleles present in the historical-source population to be lost in the bottleneck and therefore absent from the historical-destination population. On the other hand, we also expect there to be some recently arisen alleles, post-bottleneck alleles, that are only in one of the two populations because they have only recently come into existence, and there has been insufficient migration between the founding and historical-destination populations post-bottleneck for the allele to be found in both places. Thus, we can think of alleles present in one population but absent from the other as belonging to two classes. The first class are relatively older alleles that existed before the founding and by chance were lost in the historical-destination population. These are alleles that could only contribute to disease in the historical-source population. The second class are relatively younger alleles that have arisen after the bottleneck in only one of the two populations, and could be found in either the founding or historical-destination populations. Assuming the founding event is relatively recent compared to population size, we might reasonably expect the first class to considerably out number the second.

If we train our PRS in the historical-destination population, alleles contributing to disease in the historical-destination population but absent from the historical-source population – class 2 alleles found only in the historical-destination population – will exacerbate the trends described above. These are alleles contributing additive variance in the training set but with q = 0 in the test set. They will decrease the mean PRS in the test set, and decrease the correlation between PRS and phenotype. These alleles make a bad situation worse, if you think of lower correlation between PRS and phenotype as a sign the PRS is worse.

If we train our PRS in the historical-source population, all class 1 alleles and class 2 alleles found only the historical-source population will work to reverse the effects of alleles present in both populations. These alleles have q > 0 in the training population, but q = 0 in the test population, and therefore increase the mean PRS in the training population relative to the test, and contribute to V_PRS only in the training population, and therefore increase the correlation between PRS and phenotype only in the training group. Depending on how many of these alleles exist in the training set, we can imagine the mean PRS of the training population being higher, lower, or even approximately equal to the test population. Similarly we can imagine the correlation between PRS and phenotype could be higher, lower, or equal.

Thus, if we train a PRS in a historical-destination population and test it in a the historical-source population we expect the mean of the PRS to be lower in the test population, as will the correlation between PRS and phenotype. On the other hand, if we train in the historical-source population and test in the historical-destination any result is possible. The PRS means and correlation between PRS and phenotype might be higher, lower or equal in the two populations. What is certain, though, is that training in the historical-destination population and testing in the founding (historical-source) population will generally make the differences between the two populations appear to be larger than training in the founding and testing in the historical-destination. This is because class 1 alleles (those lost during the founding event) act in the opposite direction as alleles present in both populations.

2.10 Differing LD between studies

Linkage disequilibrium (LD), the non-random association between alleles at different loci, often reinforces this general pattern, but in some circumstances can reverse it. In general, LD between sites causes correlation in effect sizes. Call β_v and β_w the difference in allelic effects at sites v and w in some real population under study with linkage disequilibrium D ≠ 0 between these sites. Imagine another population, identical in every way to the first, but where sites v and w are not in LD (D = 0), and call βv˜ and βw˜ the effect sizes in this imaginary population without LD. As shown in the first paper in this series

where p_v, q_v, p_w, q_w are the allele frequencies at sites v and w. In a Falconer inspired derivation of these results we would call βv˜ and βw˜ the “true” effect sizes, and β_v and β_w the “estimated” effects in the presence of LD. We also show that the standard measure of LD, D, can be derived as a “haplotypic” covariance, and as a result the LD measure r² is, in fact, a squared correlation coefficient, so 0 ≤ r² ≤ 1 The immediate implication of this is that

We note that the equality, D=pvqvpwqw=pvqv, can only hold when q_v = q_w, and the alleles are in “perfect LD.” If the minor allele frequencies differ between the two loci, then necessarily r² < 1. If q_v > q_w, say, then −q_wq_v ≤ D ≤ p_vq_w.

In what is likely a common scenario, consider a situation where SNP w has an effect on phenotype independent of LD, β̃_v ≠ 0, but SNP v does not, β̃_v = 0. Further assume q_v ≥ q_w, and site v is included in the construction of a PRS, but site w is not, possibly because site v is included on some widely available genotyping array, and site w is not. In a Falconer construction, we would say site w is the “real” effect site, but site v is an “LD surrogate” included in the PRS. For simplicity assume that there are no other sites in LD with v and w with any LD-independent effects. Thus, if we had included site w in the PRS, and everything were perfectly estimated βw=βw˜. However, βv=Dβw˜pvqv. To understand how this influences PRSs, we must consider two situations, corresponding to positive and negative D between these sites. It is certainly possible to enumerate all four haplotypes, and move through the possible allelic configurations methodically, but an intuitive approach is probably more helpful. Recalling that D can be written as a covariance, if haplotypes containing the rare allele at both loci, A_v1A_w1, occur more often than expected by chance given their individual allele frequencies, q_wq_v, then D will be positive, because the covariance between allelic states is positive. If such haplotypes occur no more often than expected by chance, D = 0. If such haplotypes are less common than expected by chance D < 0. Put most simply, if there is positive correlation between allelic states, D is positive. If negative correlation, D is negative.

The immediate implication is that when D is positive, the sign of β_v will be the same as β_w. When D is negative the signs flip. If we subscribe to an “infinitesimal” model of complex phenotypes, we expect when |β_w| is large, β_w will generally be positive and q_w small. However, if we include site v in the PRS not w, the sign of β_v depends on the sign of D, and the general patterns we describe above hold if and only if haplotypes containing both rare alleles are found more often than expected by chance. Otherwise, β_v will more often than not have a negative sign, and all the described relationships reverse. Whether one expects the sign of D to be positive or negative is a complex question involving details of the coalescent tree in this region, and requires a somewhat detailed understanding of population sizes and how they change over time and related features of the demographic ancestry of these individuals. Importantly, when populations differ significantly in LD, the sign of D between the v and w could be opposite, resulting in the sign of β_v appearing to flip between populations [66]. There may be no easy rule of thumb here. The situation becomes even more complex when there are multiple sites with effects all of which are in LD with v, but also not included in the PRS. Multiple sites in LD will be explored more broadly in the third paper in this series.

Despite the sign of D being somewhat challenging to estimate as a general rule, several effects of including site v in the PRS rather than site w are immediately clear. First, unless q_v = q_w and the two sites are in perfect LD, then |β_v| < |β_w| because |D|pvqv<1. Including the “wrong” site in the PRS has the effect of making the effect size, β closer to zero. Similarly the additive variance due to locus v will be smaller than locus w, 2pvqv(βv)2≤2pwqw(βv)2(pvqwpwqv). This simple fact may have significant implications for understanding “missing” heritability.

Returning to the question of training a PRS in individuals from one population, but testing in a different population, we find including the “wrong” site to have effects that often accentuate the previously described patterns. On average, a bottleneck has no effects on LD; haplotype frequencies do not change, on average. However, if two sites have significant LD between them in the historical-source population, then necessarily “recombinant” haplotypes (haplotypes created only by recombination between the sites) are far more rare in the historical-source population than either of the minor alleles themselves. Being the most rare haplotype, the chance that recombinant haplotypes are completely lost during the bottleneck is greater than any other haplotype. Thus, for sufficiently small bottleneck sizes, it will not be uncommon for all recombinant haplotypes to be lost during the bottleneck, and LD to increase. Of course, conditional on the recombinant haplotypes surviving the bottleneck, LD will decline, but as a rule of thumb (which we will see in data below), LD will generally be greater in the historical-destination than in the historical-source population.

Larger LD in the historical-destination population implies that |β|’s at site v (the site included in the PRS) will be larger in the historical-destination population, as will the additive variance at site v independent of any allele frequency differences. A PRS trained in the historical-destination population should have a higher mean PRS, and better correlation between PRS and phenotype, from LD alone, reinforcing all the patterns seen from changing allele frequency described above (Section 2.8). A PRS trained in the historical-source population will have lower mean PRS than in the historical-destination and lower correlation between PRS and phenotype due to LD, unless site w is absent from the historical-destination population, in which case the patterns described in Section 2.9 hold. Of course, LD attenuates all this by a factor of Dpvqv. Thus, as a basic rule, including sites “only in LD” with the “real” effect sites will often reinforce the differences in the PRS between founding and historical-destination populations similarly to the changes due to allele frequency alone, However, if the sign of D is negative, opposite patterns can be found.

3.1 Review of challenging observations

To understand the most challenging observations in PRS analysis, begin by making the following broad-stroke assumptions about human demographic history: African populations are historical-source to all; Latino individuals are a likely mixture of historical-source and historical-destination populations; the broader European population is historical-source relative to Ashkenazi and Greenland populations. In general, rare alleles present in both training and test studies cause training in historical-destination populations and testing in historical-source to have lower PRS mean and variance in the test population than in the training population. Training in the historical-source population and testing in the historical-destination leads to more complex relationships driven by the relative contributions of rare alleles present in both populations versus rare alleles found only in the training population. “Cross-population” training studies (training studies including individuals with diverse ancestries) can benefit from larger sample sizes, but when used with SNP selection methods that are LD aware [82], can also substantially eliminate the complications described in Section 2.10. Cross-population training and SNP selection are also less likely to chose variants for inclusion in the PRS that are absent from some of the training ancestries. Thus, cross-population training is likely to result in better testing performance [83, 44] by picking better SNPs: those whose estimated effects are less influenced by LD and those less likely absent from different populations. Keeping these broad patterns in mind, cross population application of PRS seems predictable and understandable.

When PRSs for height trained in European-biased genetic studies are tested in African samples, both the PRS mean and variance appear to be considerably lower in Africans [5, 6]. Similarly, PRSs for schizophrenia trained in European-biased genetic studies have decreased mean when tested in Africans [5, 7]. PRSs derived from SNPs trained on European populations appear to underestimate (lower PRS mean) the risk of cardiovascular disease in African individuals [9]. Following a similar pattern, a PRS of only 11 SNPs in Greenlanders (historical-destination) had a V_PRS = 0.162 for LDL-cholesterol, compared to 2 million variants needed to explain 22% of the variance in LDL-cholesterol in Europeans (source) [14]. PRS models of acute lymphoblastic leukemia [13] trained in non-Latino Whites (historical-destination) explained a greater proportion of variance when tested in non-Latino Whites (historical-destination) compared to when tested in Latinos (mixture of historical-source and historical-destination). Conversely, a PRS for breast cancer trained in the broader European population (source) and tested in Ashkenazi Jewish women (historical-destination) demonstrated good discrimination (an increasing function of V_PRS) but significant overestimation of breast cancer risk (higher PRS mean) in the Ashkenazi Jewish [18]. These are all exactly the expected observations when PRS includes a significant contribution from rare alleles present in both populations and at higher average frequency in the historical-destination population.

On the other hand, when PRSs for type II diabetes and asthma are trained in multi-ethnic cohorts (historical-source and historical-destination mixture, but a substantial fraction of historical-destination individuals), the PRS mean and variance are both larger in African populations (historical-source) than any other populations tested (historical-destination) [10, 5, 11]. A PRS for acute lymphoblastic leukemia [13] trained on multi-ancestry GWAS data (historical-source and historical-destination, but historical-destination as the largest contributor) explained equal proportions of variance when tested in non-Latino Whites (historical-destination) and in Latinos (historical-source and historical-destination), which was significantly more than the variance explained by the PRS trained non-Latino Whites (historical-destination, only). Because PRSs trained in historical-source populations can include rare alleles not found in historical-destination populations, we interpret these results as suggesting that both type II diabetes and asthma have significant contributions from rare alleles only found in historical-source populations, as does acute lymphoblastic leukemia, but to a lesser extent.

The observations of Somineni et al. [12] may appear perplexing at first, until one notes that unlike most studies here effect sizes were measured on the liability scale, and PRS mean was adjusted for allele frequency differences. Thus, the differences in PRS between training and test populations were more subtle than generally observed. Nevertheless, a PRS trained in African Americans (source) yielded a higher odds ratio in the top percentile of polygenic risk individuals when tested in African Americans (source) then when tested in Europeans (historical-destination), suggesting the presence of African-specific rare alleles whose effects are most prominent in the PRS tail. On the other hand, a PRS trained in Europeans yielded only slightly higher odds ratios in the PRS tails when tested in European vs African Americans, consistent with the slightly higher V_PRS in Europeans likely contributed to by alleles present in both populations but at lower frequency in the historical-source population.

Wang et al. (2020) [20] found that linkage disequilibrium and minor allele frequency differences between ancestries could be used to account for 70-80% of the change in predictive power, a function of V_PRS, of PRSs for body mass index and type II diabetes trained in European populations and then tested in both European and African populations. Because this study only included allele frequency and LD measured in common alleles, it is not entirely surprising that there may exist effects of rare alleles not completely accounted for by common allele frequency and LD differences. On the other hand, with whole genome sequencing in admixed individuals [21] using the observed allele frequencies in the tested individuals directly, there was virtually no evidence for any difference by ancestry in effect sizes measured as β’s for any variant (although confidence intervals were often substantial). While the original authors suggest this result might reflect the presence of GxE interactions, the simple minded interpretation is that when the effects of LD and allele frequency can be fully accounted for, genetic effects measured as β’s are the same in all humans. While the examples given here are without question “cherry picked,” as it is nearly impossible not to cherry pick results from a literature with 1000s of examples, it is challenging to find observations that are entirely inconsistent with the patterns described here. We believe that a simple parsimonious explanation of the PRS literature is that β’s are very likely to be generally consistent across human populations. Allele frequency and LD can exhibit complex patterns that may vary considerably by ancestry.

3.2 Recommendations

Many of the challenges associated with PRS analysis are largely solved by simply moving to the disease liability scale and being careful to account for allele frequency differences, and consequently differences in PRS mean and variance. For a quantitative disease phenotype (like blood pressure), estimate β at all SNPs in a linear regression, or linear mixed model, in whatever is thought to be a best practice way. For a binary phenotype, estimate the odds ratio in logistic regression or other best practice manner. For the binary phenotype convert the OR’s to penetrances, and then to their allelic effects, α’s, and the difference in allelic effect, β. Thus, we have the “same” measure of effect β on the quantitative phenotype scale for a quantitative disease, or on the liability scale for binary disease. At this stage, one could thin SNPs in a typical way, select sites with large β and remove others in LD with selected sites. Alternatively, one could attempt to use the β estimates and LD regression coefficient matrix to estimate β’s, i.e., the effect size in the absence of LD. This sort of approach is fundamentally similar to [84] and related BLUP (best linear unbiased predictor) estimators of individual variant effect sizes. One would include any SNP with substantial β˜, perhaps determined by 2pqβ˜2 being above some threshold, or a similar choice criterion that is proportional to statistical confidence that β˜≠0. Select all SNPs that meet this criterion, regardless of the LD relationship between them, and use β˜ for the SNP’s effect, as this is the measure of effect size independent of LD. From a technical standpoint this procedure may be extremely challenging. Since many SNPs may have very similar LD patterns with their neighbors, the LD regression matrix may be singular and unable to be inverted. Even if technically invertible, the matrix may be very “stiff,” highly sensitive to small perturbations and numerical roundoff error. This approach also assumes that higher order LD is insignificant compared to pairwise LD. Finally, if alleles contributing to effect in an LD-independent way are not genotyped in the training study, β˜=0 for all SNPs actually typed, but there are untyped SNPs in LD with these causing |β| > 0, this approach could lead to extremely challenging to interpret results, and is unlikely to fully resolve the challenges LD induces in PRS analysis. The availability of whole genome sequencing among individuals in the training study combined with LD estimated directly from these subjects will certainly help to alleviate these problems. There is much room for development of new, robust estimation routines at this stage. Cross-population training studies can assist in limiting the effects of LD. The greater the variation in LD among the training individuals, the less stiff the LD matrix will be. Put more simply, variation in LD always assists in fine-mapping variants [82].

Nevertheless, if we assume that the estimates are “reliable,” then the SNPs selected and their estimated effect sizes in a training study are available to be applied to any test study, regardless of allele frequency or LD differences, under the simple assumption that β˜’s are the same (one of the few measures that could be the same between study populations). However, to do so and avoid the complications and confusions associated with differing PRS means, construction of the PRS score for individual j in the test population proceeds naturally using the allele frequencies q_v in the test population,

Of course, PRS_j is fundamentally an individual measure; it is a property of person j. While we might abstractly think of person j as coming from a particular population, for many people their recent ancestry is likely to be complex, and there is no single value of q_v that can be said to reflect their “population.” On the other hand, for any individual with data available from whole genome sequencing, or modern genotyping arrays, there exist several robust methods to infer ancestry [85] that usually report measures that can be viewed as an estimate of the proportion of individual j’s ancestry that comes from given reference populations. Using these estimates as weights combined with reported allele frequencies from gnomAD [75], a “personalized” estimate of q_v can be used, q_vj, that might be thought of as the allele frequency among individual j’s recent ancestors. In such a manner, with any robust estimates of β˜, a PRS can be constructed for individual j regardless of j’s personal ancestry, that should have mean ≈ 0, and be largely comparable between any individuals j.

Evaluation of a PRS should avoid comparing training to test studies as any reasonable evaluation technique will be a function of V_PRS, which is in turn a function of allele frequencies, and potentially personalized ancestry of the individuals, and therefore not directly comparable if the training and test studies include individuals with differing ancestries. Instead the performance of the PRS in the test population should most naturally be compared to its expectations assuming the β˜ values are the same in all individuals. For a quantitative phenotype with total variance V_P, the squared correlation, r², between the phenotype of individual j and PRS_j is expected to be VPRSVP. Thus, one could test whether r² is bigger than zero (establishing that at least some of the SNPs included in the PRS have effects similar to the training set), or more quantitatively if r² is significantly different from VPRSVP. A failure to reject equality implies there is no evidence for differing genetic effects between the training and test studies.

For a binary phenotype, the quantity most interesting to the medical geneticist is the probability an individual will ultimately develop disease as a function of their PRS. If D_j is an indicator that individual j will develop disease such that D_j = 1 if the person ultimately develops disease and D_j = 0 otherwise, then

where t is the disease threshold and Φ is a standard normal cumulative distribution. If disease phenotype were known, one could calculate measures such as ROCs for various PRS thresholds etc. Of course, as is obvious from the above, the expectation for these measures is a function of V_PRS which is a function of allele frequency and ancestry, and thus not comparable between studies. A natural evaluation technique might be to compare the mean PRS in cases and in controls, or area under the curve (AUC) of the ROC, and ask if the difference in means (AUC) is greater than zero (0.5) - establishing that at least some of the β’s are similar between training and test set, or more quantitatively, if the difference in means (AUC) is consistent with V_PRS. Failure to reject equality would argue there is no substantial evidence for a difference in genetic basis of disease in the two studies under comparison.

3.3 Application of PRS recommendations

To understand the effect of these recommendations we use training data from [86] with allele frequency and effect size on Crohns Disease (CD) reported in the GWAS Catalog [87]. The training study consisted primarily of European ancestry individuals (86, 640) meta-analyzed with a much smaller number of East Asian Individuals (9, 846). As the test study we use data from [12] with 1, 335 African Americans with CD and 1, 644 control African Americans. For each individual in the test study, we first estimate the global fraction of their ancestry which derives from each continent using the published gnomAD 4 ancestry PCs [75]. We next calculate a “personalized" allele frequency at each site for each individual in the test study by weighting the gnomAD published allele frequencies for each continental group by the estimated fraction of that individual’s ancestry coming from the continent. Assuming an overall prevalence of CD of 0.00061, we convert reported effect sizes to β’s on the liability scale via Equation 90. We construct the PRS with Equation 130, using the individual’s personalized allele frequency at all sites. The V_PRS is calculated using the average over all individuals of the personalized allele frequencies at each site. The associated penetrance of each individual is found via Equation 134, and then converted to an OR relative to a randomly chosen individual. To compare these recommendations to the “standard” methods of presenting PRS, we calculate the a PRS with Equation 44 and the associated OR with Equation 47. First we consider 128 genome-wide significant SNPs identified by the training study authors. All of these SNPs have been identified after fine-mapping and represent the authors’ best estimate of the “casual” variants in each LD region. If these SNPs are directly contributing to phenotype, and the β’s are the same between ancestries, the only difference between training and test studies should be due to allele frequency differences.

Values in Table 1 for the training study use only the reported summary statistics, and should be thought of as akin to their expected values given the reported effect sizes. Values for the test study are the observed values from the test study. Given the allele frequencies at these sites, the V_PRS in the test study was estimated as 0.047, very close to the observed values for these individuals, and approximately 10% smaller than the training study. When OR is estimated in the recommended fashion, interpretable results occur. The OR in controls is approximately 1 on average, and 1.3 in cases, suggesting the combined effect of these 128 variants which explain 4% to 5% of the variance have combined to increase the odds of the average case by a factor of 1.3. There is no difficulty in interpreting differences between the training and test study. Using the standard approach leads to values for the PRS and OR that are not immediately interpretable. While the standardly calculated OR in the training study appears slightly different than the OR in test controls, it is likely close enough to avoid any serious interpretation challenges.

Table 1 128 Genomewide Significant SNPs selected after fine-mapping. Means are over individuals studied.

Using more SNPs with less confident estimates of effect size in the training study creates far more difficult interpretation challenges with the standard method. To see this we now select any SNP with a reported p < 10^–3 thinning with the simplest possible algorithm. Any selected SNP within a megabase of any other selected SNP is placed in a set. Each set is thinned to a single member by selecting the SNP with the lowest p-value. This rather rudimentary algorithm results in the selection of 574 SNPs (Table 2). Because the p-value threshold is well above multi-test corrected significance, we expect many of the selected SNPs to have no actual effect on phenotype. Their estimated β’s are pure noise. Also, because there was no careful attempt to account for the effects of LD, we expect the selected SNPs to sometimes be in LD with the actual effect sites in complex ways that might differ between studies. Given the personalized allele frequencies, the V_PRS was estimated to be 0.119 in the training study, a value certainly inflated by the inclusion of variants not actually contributing to the trait.

Table 2 574 SNPs with p < 10^–3. Means are over individuals studied.

Here inclusion of sites with less evidence of effect and poor accounting for LD leads to biased estimates using our recommended procedure. The control mean is clearly larger than 0. It should be slightly less than 0. The observed V_PRS in both cases and controls appears biased upward. The estimated OR’s also appear biased upwards. Nevertheless, the effect of including 574 SNPs can be directly compared to using only 128. Using more SNPs increased the difference in PRS and OR means between cases and controls. This strongly suggests that some of the newly included sites contribute to phenotype in a meaningful way. A simple interpretation ensues: some but not all of the newly included sites likely contribute to phenotype, but their effect sizes may be poorly estimated, and LD incompletely accounted for. On the other hand, the standard approach gives values that are exceedingly hard to interpret. The mean PRS in the training study is substantially lower than the control mean from the test study. In fact the difference between the training population and test controls is larger than the difference between test cases and test controls. OR’s show a similar pattern. More pointedly, the OR’s when using 574 SNPs seem to be on a wholly different scale than when using 128 sites. This is presumably because those OR’s are relative to a hypothetical individual with A₀A₀ genotype at all loci. It is easy to believe that any investigator might find interpretation of the standard PRS values very challenging.

3.4 Combining PRS with environmental factors

One could further elaborate on this analysis by explicit inclusion of known environmental factors, to either increase the predictive power of a PRS or perhaps to test for a PRS by environment interaction [88, 89]. As discussed in [23], testing for interactions for continuously distributed factors is a complex estimation problem often dependent on the precise distribution of the continuous factors and their scale relative to each other. However, for any combination of discretely distributed factors, testing for interaction is straightforward and not dependent on the phenotypic or factor distribution. Here we illustrate a methodology for discretizing a PRS and combining with a binary environmental variable, but this basic procedure can be extended to any number of categorical environmental variables.

Assume a binary environmental variable E, with two states E₀ and E₁, with average effect on phenotype P of ϵ₀ and ϵ₁, respectively for the two states. Estimate the effect of substituting environment E₁ for E₀, β_E = ϵ₁ — ϵ₀, in a standard way (linear regression or linear mixed model, say), accounting for population structure and other important covariates. For binary phenotypes, estimate the OR of E₁ to E₀ in a logistic regression / mixed model again accounting for the important covariates, and convert this OR_E to β_E measured on the liability scale. Thus, for both a continuous or binary phenotype we have the effect measured as β_E on the phenotypic / liability scale. If f_E0 and f_E1 = 1 – f_E0 are the frequencies individuals are in state E₀ and E₁, then the variance due to this environmental factor is Ve=fE0fE1βE2. Calculate the combined PRS + Environment score, PRSE_j for individual j as

where e_j is 0 if individual j has environmental state E₀, and 1 otherwise. Of course for a binary phenotype

The simplest method to test for a PRS by environment interaction is to create discrete PRS bins. For sufficiently large sample sizes, one might considered PRS percentile (100) bins as described above. For smaller sample sizes quartile (4) or decile (10) bins might be more appropriate. There are fundamentally two sources of genetic by environment interaction. The first is correlation between genotypic state and environmental state. This can be directly tested in this framework by asking if the frequency of state E₁ differs between PRS bins. This test could be asked comparing, perhaps, the first and last bin, or as an “omnibus-test” by asking if the observed variance in the frequency of E₁ across bins was greater than expected by chance. If there is significant variation in environmental state between PRS bins, then an interaction exists.

An interaction induced by a correlation between genetic and environmental states is not the only form of interaction that is possible. Regardless of whether or not correlation in state exists between PRS and environment, we say a gene by environment interaction exists if the mean phenotype of an individual with a combination of PRS and environment states differs from the sum of the mean given the environment plus the mean given the PRS. Adapting the notation from [23], and assuming percentile bins, and environmental state E_m, m ∈ {0, 1},

where fE1i is the frequency of environment E₁ in PRS bin i. If there is no detectable variation in the frequency of environments across bins, then we assume fE1i=fE1 for all bins i. VdIgePRS,E is the deviation variance caused by non-additive interactions between PRS and environment E. One can test for the significance of VdIgePRS,E in some standard way. Of course, one would naturally estimate VdIgePRS,E in some standard way. Of course, one would naturally estimate βE=ϵ1−ϵ0 the average phenotype given the environment and PRS bin, in the same manner used to estimate β_E = ϵ₁ – ϵ₀, with a linear or logistic regressions controlling for presumed important covariates. As usual, convert OR’s to effect sizes on the liability scale for binary phenotypes. In this manner we can observe and quantify interactions between PRS and environment whether induced by correlation in state, deviations from additivity, or both for quantitative or binary phenotypes.