rs12409788

rs1576517

rs151240271

rs12752436

rs76864380

rs6586443

rs35213023

rs34910942 * Choose first CHS (Southern Han Chinese) population and press `Calculate`. It shows that only two of these SNPs have MAF>5% and four are monomorphic in the Table shown, i.e., show only one allele in these data. There are only 4 haplotypes present in the population out of theoretically possible $2^5=32$. * Change population to FIN. There are 5 haplotypes present and one SNP is monomorphic. Compare SNPs 1 and 6. They have the same pattern of variation across the haplotypes: Whenever one has **T** also the other has **T**. They are completely associated with each other and, consequently, in a GWAS they would give exactly the same results. Also triple 5, 8 and 9 is similarly fully correlated. * Change population to LWK (Kenya). Now we see 8 haplotypes, no monomorphic SNPs and only the SNPs 5 and 8 are completely linked with each other. This example shows how some variants in a narrow region of 2500 bps are highly correlated with their neighbors when it comes to the haplotype patterns that exist in any one population. Typically, African populations show more diversity than Asian and European populations due to more recent genetic bottlenecks that have reduced variation in the history of the latter two. But also in African populations, the haplotype patterns are very much restricted compared to what they would be if the alleles at nearby loci would be anywhere close to independence at the population level. The human genome has a **haplotype block structure** (slide 5), which suggests that the recombination process shuffles the chromosomes mainly at only a relatively small number of sites on the genome, so called *recombination hotspots*. Between these hotspots, there are larger regions of genome with low recombination rates (slide 6). Consequently, * Phasing, especially in the low recombination areas, can be done accurately by collecting a large number of genotyped individuals and using a probability model to estimate which are the most likely haplotypes in this population. Since only a few haplotypes exist in the population, we can get a fairly reliable phasing for individuals by sharing haplotype information among them, and we simultaneously also get accurate estimates for the population frequencies of the haplotypes. Existing high quality reference haplotype data from the target population also improves phasing. Phasing software: [Beagle5](https://faculty.washington.edu/browning/beagle/beagle.html), [Eagle2](https://data.broadinstitute.org/alkesgroup/Eagle/), [SHAPEIT3](https://jmarchini.org/shapeit3/). * SNPs in the same haplotype block tag each other well because each allele of a particular SNP sits on only a limited number of haplotypes in that block. In particular, by typing the set of the most informative tag-SNPs, we can already cover also a lot of neighboring genetic variation that we don't directly genotype. This was a property that boosted the first GWAS: in order to see statistical associations with phenotypes, it is enough that we manage to genotype any tag-SNP of a causal SNP. * The exact causal SNPs (or other variants) are difficult to pinpoint based on the statistical information alone since so many SNPs are so highly correlated with each other. This is the **fine-mapping problem**. ### 7.2 Linkage disequilibrium (LD) We know from empirical data that the human genome has a haplotype block structure. This correlation structure arises because when a new variant (assume to be SNV) is introduced to the population as a mutation, it emerges on one of the existing haplotypes of the population. Because haplotypes are inherited as large segments of grandparental haplotypes, with often just 1 or 2 recombination events per chromosome, the new mutation will be, for many subsequent generations, passed on to the next generation with the same alleles that were close by to it on its orginal haplotype background. However, as more generations pass, recombinations cut down the original haplotype background to smaller and smaller shared segments between the descendants of the original haplotype, and hence the correlation at the population level between the variant and its neighbors tends to decrease over time (slides 8-9). Consider any two SNPs at distinct sites of the genome. Suppose that we had perfect haplotype information. Then we could observe how often their alleles 1 are on the same haplotype in the population and compare that to the expected proportion if the two variants were independent. If allele 1 frequencies in population are $f_k$ and $f_l$ at SNPs $k$ and $l$, and if the frequency of the haplotype 1-1 is $f_{kl}$, then we define * disequilibrium coefficient $D_{kl}=f_{kl}-f_k\cdot f_l$, which, theoretically, is 0 if the two loci are independent, also said to be in **linkage equilibrium**, LE. If the loci are not independent ($D_{kl}\neq 0$), then we say that the loci are in **linkage disequilibrium** (LD). Naturally, we never know the frequencies exactly in population and we use sample estimates to compute $D_{kl}$. Hence we will always observe some non-zero value for $D_{kl}$ in our sample. The interest is in the magnitude of $D_{kl}$. We can measure the amount of LD also by using the Pearson's correlation coefficient $r_{kl} = \textrm{cor}(a_k,a_l)$, where $a_k$ and $a_l$ are the indicators of the alleles at SNPs $k$ and $l$ on the same haplotype. It can be shown that $r$ and $D$ are related through * $r_{kl} = \frac{D_{kl}}{\sqrt{f_k(1-f_k)f_l(1-f_l)}}.$ Note that, theoretically, both $r$ and $D$ are 0 if and only if the loci are in linkage equilibrium. In what follows, we are interested in $r$, because it determines the statistical consequences of LD on GWAS analyses. A statistical way to think about LD is that if SNPs $k$ and $l$ are in LD, then by observing for a particular haplotype its allele at SNP $k$, we also gain more information about its allele at SNP $l$ than what we had solely based on the population allele frequencies at SNP $l$. #### 7.2.1 Generating two-locus haplotype data Let's follow [Vukcevic et al. (2011)](https://onlinelibrary.wiley.com/doi/full/10.1002/gepi.20576) who derived formulas for the 4 haplotype frequencies in the population, given the MAFs at the two SNPs (1 and 2), and LD as measured by $r = r_{12}$ between the SNPs. (More comprehensive treatment of the subject is available in [Damjan Vukcevic's D.Phil thesis, Oxford, 2009](https://ora.ox.ac.uk/objects/uuid:8f89593e-a4ab-4df0-b297-74194be7891c).) Suppose that SNP1 has alleles **a** (minor) and **A** (major) and SNP2 has alleles **b** (minor) and **B** (major) and that $f_a \leq f_b.$ First, in order that the allele frequencies and the given level of LD, as measured by $r$, correspond to a haplotype distribution, the following inequalities must hold: $$ -\sqrt{\left(\frac{f_a}{1-f_a}\right) \left(\frac{f_b}{1-f_b}\right)} \leq r \leq \sqrt{\left(\frac{f_a}{1-f_a}\right) \bigg/ \left(\frac{f_b}{1-f_b}\right)}. $$ For example, correlation between two SNPs can be really high ($r\approx 1$) only if both have very similar MAFs, and in order for correlation to be highly negative ($r\approx -1$), both MAFs must be close to 0.5. This is because if we were to assign the minor allele at both SNPs to haplotypes, then the only way that we can make correlation $\approx 1$ is that **a** and **b** are almost always on the same haplotype and **A** and **B** are almost always on the same haplotype. This is only possible if $f_a \approx f_b$ because otherwise there are either some extra **a**s that don't find any **b** to pair with and therefore must pair with **B**, or there are some extra **b**s that don't find an **a** to pair with and therefore must pair with **A**, both of which will lead to a reduction of correlation from $\approx 1$. To get a high negative correlation for the minor alleles, each minor allele must be paired up with the major allele at the other locus, and this is possible only if MAF is close to 0.5 at both loci. Note that if we change the allele coding at one locus from allele 1 being the minor allele to allele 1 being the major allele, that simply changes the sign, but not the magnitude, of the corresponding $r$. The haplotype frequencies, and samples from them, can then be computed using the R-function `geno.2loci()` given below. ```{r} geno.2loci <- function(n, r, mafs, return.geno = TRUE){ #INPUT: # n, individuals # r, correlation coefficient between the alleles on the same haplotype of the two loci # mafs, MAFs of the two loci #OUTPUT: # if return.geno=TRUE: n x 2 matrix of GENOTYPES of n individuals at 2 loci # if return.geno=FALSE: (2n) x 2 matrix of HAPLOTYPES of n individuals (2n haplotypes) at 2 loci stopifnot( r >= (-1) & r <= 1 ) stopifnot( length(mafs) == 2 ) stopifnot( all(mafs > 0) & all(mafs <= 0.5) ) stopifnot( n > 0) #Label SNPs and alleles so that a and b are minor alleles and freq a <= freq b. #At the end, possibly switch the order of SNPs back to the one given by the user. f.a = min(mafs) # maf at SNP1 f.b = max(mafs) # maf at SNP2 #With these parameters, the possible LD coefficient r has values in the interval: r.min = max( -1, -sqrt(f.a/(1-f.a)*f.b/(1-f.b)) ) r.max = min( 1, sqrt(f.a/(1-f.a)/f.b*(1-f.b)) ) #c(r.min,r.max) #Check that r is from this interval if(r < r.min | r > r.max) stop(paste0("with these mafs r should be in (",r.min,",",r.max,")")) # Alleles SNP1: A (major) and a (minor); SNP2: B (major) and b (minor). # Compute conditional probabilities for allele 'a' given allele at locus 2: q0 = f.a - r*sqrt(f.a*(1-f.a)*f.b/(1-f.b)) #P(a|B) q1 = f.a + (1-f.b)*r*sqrt(f.a*(1-f.a)/f.b/(1-f.b)) #P(a|b) #Compute the four haplotype frequencies: f.ab = f.b*q1 f.aB = (1-f.b)*q0 f.Ab = f.b*(1-q1) f.AB = (1-f.b)*(1-q0) f = c(f.ab,f.aB,f.Ab,f.AB) f #These are the haplotype frequencies in the population. haps = matrix(c(1,1,1,0,0,1,0,0), nrow = 4, ncol = 2, byrow = T) #4 haplotypes in the population. #Generate data for n individuals where each individual is measured at these two SNPs: hap.ind = sample(1:4, size = 2*n, replace = T, prob = f) #There are 2*n haplotypes, 2 for each individual if(mafs[1] > mafs[2]) haps = haps[,2:1] #Whether to change the order of loci? #Either make genotype matrix by summing the two haplotypes for each individual... if(return.geno) X = haps[hap.ind[1:n],] + haps[hap.ind[(n+1):(2*n)],] if(!return.geno) X = haps[hap.ind,] #...or return haplotypes as such. return(X) } ``` Let's check that this seems to work and gives a correlation of $r$. ```{r} params = matrix(c(0.3,0.3,0.5, #each row has: r, MAF at SNP1, MAF at SNP2 0.8, 0.4, 0.3, -0.5, 0.34, 0.4), byrow = T, ncol = 3) n = 10000 for(ii in 1:nrow(params)){ X = geno.2loci(n, r = params[ii,1], mafs = params[ii,2:3], return.geno = FALSE ) #haplotype level data print(c(cor(X[,1],X[,2]), colSums(X)/nrow(X))) } #and compare to params above ``` It works, so let's use it to demonstrate how haplotypes with different LD look like. Start from high LD (r=0.9), then intermediate LD (r=0.66) and finally no LD (r=0). MAFs at the two SNPs must be similar for the high LD case. (Here we use 0.4.) ```{r, fig.width = 9} n = 50 rs = c(0.9, 0.666, 0.0) #correlation values mafs = matrix(c(0.4,0.4, 0.4,0.5, 0.4,0.4), byrow = T, ncol = 2) cbind(r = rs, mafs) #check that we put correct values in matrix par(mfrow = c(3,1)) par(mar = c(2,4,4,1)) for(ii in 1:nrow(mafs)){ X = geno.2loci(n, r = rs[ii], mafs = mafs[ii,], return.geno = F) #return haplotypes X = X[order(X[,1]),] #order so that haps with major allele at SNP1 come first image(X, yaxt = "n", xaxt = "n", xlab = "", ylab = "", col = c("floralwhite","dodgerblue"), cex.main = 1.3, main = paste0(2*n," haplotypes, observed r=",signif(cor(X[,1],X[,2]), 3) )) abline(h = 0.5, lwd = 1.5, lty = 2, col = "black") axis(2, at = c(0,1), labels = c("SNP1","SNP2"), cex.lab = 1.3 ) } ``` In the first case, we see that allele **A** (white) at SNP1 almost always indicates allele **B** (white) at SNP2, and ALSO allele **a** (blue) indicates allele **b** (blue) at the other SNP. Thus, knowing the allele at either of the SNPs predicts the allele at the other SNP with very high accuracy. These two SNPs are in **high LD**: population parameter was $r=0.9$, the observed value from the sample is shown in the title. In the second case, MAF at SNP1 is lower than at SNP2 and even though blue at SNP1 predicts well that SNP2 also has blue allele, also many white alleles at SNP1 go together with a blue allele at SNP2. Hence LD, as measured by $r$, is less than in the first case. This case demonstrates that SNPs with very different MAFs cannot have very high correlation because there is no way to make all blue match blue AND all white match white when there are different numbers of blue and white at the two SNPs. These SNPs are clearly in LD with each other, although not as strongly as the SNPs in the first case. In the third case, the SNPs are independent in the population ($r=0$ in simulation), and also in the observed data there is little correlation between the loci. Thus, by knowing the allele at SNP1 does not help predicting what is the allele at SNP2 on the same haplotype. The best predicton we can do statistically is to guess the allele at SNP2 based on the population allele frequencies, and ignore which allele was observed at SNP1, since it does not give additional information over the allele frequencies in population. These SNPs are **not in LD** with each other. Let's also check what happens to the estimate of $r$ if we estimate it at the genotype level, not at the haplotype level, as we have done so far. Thus, we'll generate haplotypes as previously, but before computing $r$, we collapse the haplotypes into genotypes (by `return.geno = TRUE`) and then compute correlation between the genotypes at the two loci. ```{r} params = matrix(c(0.3,0.3,0.5, #each row has: r, MAF at SNP1, MAF at SNP2 0.8, 0.4, 0.35, -0.5, 0.34, 0.4), byrow = T, ncol = 3) n = 10000 for(ii in 1:nrow(params)){ X = geno.2loci(n, r = params[ii,1], mafs = params[ii,2:3], return.geno = TRUE ) print(c(cor(X[,1],X[,2]), colSums(X)/2/nrow(X))) } #and compare to params above ``` This seems to estimate the haplotype level correlations even though the input data are genotypes. This is good to know as the data we observe in practice are at the level of genotypes and phasing it to haplotypes is not always straightforward. **Example 7.2.** We can also get LD estimates for human populations from [LDlink](https://ldlink.nci.nih.gov/?tab=home). * Choose LDhap and consider two SNPs rs7837688, rs4242382 in the CEU population. You'll see the haplotype distribution in 1000G CEU data. There are 3 haplotypes (out of all 4 possible) and one is quite rare (1%). So these are almost in perfect LD where allele at one SNP tells the allele at the other. * While you could compute $r$ from the info given by LDhap, you can also query it directly from LDpair tab, which gives $r^2=0.8791$, i.e., $r=0.938.$ (Slide 10.) * Look at the same SNPs in LWK population. All 4 haplotypes are present and there is no LD in terms of $r^2$. Often African populations have less LD than non-African since the latter have been through more recent bottlenecks. Such bottlenecks increase LD because when there is a relatively few haplotypes that found the subsequent generations then those haplotypes will determine the LD patterns. (Slide 11.) * Go to LDproxy tab and look rs7837688 in the Chinese CHB population. You'll see the Figure and list of the variants that are in highest LD with rs7837688. Note also how the high LD is concentrated within a haplotype block that is between recombination hotspots (region between high peaks of recombination rate cM/Mb). (Slide 12.) * LDmatrix (slide 13) can be made for a user given set of SNPs. **Example 7.3.** Let's plot an LD matrix for 74 variants around *APOE* gene on chr 19 using 1000 Genomes Finnish samples (coordinates GRCh37). ```{r,fig.width = 6} path = "https://www.mv.helsinki.fi/home/mjxpirin/GWAS_course/material/APOE_1000G_FIN_74SNPS." haps = read.table(paste0(path,"txt")) info = read.table(paste0(path,"legend.txt"),header = T, as.is = T) dim(haps) #rows haplotypes, cols SNPs info[1:3,] #info for first three SNPs haps[1:3,] #data for first three haplotypes R = cor(haps) #LD matrix is the pxp correlation matrix n.cols = 100 layout(matrix(c(1,2), nrow = 1), widths = c(9,1)) par(mar = c(2,2,3,1)) #plot LD matrix using image image(abs(R), col = heat.colors(n.cols), breaks = seq(0, 1, length = n.cols + 1), asp = 1, xaxt = "n", yaxt = "n", bty = "n", main = paste("|r| of ", dim(R)[1],"SNPs around APOE")) par(mar = c(2,1,3,1)) #plot scale for colors plot.window(xlim = c(0,1), ylim = c(0,n.cols)) points(x = rep(1,n.cols + 1),y=(0:n.cols), col = heat.colors(n.cols + 1), pch = 15, cex = 2) axis(4, at = c(0, n.cols / 2, n.cols), labels = c(0, 0.5, 1) ) ``` LD matrices from two Finnish cohorts show highly similar LD structure to each other on slide 14. ### 7.3 Effect of LD on GWAS results Let's then experiment with a GWAS setting where SNP1 is a causal variant with causal effect size of $\lambda_1 = 0.2$ and MAF = 0.2 whereas SNP2 does not have a causal effect ($\lambda_2=0$) and has MAF = 0.4, on a quantitative phenotype whose variance is 1 in population. LD between the SNPs is $r_{12}=0.6$. So far on this course we have used $\beta$ to denote the effect size of allele 1. Reason why we now use $\lambda$ rather than $\beta$ is because we want to separate the **causal effect** ($\lambda$) from the **marginal effect** ($\beta$), as will become clear soon. Let's make 1000 simulations of such setting, using $n=1000$ individuals in each simulation, and let's do the typical single-SNP GWAS on each variant and collect the effect estimate $\widehat{\beta}_l$ its SE and P-value for both variants $l=1,2.$ ```{r, fig.height=4} n.iter = 1000 n = 1000 r = 0.6 mafs = c(0.2, 0.4) lambda = c(0.2, 0) #causal effects of each SNP res = matrix(NA, ncol = 6, nrow = n.iter) #6 cols: beta1, SE1, P1; beta2, SE2, P2 colnames(res) = c("beta1","SE1","pval1","beta2","SE2","pval2") for(iter in 1:n.iter){ X = geno.2loci(n, r, mafs, return.geno = T) #generate 2-locus genotypes y = X %*% lambda + rnorm(n,0, sqrt(1-var(X %*% lambda)) ) #var(y) = 1, res[iter, 1:3] = summary(lm(y ~ X[,1]))$coeff[2, c(1,2,4)] #collect beta, SE, P-val of SNP1 res[iter, 4:6] = summary(lm(y ~ X[,2]))$coeff[2, c(1,2,4)] #collect beta, SE, P-val of SNP2 } ``` Let's plot effect estimates $\widehat{\beta}_1$ and $\widehat{\beta}_2$ against each other and as separate distributions. Let's add vertical lines at the expected mean values: $\lambda_1$ for SNP1 and $r_{12}\,\lambda_1\,\sqrt{\frac{f_1(1-f_1)}{f_2(1-f_2)}}$ for SNP2 (explained later where this comes from). ```{r, fig.height = 5, fig.width = 9} par(mfrow = c(1,2)) plot(res[,"beta1"], res[,"beta2"], xlab = expression(beta[1]), ylab = expression(beta[2]), main = paste("r =",signif(cor(res[,"beta1"], res[,"beta2"]), 3))) abline(0,1) plot(density(res[,"beta1"]), xlab = expression(beta), col = "black", main = "", xlim = c(0,0.4), ylim = c(0,8)) lines(density(res[,"beta2"]), col = "red") abline(v = lambda[1], lty = 2) abline(v = lambda[1]*r*sqrt(mafs[1]*(1-mafs[1])/mafs[2]/(1-mafs[2])), lty = 2) legend("topright", col = c("black","red"), lwd = 1, legend = c(1,2) ) ``` We see that there is a clear correlation in GWAS effect estimates between SNPs 1 and 2, and that correlation seems to match well with the LD correlation between the SNPs ($r_{12}=0.6$). The distributions of the effect estimates shows that $\widehat{\beta}_1$ is centered on the true causal effect $\lambda_1$ and $\widehat{\beta}_2$ is centered at value $r_{12}\,\lambda_1 \sqrt{\frac{f_1(1-f_1)}{f_2(1-f_2)}}$ (and not at 0, which would be the true *causal* effect of SNP2; derivation on slide 16). Let's look at the P-values: ```{r, fig.height = 5, fig.width = 9} par(mfrow = c(1,2)) plot(-log10(res[,"pval1"]), -log10(res[,"pval2"]), xlab = "-log10P1", ylab = "-log10P2") abline(0,1) plot(density(-log10(res[,"pval1"])), xlab = "-log10", col = "black", main = "", ylim = c(0,0.35)) lines(density(-log10(res[,"pval2"])), col = "red") legend("topright", col = c("black","red"), lwd = 1, legend = c(1,2) ) ``` Also the P-values are correlated and SNP2 has P-values that are clearly much smaller than under the null even though SNP2 does not have a causal effect on the phenotype. We also see that there are a few cases where SNP2 actually has lower P-value (higher -log10 P-value) than SNP1 even though SNP1 is the causal variant that is driving the genotype-phenotype association here. Consequently, * When SNPs 1 and 2 are LD-friends (in LD with each other), and only SNP1 is causal, then the causal effect of SNP1 will leak to the marginal effect of SNP2 in the sense that the standard one-SNP-at-a-time analysis at SNP2 is estimating marginal effect $$\beta_2 = r_{12} \lambda_1 \sqrt{\frac{f_1(1-f_1)}{f_2(1-f_2)}},$$ rather than the true causal effect of SNP2 (which would be 0). * This means that if we had not genotyped SNP1 but only genotyped SNP2, we would still see the association signal in this region, but its effect size would be reduced by a factor of $r_{12}\sqrt{\frac{f_1(1-f_1)}{f_2(1-f_2)}}$, and its NCP would by reduced by a factor of $r_{12}^2$, compared to those values that we would estimate at the true causal variant. * Hence, when we see peaks on Manhattan plots, we cannot be sure whether we have the causal variation itself included in our analysis, or whether we are just tagging it by our SNPs. In principle, the possible hidden causal variation does not need to be a single SNV, it can be a structural variant, or a haplotype-level effect, for example. * The top variant in terms of P-value does not need to be the causal variant because of tagging and statistical sampling effects. In particular, when there are many variants in high LD with the causal(s) then the association region looks messy (slide 15) and we can't know for sure which are the causal variants by looking at the marginal effect sizes and P-values. #### 7.3.1 Relation between causal effect $\lambda$ and marginal effect $\beta$ Consider a region with $p$ SNPs with their causal effects in vector $\pmb{\lambda} = (\lambda_1,\ldots,\lambda_p)^T.$ This means that the region contributes to the phenotype by $\sum_{l=1}^p x_{il}\lambda_l.$ When we do a standard GWAS, that includes one variant at a time in the regression model, we estimate the marginal effects $\beta_l$ rather than the causal effects $\lambda_l$. Above we saw how the non-causal variant ($\lambda_2=0$) had non-zero marginal effect size ($\beta_2\neq 0$), whose value we could predict from the causal effect $\lambda_1$, $r_{12}$ and allele frequencies at the two SNPs. This relationship generalizes to an arbitrary number $p$ of SNPs in a region and the marginal effect at SNP $l$ is $$\beta_l = \sum_{k=1}^p r_{kl} \lambda_k \sqrt{\frac{f_k(1-f_k)}{f_l(1-f_l)}} = \lambda_l + \sum_{k\neq l} r_{kl} \lambda_k \sqrt{\frac{f_k(1-f_k)}{f_l(1-f_l)}},$$ where the first sum is over all variants in the region, and, in the second formulation, we have separated the variant's own causal effect from the sum (as $r=1$ between the variant and itself) (Slide 17). Thus, if $l$ is not in LD with any other (causal) variant, then the marginal effect $\beta_l$ equals to the causal effect $\lambda_l$. Otherwise, the marginal effect at SNP $l$ combines the causal effect of $l$ with the causal effects from all other variants, weighted by their LD with the SNP $l$. In particular, when a non-causal SNP $l$ is in high LD with one or more causal variants, then we estimate its marginal effect as non-zero even when the SNP has no biologically causal effect on the phenotype. This phenomenon explains why Manhattan plots show piles of nearby SNPs that together reach highly significant P-values: a few of them may be truly causal, but most of them are only in LD with those few causals, and most of the non-zero marginal effects simply reflect the causal effects of their LD-friends. We want the formulas to be compact and easy to remember, and we are not happy with the above formula with this respect, so we introduce some scaling to help us here. #### Scaled genotypes and scaled effects. Suppose that MAF at SNP $l$ is $f_l$ and individual $i$ has genotype $x_{il}$. We define **scaled genotype** of individual $i$ as $$x_{il}^* = \frac{x_{il}-2f_l}{\sqrt{2f_l(1-f_l)}}.$$ Thus, while $x_{il}$ had three possible values 0, 1 or 2 independent of MAF, also $x_{il}^*$ has three possible values, but these values are chosen so that the average scaled genotype in the population is 0, its standard deviation and variance is 1, and, consequently, the numerical difference between two genotype classes gets larger as MAF gets smaller. This is shown in the Table below, where the columns present the numerical values of the scaled genotypes corresponding to the genotypes 0,1 and 2, for four values of MAF $f$. |genotype | $f=0.5$ | $f=0.10$|$f=0.01$|$f=0.001$| |---------|---------|---------|---------|---------| |0 |-1.41 | -0.47 |-0.14 |-0.045 | |1 |0 |1.89 | 6.96| 22.33| |2 |1.41 | 4.24| 14.07| 44.70| We further define the scaled effects $\beta_l^*$ by running the (linear or logistic) regression model using the scaled genotypes: $$Y \sim \mu^* + X^*\beta^*.$$ This model yields exactly the same P-value for $\beta^*_l$ as the regression using original genotypes gives for $\beta_l$, but the effect size that is being estimated is now $\beta^*_l = \beta_l \sqrt{2f_l(1-f_l)}$ and its standard error is $\textrm{SE}^*_l = \textrm{SE}_l \sqrt{2f_l(1-f_l)}$ (so z-score $\beta_l/\textrm{SE}_l = \beta^*_l/\textrm{SE}^*_l$ doesn't change). **Allelic effect** $\beta_l$ describes the change in phenotype per each copy of allele 1, whereas **scaled effect** $\beta_l^*$ describes change in phenotype per each standard deviation unit of the genotype variable in the population. This linear transformation does not change the statistical model fit ($R^2$ or P-value) at all. In particular, we need not run any new GWAS in order to estimate the scaled effects, but we can simply transform the existing allelic effect estimates $\widehat{\beta}_l$ to the scaled effect estimates $\widehat{\beta^*_l}$ by multiplying by $\sqrt{2f_l(1-f_l)}.$ Note also that, for the quantitative traits, the square of the scaled effect ${\beta_l^*}^2 = 2f_l(1-f_l)\beta^2$ is directly the phenotypic variance explained by the SNP. We do the same with causal effects where **scaled causal effect** $\lambda_l^* = \lambda_l \sqrt{2f_l(1-f_l)}.$ With scaled effects, the relationship between marginal and causal effects is simply $$\beta^*_l = \sum_{k=1}^p r_{kl}\lambda_k^* = \lambda_l^* + \sum_{k\neq l} r_{kl}\lambda_k^*,$$ where $r_{kl}$ is the correlation between SNPs $k$ and $l$, and that correlation is exactly the same whether we consider original genotypes or the scaled genotypes. Thus we have simplified our earlier formula by replacing the complicated allele-frequency-ratio-square-root by a few *s. It is clear that when we are interested in the biology behind the phenotype, then we are only interested in the causal variants and their causal effects. How can we estimate these causal effects? So far our GWAS results have only been about marginal effects $\beta$. #### 7.3.2 Joint model Let's write the phenotype as an additive combination of the (causal) effects of all $p$ variants in the genome: $$y_i = \mu + \pmb{x}_i\cdot \pmb{\lambda} + \varepsilon_i = \mu + \sum_{l=1}^p x_{il}\lambda_l + \varepsilon_i,$$ where $\pmb{x}_i$ is the row $i$ of the full genotype matrix $\pmb{X}$ that contains the $p$ genotypes of individual $i$, and $\pmb{\lambda}$ is the vector of the (unknown) causal effects of the $p$ variants. For binary phenotype, $y_i$ is to be interpreted as the log-odds of the disease for individual $i$ (like in logistic regression). The difference between this **joint model** of all variants and the **marginal model** of a single variant ($y_i = \mu + x_{il}\beta_l + \varepsilon_i$) is that the joint model accounts for the effects of all variants on the phenotype simultaneously whereas the marginal model only considers one variant at a time. Consequently, an estimator of an effect from the marginal model will estimate the marginal effect $\beta_l$ whereas an estimator of the joint model is estimating the causal effects $\pmb{\lambda}$ simultaneously. A standard GWAS applies $p$ separate marginal regression models, one for each variant, and if we collect the marginal effects into a single vector $\pmb{\beta}=(\beta_1,\ldots,\beta_L)^T$, then by the formulas above, we have a linear relationship between the marginal and causal effects: $$\pmb{\beta^*} = \pmb{R}\pmb{\lambda^*} \quad \textrm{ or equivalently } \quad \pmb{\lambda^*} = \pmb{R}^{-1} \pmb{\beta^*},$$ where $\pmb{R}$ is the $p\times p$ **LD-matrix** whose element $(l,k)$ is the pairwise correlation $r_{lk}$ of variants $l$ and $k$. We give this formula for scaled effects rather than the allelic effects because the same formula for allelic effects involves ratios of allele frequencies for each pair of variants, which would complicate the notation compared to this version where the differences in allele frequencies have been scaled away. But remember that we can always get from the scaled effect of variant $l$ to the allelic effect of $l$ by dividing by $\sqrt{2f_l(1-f_l)}$. These considerations suggest two ways to estimate the causal effects $\pmb{\lambda}$. 1. Use a multiple regression model where phenotype is regressed on multiple genotypes simultaneously. In R, this is done by simply using a matrix $X$ of all genotypes in the formula like `lm(y~X)` or `glm(y~X, family="binomial")`. 2. Take the marginal GWAS effects and transform them by the inverse of the LD matrix: $\pmb{\lambda^*} = \pmb{R}^{-1} \pmb{\beta^*}$. Let's try the first one with the settings of our previous example, where we had MAFs 0.2 and 0.4, causal effects $\pmb{\lambda} = (0.2, 0)$ and $r_{12}=0.6$ in the population. Let's fit both the marginal models and the joint model and plot the estimates of $\beta$s and $\lambda$s. ```{r, fig.height = 5, fig.width = 9} n.iter = 500 n = 1000 r = 0.6 mafs = c(0.2, 0.4) lambda = c(0.2, 0) #causal effects of each SNP res = matrix(NA, ncol = 12, nrow = n.iter) colnames(res) = c("beta1","SE.b1","pval.b1","beta2","SE.b2","pval.b2", "lambda1","SE.l1","pval.l1", "lambda2","SE.l2","pval.l2") for(iter in 1:n.iter){ X = geno.2loci(n, r, mafs, return.geno = T) #generate 2-locus genotypes y = X %*% lambda + rnorm(n, 0, sqrt(1 - var(X %*% lambda)) ) #var(y) = 1, res[iter, 1:3] = summary(lm(y ~ X[,1]))$coeff[2, c(1,2,4)] #collect beta, SE.b1, P.b1 of SNP1 res[iter, 4:6] = summary(lm(y ~ X[,2]))$coeff[2, c(1,2,4)] #collect beta, SE.b1, P.b1 of SNP2 joint.coeff = summary(lm(y ~ X))$coeff res[iter, 7:9] = joint.coeff[2, c(1,2,4)] #collect lambda1, SE.l1, P.l1 of SNP1 from joint model res[iter, 10:12] = joint.coeff[3, c(1,2,4)] #collect lambda2, SE.l2, P.l2 of SNP2 from joint model } par(mfrow = c(1,2)) #1st row (2 panels): #Plot marginal effects as separate distributions #Add blue lines to the true causal effects lambda #Add black line to the expected marginal effects of beta2: #r*sqrt(maf1*(1-maf1)/maf2/(1-maf2))*lambda1 for(ii in 1:2){ eff = c("beta","lambda")[ii] title.txt = c("marginal model", "joint model") plot(density(res[,paste0(eff,"1")]), xlab = eff, col = "black", main = title.txt[ii], ylim = c(0,9), xlim = c(-0.2,0.5)) lines(density(res[,paste0(eff,"2")]), col="red") abline(v = lambda[1], lty = 2, col = "blue") abline(v = lambda[2], lty = 2, col = "blue") abline(v = lambda[1]*r*sqrt(mafs[1]*(1-mafs[1])/mafs[2]/(1-mafs[2])), lty = 2, col = "black") legend("topright", col = c("black","red"), lwd = 1, legend = c(1,2) ) } ``` We see that the estimates from the joint model are centered on the true causal effects (blue lines, right panel), whereas the marginal effect of the non-causal SNP2 is estimating the value $r_{12}\,\lambda_1 \sqrt{\frac{f_1(1-f_1)}{f_2(1-f_2)}}$ (black line, left panel). It also seems that the distributions of estimates are wider for $\lambda$s than for $\beta$s. That is because when the two variables are correlated, the regression model tries to figure out how to divide the total effect between the two correlated variables, which increases uncertainty about how large the effect should be *for any one of them separately*. This uncertainty becomes larger the higher the magnitude of the correlation. At the extreme, for two perfectly correlated variables, there is no way to separate their effects statistically; the standard regression model breaks down if you try to include the same variable twice there. The price to pay for an unbiased estimate of the causal effect by the joint model is an increased uncertainty of the effect estimate compared to the marginal model. It can be shown that $$\textrm{SE}\left(\widehat{\lambda_1}\right) = \frac{\sigma_\varepsilon}{\sqrt{2n\left(1-r_{12}^2\right)f_1(1-f_1)}} \approx \frac{\textrm{SE}\left(\widehat{\beta_1}\right)}{\sqrt{\left(1-r_{12}^2\right)}}.$$ Thus, when SNPs are independent ($r_{12}=0$), the precision of $\widehat{\lambda}$s is the same as for $\widehat{\beta}$s and when correlation approaches $\pm 1$, then the joint model becomes more and more unsure how the joint effect should be split between the two candidates and SE of $\widehat{\lambda}$s grows towards infinity. In practice, it is not useful to include very highly correlated variables in a regression model simultaneously because their effects cannot be separated from each other (except with extreme amounts of data). It follows that if we look at the P-value from the Wald statistic of SNP1, then we would expect to see lower P-values for $\beta_1$ than for $\lambda_1$ because while both estimate the correct causal effect 0.2, SE of $\widehat{\beta}_1$ is lower than SE of $\widehat{\lambda}_1$, which leads to lower P-values for $\beta_1$ than for $\lambda_1$ on average. This is not a problem in standard GWAS, where we first identify the regions where variants show high association based on the marginal models, ($\widehat{\beta}$s and their P-values), and then we turn to a more detailed analysis of those regions using some versions of a joint model. A joint (linear or logistic) regression model of several SNPs gives unbiased estimates of the causal effects of the SNPs assuming that ALL causal variants of the region are included in the model. This is a nice property, but in the same time these standard regression models are not suitable for large scale variable selection among thousands of candidate variables and they do not handle well highly correlated variables. Therefore, we can only use them for joint models of a handful of variants at a time. Next we will look how this is typically done for a GWAS region by stepwise forward selection. And later we will see how more recent fine-mapping models combine sparse causal models with efficient search algorithms. ### 7.4 Stepwise forward selection Let's aim to reduce the set of $p$ SNPs in the region into a (much) smaller subset $S$ of SNPs that can already statistically explain the association pattern that we see in the region. Technically, this means that we want to find as small a subset $S$ as possible, for which the joint regression model $$Y \sim \pmb{X}_S$$ explains the data $Y$ already so well that any extended model $$Y \sim \pmb{X}_S + X_a,$$ with one (or more) additional SNP(s) $a \notin S$ included, would not lead to a statistically better fit to the observed data. (Here $\pmb{X}_S$ is a matrix whose columns contain the genotypes of the SNPs included in set $S$.) There could be many versions of such a set $S$, and here we start by trying to find one of them. We will demonstrate this with two data sets generated using the *APOE* region data on 186 Finnish haplotypes from the 1000 Genomes project, and we call these data sets as APOE.1 and APOE.2. **Example 7.4: APOE.1 data.** Let's generate 1000 individuals by sampling their haplotype pairs from the reference pool with replacement and adding them up to get individual level genotypes. This reflects what we would observe from genotyping. ```{r} path = "https://www.mv.helsinki.fi/home/mjxpirin/GWAS_course/material/APOE_1000G_FIN_74SNPS." haps = as.matrix(read.table(paste0(path,"txt"))) dim(haps) #rows haplotypes, cols SNPs n = 1000 X = haps[sample(1:nrow(haps), size = n, repl = T),] + haps[sample(1:nrow(haps), size = n, repl = T),] #n individ. ``` Let's then choose three causal variants: 7, 37 and 51, and give them all a causal effect of $\lambda = 0.3$. Let's plot the causal and marginal effects for the region. ```{r,fig.height=4, fig.width=9} c.ind.1 = c(7, 37, 51) # causal SNP indexes for APOE.1 data p = ncol(X) lambda = rep(0, p) lambda[c.ind.1] = 0.3 f = colSums(X)/2/nrow(X) # allele 1 freqs f.sc = sqrt(2*f*(1-f)) #scaling factor for genotypes R = cor(X) # R is computed from genotypes beta = (R %*% (lambda * f.sc)) / f.sc # from lambdas to betas par(mfrow = c(1,2)) par(mar = c(4.5, 4.5, 2, 0.5)) plot(lambda, xlab = "SNP", ylab = expression(lambda), col = "dodgerblue", pch = 19, ylim = c(-0.5,0.5), main = "causal effects") plot(beta, xlab = "SNP", ylab = expression(beta), col = "red", pch = 1, ylim = c(-0.5,0.5), main = "marginal effects") ``` We see that many marginal effects are larger than any causal effects because some SNPs are in high LD with several causals and therefore their marginal effects are a weighted sum of the causals. In particular, some non-causals have larger marginal effects than some causals, and therefore it could happen that the top-SNP in this region is not a causal variant. Let's generate a trait with variance 1, do a standard marginal GWAS model, and highlight the causal variants in the Manhattan plot of this region. ```{r,fig.height = 4, fig.width = 9} set.seed(18) #21 y.1 = X %*% lambda + rnorm(n, 0, sqrt(1 - var(X %*% lambda)) ) #trait for APOE.1 data res.1 = apply(X, 2, function(x){summary(lm(y.1 ~ x))$coeff[2,c(1,2,4)]}) #3 rows: beta,SE,pval par(mar=c(4.5, 4.5, 2, 1)) plot(-log10(res.1[3,]), xlab = "SNP", ylab = "-log10 P", main="APOE.1 marginal analysis") points(c.ind.1, -log10(res.1[3,c.ind.1]), cex = 2, lwd = 1.4, col = "dodgerblue") #highlight causals ``` The top variant is one of the causals (37) but there are non-causals also nearby in P-values. Could we fit a joint model for all the variants in the region and get a clearer picture of which might be the causal ones? ```{r, fig.height = 4, fig.width = 9} joint.coeff = summary(lm(y.1 ~ X))$coeff[-1,c(1,2,4)] #joint model isn't able to return results for all SNPs nrow(joint.coeff) #how many coefficients we got out of 74? par(mar = c(4.5,4.5,2,1)) plot(-log10(joint.coeff[,3]), xlab = "SNP", ylab = "-log10 P", main = "APOE.1 joint analysis") ring.i = which(row.names(joint.coeff) %in% paste0("XV",c.ind.1)) # indexes of the causals present points(ring.i, -log10(joint.coeff[ring.i,3]), cex = 2, lwd = 1.4, col = "dodgerblue") ``` Unfortunately, the standard linear regression is unable to compute statistics for 16 out of 74 variants because of high correlations -- we miss also one of the causals -- and the joint model gives next to useless results for the remaining SNPs as the effects are split between all correlated SNPs simultaneously, whereas in reality there are only three causals here. Technically, the joint model would work, if we had a large enough sample size, but for many SNPs with high correlations, those sample sizes would need to be unrealistically large. ```{r, eval = F, include = F} #Write files for FINEMAP data.name = "APOE_1" FM.path.1 = paste0("FINEMAP/FINEMAP_",data.name) FM.path.2 = paste0("FINEMAP_",data.name) write(R, ncol = ncol(R), file = paste0(FM.path.1,".ld")) write.table(cbind(rsid = as.numeric(1:p), chromosome = rep(19,p), position = info[,2], allele1 = info[,3], allele2 = info[,4], maf = apply( cbind(info[,5], 1-info[,5]), 1, min), beta = res.1[1,], se = res.1[2,]), file = paste0(FM.path.1,".z"), quote = F, row.names = F) cat("z;ld;snp;config;cred;log;n_samples\n", file = paste0(FM.path.1,".master")) cat(paste0(FM.path.2,c(".z;",".ld;",".snp;",".config;",".cred;",".log;"),collapse = ""),n,"\n", sep ="",file = paste0(FM.path.1,".master"), append = T) ``` A more robust way is to use **stepwise forward selection** (also called **iterative conditioning**, or **conditional analysis**) to build iteratively a joint model for a small subset $S$ of SNPs by choosing one SNP at a time to be included in $S$. This approach is made popular by [GCTA](https://cnsgenomics.com/software/gcta/#Overview)'s Conditional & joint (COJO) analysis of GWAS results. The algorithm works as follows ``` initially S is empty repeat until all P-values outside S > threshold add SNP with the lowest P-value to S update P-values of all SNPs 'l' outside S using joint model Y ~ X.S + X.l end repeat ``` After the algorithm finishes, $S$ contains a set of SNPs whose P-values in the joint model with their predecessors in $S$ are below the given threshold, and there are no additional SNPs with this property outside $S$. Idea is that $S$ will only contain one of the SNPs that are tagging similarly the causals because when two such SNPs are simultaneously in the model, then they don't have low P-values anymore, which is saying that the second one is not needed *after the first one is already in the model*. This *conditioning* on the current contents of $S$ gives the procedure the name of iterative conditioning / conditional analysis. The name we use here, *stepwise forward selection*, means that we proceed *forwards* from the empty model, step-by-step, where at each step one SNP can be selected to be included in S. (*Backward selection* would start with all variables in the model and would proceed by eliminating them one-by-one. This would not be a feasible strategy when there are a lot of candidate variables available.) On the other hand, if there are two causal variants in the region, then the P-value of a SNP tagging primarily the second causal variant will remain low after we condition on the top-SNP that tags primarily the first causal signal. Hence both SNPs may be included in $S$. Sometimes the P-value of the second SNP may even get much lower after the first SNP is included in the model, if the two SNPs were masking each other in their marginal effects. Slide 18 has an example of such a case from a [Parkinson's disease locus 4q22/*SNCA*](http://hmg.oxfordjournals.org/content/20/2/345.long). Thus, when stepwise selection includes several variants in $S$, it suggests that either there are multiple causal variants in the region, or at least we need more than one variant to tag well a single causal variant in the region, e.g., because the causal variant is not genotyped. Let's write a function that does the stepwise selection and produces a plot of P-values after each step * showing the newly chosen SNP in red, * showing the preceding SNPs in $S$ in gray, * highlighting the true causals with blue. ```{r} stepwise.fwd <- function(y, X, p.thresh = 1e-4, plot.path = T, ring.i = NULL){ #Does stepwise forward selection where at most one variant is included in S at each iteration ## INPUT # y, quantitative trait for n individuals # X, genotypes for n individuals x p SNPs # p.thresh, P-value threshold that is used for deciding whether to include variant in S # plot.path, TRUE / FALSE, whether produces a plot of every step of the algorithm # ring.i, set of indexes of SNPs that should always be highlighted in plots by ringing ## OUTPUT # chosen.i, indexes of chosen variants in S, # chosen.p, the P-values at the iteration when each chosen variant was chosen, # next.p, the smallest P-value left outside S when finished. p = ncol(X) #number of SNPs col.new = "red" col.old = "gray" col.ring = "dodgerblue" pval = apply(X, 2, function(x){summary(lm(y ~ x))$coeff[2,4]}) #start from marginal P-values cols = rep("black", p) pchs = rep(1, p) chosen.i = c() #collect here the chosen SNPs chosen.p = c() #collect here the P-values from the iteration where the choice was made if(plot.path) { par(mfrow = c(1,3)) par(mar = c(1,2,4,0.1)) txt = "Initial state" plot(-log10(pval), xaxt = "n", xlab = "", ylab = "-log10P", main = txt, col = cols) if(!is.null(ring.i)) points(ring.i, -log10(pval[ring.i]), cex = 2, lwd = 1.4, col = col.ring) #highlight abline(h = -log10(p.thresh), lty = 2, col = "red") } ii = 0 #iteration index while( min(pval) < p.thresh ){ #continue as long as something is below the threshold ii = ii + 1 chosen.i[ii] = which.min(pval)[1] # add 1 SNP with min P-val to the chosen ones... chosen.p[ii] = pval[chosen.i[ii]] #... and store its P-value at this iteration. #test other SNPs except the already chosen ones -- and include the chosen ones as covariates tmp = apply(X[,-chosen.i], 2, function(x){summary(lm(y ~ x + X[,chosen.i]))$coeff[2,4]}) pval[-chosen.i] = tmp # we have P-values for other SNPs except already chosen ones in 'tmp' if(plot.path) { pval[chosen.i] = chosen.p # add causals with P-values temporarily for plotting cols[chosen.i] = col.old # earlier chosen remain gray -- cols[chosen.i[ii]] = col.new # -- the one newly chosen will be red pchs[chosen.i[ii]] = 19 #solid for chosen ones txt = paste0(ii,": chose ",chosen.i[ii]) plot(-log10(pval), xaxt = "n", xlab = "", ylab = "-log10P", main = txt, col = cols, pch = pchs) if(!is.null(ring.i)) points(ring.i,-log10(pval[ring.i]), cex = 2, lwd = 1.4, col = col.ring) abline(h = -log10(p.thresh), lty = 2, col = "red") } pval[chosen.i] = 1 #mark chose ones by P-value = 1 for the while-loop condition } return(list(chosen.i = chosen.i, chosen.p = chosen.p, next.p = min(pval))) } ``` We use a fairly liberal threshold of P=1e-4 here for two reasons: 1. We want the set $S$ to fully explain the association pattern in the region so that no peaks of P-values remain. 2. When we know that there is a robust GWS association in the region, then the variants in the region have a higher prior probability to be associated than random variants in the genome. There is no consensus value to use at stepwise conditioning. Sometimes it is run with the GWS-level of 5e-8, and sometimes continued until P > 0.05. Luckily, the path of the algorithm doesn't depend on the threshold used, so one can also cut the path later to be more stringent as needed for any particular application. ```{r,fig.height=3, fig.width=9} step.1 = stepwise.fwd(y.1, X, p.thresh = 1e-4, plot.path = T, ring.i = c.ind.1) ``` Here the stepwise selection worked well in the sense that it managed to pick exactly the three causal variants. Note also how the P-values of true causal variants 7 and 51 get smaller in steps 1 and 2, respectively, compared to the P-values from the previous iteration. This shows how conditioning on the other causal variants can reveal the true causal effects more clearly. **Example 7.5: APOE.2 data.** Stepwise forward selection doesn't always work quite so nicely: ```{r,fig.height=3, fig.width=9} set.seed(20) c.ind.2 = c(23,36,51) lambda = rep(0, p) lambda[c.ind.2] = 0.4 y.2 = X %*% lambda + rnorm(n, 0, sqrt(1-var(X %*% lambda)) ) res.2 = apply(X, 2, function(x){summary(lm(y.2 ~ x))$coeff[2,c(1,2,4)]}) step.2 = stepwise.fwd(y.2, X, p.thresh = 1e-4, plot.path = T, ring.i = c.ind.2) ``` Let's see how the chosen variants 20 and 66 are related to the causal ones 23, 36 and 51: ```{r} ind = c(step.2$chosen.i, c.ind.2) R[ind,ind] ``` It seems that none of the causals is very highly tagged by either of the chosen ones, but all causals show considerable correlation with the chosen ones. Let's check whether the chosen set explains the phenotype similarly to the causal variants. ```{r} lm.chosen = lm( y.2 ~ X[,step.2$chosen.i]) #chosen SNPs 20 and 66 summary(lm.chosen) lm.causal = lm( y.2 ~ X[,c.ind.2]) #causal SNPs 23, 36, 51 summary(lm.causal) ``` The chosen variants explain about 12.4% (R^2-adj.) of the trait variance whereas the three causals together explain about 15.0%. Note that we compare the adjusted-R^2 measure which attempts to take into account the number of parameteres in the models whereas multiple-R^2 tends to always favor model with more predictors. Another way to compare models with possibly different numbers of parameters is the [Bayesian Information Criterion](https://en.wikipedia.org/wiki/Bayesian_information_criterion) (BIC) that penalizes models for additional parameters. In short, a lower value of BIC suggests a better model. ```{r} cbind(chosen = BIC(lm.chosen), causal = BIC(lm.causal)) ``` We see that also BIC favors the true causal model over the one chosen by the stepwise selection (23 units in BIC is a clear difference). What if we put all 5 variants in the same model? ```{r} lm.joint = lm( y.2 ~ X[,ind]) #all 5 SNPs summary(lm.joint) ``` We see that, in this joint model, the two non-causals (20,66) picked by conditional analysis are not needed in the model once we have the causals there. (This model explains 15.1% while model with 3 causals already explained 15.0%.) We conclude that even though the stepwise selection picked a model that had a reasonable explanatory power, the model comparison shows that the chosen model has a statistically worse fit than the true causal model. Let's next move to recent work on fine-mapping methods to see how we can improve the search and quantification of causal variants. ```{r, eval = F, include = F} #Write files for FINEMAP res = apply(X, 2, function(x){summary(lm(y.2 ~ x))$coeff[2,1:2]}) data.name = "APOE_2" FM.path.1 = paste0("FINEMAP/FINEMAP_",data.name) FM.path.2 = paste0("FINEMAP_",data.name) write(R, ncol = ncol(R), file = paste0(FM.path.1,".ld")) write.table(cbind(rsid = 1:p, chromosome = rep(19,p), position = info[,2], allele1 = info[,3], allele2 = info[,4], maf = apply( cbind(info[,5], 1-info[,5]), 1, min), beta = res.2[1,], se = res.2[2,]), file = paste0(FM.path.1,".z"), quote = F, row.names = F) cat("z;ld;snp;config;cred;log;n_samples\n", file = paste0(FM.path.1,".master")) cat(paste0(FM.path.2,c(".z;",".ld;",".snp;",".config;",".cred;",".log;"),collapse = ""),n,"\n", sep ="",file = paste0(FM.path.1,".master"), append = T) ``` ### 7.5 Fine-mapping When we **fine-map** a GWAS region, our goal is to identify the set of causal variants $C$. Since high levels of LD among candidate variants often makes it impossible to give a definite answer based purely on statistical evidence, our answers will be probabilistic. In particular, we work with * probability for variant $l$ to be one of the causal variants in the region, $p_l = P(l\in C \,|\, \textrm{Data})$. This is also called **posterior inclusion probability (PIP)**. * probability for a set $S$ of variants to be the **causal set / causal configuration** of variants in the region, $p_S = P(S = C\, |\, \textrm{Data})$. * **credible sets** of causal variants which have the property that, with some fixed coverage, such as 95% probability, the credible set contains one of the causal variants or all of the causal variants, where the choice between these two definitions depends on the context. * probability distribution that the region contains $1,2,\ldots K$ causal variants up to some value of $K$, typically $\leq 20$. Note the difference between such quantities and the results from stepwise selection, which only included one set of variants, but did not try to estimate how likely each of those variants is to be a causal variant *given the LD patterns and other variants in the region*. Even if variant has very low marginal P-value, if it shares the signal with 10 other highly correlated variants, then its PIP will likely be $\leq 0.1$ and its credible set will contain all its 10 LD-friends, neither of which property can be seen from the output of stepwise selection. A recent review of fine-mapping methods is given by [Schaid et al.](https://www.nature.com/articles/s41576-018-0016-z) A common theme is a use of a Bayesian model that allows quantifying the probabilities of causality for each set and each variant, under the assumptions that all causal variants are included in the analysis and that there are no interaction effects between them. Slides 21-30 have some more details of the model behind a widely-used [FINEMAP](http://www.christianbenner.com/) software, written by Christian Benner at FIMM, Helsinki. Let's see what happens when FINEMAP is run on our two data sets APOE.1 and APOE.2 where we observed varying success with stepwise selection. #### FINEMAP on APOE.1 FINEMAP gives three output files: * **.snp** includes summary data (probabilities and causal effect sizes) for each SNP. * **.config** includes list of most probable sets of causal variants, with their probabilities and variance explained. * **.cred3** includes credible sets for 3 causal variants corresponding to the top-configuration in .config file. Let's read in FINEMAP's results for individual SNPs and print some columns. ```{r} fm.path = "https://www.mv.helsinki.fi/home/mjxpirin/GWAS_course/material/FINEMAP/FINEMAP_APOE_" fm.1.snp = read.table(paste0(fm.path,"1.snp"), as.is = T, header = T) fm.1.snp[1,] fm.1.snp[1:10,c(1,11,13:16)] ``` Here `prob` is the PIP, i.e., probability that this is one of the causal variants, `mean` is the estimated allelic causal effect $\lambda$, averaged over all possible causal sets and `mean_incl` is the estimate of $\lambda$ over those causal sets in which this SNP is included. The three simulated causals were 7, 37 and 51. The first is labelled as a definite causal with PIP ~ 100%. The variant 51 cannot be as perfectly singled out from the other correlated variant 52, but gets PIP of ~ 50%. The variant 37 is sharing its effect with several other correlated variants and has a PIP of only 11%. The causal effects in `mean_incl` seem estimating the true value of 0.4 given their uncertainty as measured by posterior SD. Let's plot the marginal GWAS P-value and FINEMAP's PIPs, and highlight the true causal variants. ```{r,fig.height = 4, fig.width = 9} par(mfrow = c(1,2)) par(mar = c(4.5,4.5,3,1)) plot(-log10(res.1[3,]), xlab = "SNP", ylab = "-log10 P", main = "MARGINAL GWAS APOE.1") points(c.ind.1, -log10(res.1[3,c.ind.1]), cex = 2, lwd = 1.4, col = "dodgerblue") plot(fm.1.snp$index, fm.1.snp$prob, xlab = "SNP", ylab = "prob of causality",main = "FINEMAP APOE.1") i = which(fm.1.snp$index %in% c.ind.1) points(fm.1.snp[i,"index"], fm.1.snp[i,"prob"], cex = 2, lwd = 1.4, col = "dodgerblue") ``` FINEMAPping has given a clearer picture here as there are less variants with interesting PIPs after fine-mapping than there are in the marginal P-values. Note how FINEMAP is most certain that variant 7 is a causal variant whereas the marginal P-value was smaller for variant 37 than variant 7. A likely reason for this difference is that 37 has more LD-friends with which it is forced to share some probability of being causal, whereas 7 is less correlated with others, and hence more clearly a causal itself. Such information cannot be extracted from the marginal P-values. The results in .snp file are SNP-wise summaries over all candidate causal configurations that FINEMAP has evaluated. By default .config file reports top 50,000 of these configs. Let's show a few top ones, with posterior probabilities and estimated proportion of variance explained by the config ($h^2$, heritability). ```{r} fm.1.config = read.table(paste0(fm.path,"1.config"), as.is = T, header = T) fm.1.config[1:10,c(2:4,8:9)] ``` Let's also see the credible sets given by FINEMAP. ```{r} fm.1.cred = read.table(paste0(fm.path,"1.cred3"), as.is = T, header = T) fm.1.cred[,-1] ``` For credible sets (with 3 causal variants), FINEMAP takes the top causal 3-SNP config (here 7,51,54), and then asks, for each variant in that set, which are the other candidates that could possibly replace that variant in this causal config. For example, column `cred_1` shows that the credible set $A_1 = \{7\}$ includes only SNP 7 as it already covers 0.99999 of the posterior probability of being in the causal set with 51 and 54. In other words, the sum of probabilities of any other variant replacing 7 in this causal set, that additionally contains SNPs 51 and 54, is < 0.00001. In $A_{2}=\{51,52\}$ the correct causal variant 51 is the the best candidate but there is also a second candidate, 52, having slightly smaller probability than 51 to be the second causal variant. For the third credible set $A_{3}$, the picture is less clear and we need 16 variants before we reach 95% credible set coverage. The true causal variant, 37, is the 3rd in the list. Credible set collects the most potential causal variants and is a place to search for known functional variants that could be responsible of the association signal. If the credible set is small enough, it also provides candidates to be taken further to experimental studies that could then unequivocally identify the causal variant among these candidates through biological function. The size of the credible set depends mainly on the GWAS sample size and the LD-structure, where large sample size and lower levels of LD both decrease the size of the credible sets. #### FINEMAP on APOE.2 On APOE.2 data causals were 23, 36 and 51. but due to LD, there were many other variants that had much lower P-values than the true causals. Conditional analysis chose 20 and 66, which was a clearly worse fit than the true causal set. Let's see FINEMAP's results. ```{r,fig.height = 4, fig.width = 9, echo = F} fm.2.snp = read.table(paste0(fm.path,"2.snp"), as.is = T, header = T) #fm.2.snp[1:10,c(1,11,17:20)] pchs = rep(1,p) par(mfrow = c(1,2)) par(mar = c(4.5,4.5,3,1)) plot(-log10(res.2[3,]), xlab = "SNP", ylab = "-log10 P", pch = pchs, main = "MARGINAL GWAS APOE.2") points(c.ind.2, -log10(res.2[3,c.ind.2]), cex = 2, lwd = 1.4, col = "dodgerblue") plot(fm.2.snp$index, fm.2.snp$prob, xlab = "SNP", ylab = "prob of causality", pch = pchs, main = "FINEMAP APOE.2", ylim = c(0,1)) i = which(fm.2.snp$index %in% c.ind.2) points(fm.2.snp[i,"index"], fm.2.snp[i,"prob"], cex = 2, lwd = 1.4, col = "dodgerblue") ``` We see that FINEMAP gives PIPs between 20% and 50% to the true causals, and that there are two non-causals getting even higher PIPs. Even if the true causals are not the absolute top PIPs, FINEMAP results are still much more robust for the true causal variants than the marginal P-values. Let's check the top 15 configs. ```{r, echo = F} fm.2.config = read.table(paste0(fm.path,"2.config"), as.is = T, header = T) fm.2.config[1:15,c(2:4,8:9)] ``` The top FINEMAP config is (24,42,51) whereas the true causal config (23,36,51) has rank 12. Let's compare statistical evidence between the top FINEMAP config, the true causal config and the conditional analysis config (20,66). ```{r} lm.fm = lm(y.2 ~ X[,c(24,42,51)]) lm.causal = lm(y.2 ~ X[,c(23,36,51)]) lm.step = lm(y.2 ~ X[,c(20,66)]) data.frame(finemap = BIC(lm.fm), true_causal = BIC(lm.causal), stepwise = BIC(lm.step)) ``` The top config from FINEMAP happens to be slightly better description of the data than the true causal config, which is a result of statistical sampling effects in this finite sample data set. We also see that the conditional analysis has failed badly here compared to FINEMAP in terms of explaining the data (26 units difference in BIC). We conclude that FINEMAP has made an efficient search throught the space of configurations and identified a config with more statistical evidence than the true causal config. Given the data available, we can't expect to get better fine-mapping results from statistical evidence alone. #### 7.5.1 Fine-mapping assuming one causal variant Sometimes we do not have access to the original GWAS data and therefore we cannot accurately estimate LD that corresponds to the GWAS summary statistics. This means that we can't do full fine-mapping with FINEMAP or other similar methods that require LD information. In these cases, it is still possible to compute, for each variant in the region, the posterior probability of causality under the assumption that **there is at most one causal variant in the region and it is included in the study**. The approach was introduced by [Maller et al. 2012](https://www.nature.com/articles/ng.2435), together with the idea of credible sets. Let's mark by $H_l, l=1,\ldots,p$ the hypothesis that says $l$ is the causal variant and let $H_0$ be the null hypothesis that says there is no causal variant in the region. Mark by $H_C = H_1 \cup \ldots \cup H_p$ the hypothesis that says that there is exactly one causal variant in the region. We assume that, *a priori*, each variant is equally likely to be the causal variant, i.e., $P(H_l\,|\,H_C) = \frac{1}{p}$. It follows that the regional Bayes factor is $$\textrm{BF}_{\textrm{reg}} = \frac{P(\textrm{Data}\, | \, H_C)}{P(\textrm{Data}\, | \, H_0)}= \frac{\sum_{l=1}^p \frac{1}{p} P(\textrm{Data}\, | \, H_l)}{P(\textrm{Data}\, | \, H_0)}= \frac{1}{p} \sum_{l=1}^p \textrm{BF}_l,$$ where $\textrm{BF}_l$ is the single variant Bayes factor that does not depend on other variants or the LD structure of the region. In practice, it can be computed as Approximate Bayes Factor (ABF) from GWAS summary statistics, as explained in Section 4 of the course material. The posterior probability for variant $l$ is then $$P(H_l \, | \, \textrm{Data}, H_C) = \frac{P(H_l\, | \, H_C) P(\textrm{Data} \, | \, H_l)} {P(\textrm{Data} \, | \, H_C)} = \frac{1}{p} \frac{\textrm{BF}_l}{\textrm{BF}_{\textrm{reg}}} =\frac{\textrm{BF}_l}{\sum_{k=1}^p \textrm{BF}_k} \propto \textrm{BF}_l.$$ Thus, the posterior probability that SNP $l$ is the single causal variant in the region, under the assumption that there is exactly one causal variant, is proportional to the single-SNP Bayes factor of SNP $l$. In practice, we get the posterior probabilities by computing ABF for each variant and then normalizing them to sum to one. This derivation is independent of LD in the region because, under the assumption of only a single causal SNP, the marginal effect at the causal SNP is the causal effect itself, independent of the LD with other variants. With significantly associated but also highly correlated SNPs, we are left with the fact that any of them could be the causal. This will manifest through high single-SNP and regional BFs, but with the posterior distributed fairly evenly across multiple SNPs. We have assumed that the causal SNP is included in the study. When this does not hold, the above method is still applicable if there is a good surrogate SNP for the causal effect. In the presence of multiple causal SNPs, this approach is no longer optimal. It will tend to pick out the SNP with the best marginal effect, which may or may not be one of the causal SNPs, as we have seen before. **Example 7.6: APOE.3 data.** Let's make yet another phenotype data set for the *APOE* region genotype data that we currently have in `X` matrix. Now we make one causal variant, SNP 42. ```{r, fig.width = 9, fig.height = 4} set.seed(6) c.ind = c(42) lambda = rep(0, p) lambda[c.ind] = 0.3 y.3 = X %*% lambda + rnorm(n, 0, sqrt(1-var(X %*% lambda)) ) res.3 = apply(X, 2, function(x){summary(lm(y.3 ~ x))$coeff[2,c(1,2,4)]}) #runs GWAS in the region tau = 0.2 #For computing ABFs b.est = res.3[1,] se = res.3[2,] log.abf = dnorm(b.est, 0 , sqrt(tau^2 + se^2), log = T) - dnorm(b.est, 0, se, log = T) abf = exp(log.abf) posterior = abf/sum(abf) #Let's see which variants have probability > 1% ind = which(posterior > 0.01) ind = ind[order(posterior[ind], decreasing = T)] #from largest to smallest data.frame(SNP = ind, prob = posterior[ind], Pvalue = res.3[3,ind]) #Let's plot the original -log10 P-values and posterior probabilities par(mfrow = c(1,2)) cols = rep("gray", p) cols[ind] = "indianred2" plot(-log10(res.3[3,]), xlab = "SNP", ylab = "-log10(P-value)", main = "APOE.3", cex = 0.7, col = cols, pch = 19) plot(posterior, xlab = "SNP", ylab = "posterior probability", main = "", cex = 0.7, col = cols, pch = 19) ``` We see that the 3 smallest P-values also have the largest posterior probabilities of being causal. The SNP 42 (78%) has about 4 times higher probability than SNP 36 (18%) and 50 times than SNP 30 (1.6%) of being the causal variant in the region (when we assume that there is only one causal variant). To determine the 95% credible set, we would include the SNPs with the highest posteriors until their cumulative sum is > 95%. Here this set would contain only SNPs 42 and 36 as their sum is already 78% + 18% = 96%. If we would have several conditionally independent association signals in the region, as returned by the stepwise forward search, then we could fine-map each of them separately. For this, we would compute the association statistics (beta and SE) by including the SNPs representing the other signals in the region as covariates, and then carry out the fine-mapping as above by assuming that there is only one causal variant left after conditioning on the other signals in the region (in Exercise 5.2 we do that for one signal conditioned on two others).