CT-LDSC

Cross Trait Linkage Disequilibrium Score Regression (CT-LDSC)¹ is an extension of Linkage Disequilibrium Score Regression (LDSC)² that estimates genetic correlation.

Derivation for Quantitative Traits

Data Generating Model

Recall that LDSC assumes a linear data-generating model. CT-LDSC assumes two linear data-generating models, one for each of two traits:

\[ \begin{align} y_1&= Y\beta_1 + \delta \label{dg1}\\ y_2&= Z\beta_2 + \epsilon \label{dg2} \end{align} \]

where:

There are \(M \gg 0\) genetic variants.
There are \(N_1 \gg 0\) individuals in the trait 1 GWAS and \(N_2 \gg 0\) individuals in the trait 2 GWAS.
There are \(N_s\) individuals included in both GWAS. Without loss of generality, assume these \(N_s\) individuals are listed first in the lists of participants in both studies.
\(y_1\in \mathbb{R}^{N_1}\) and \(y_2 \in \mathbb{R}^{N_2}\) are the phenotype vectors for the trait 1 and trait 2 GWAS, respectively.
\(Y\in\mathbb{R}^{N_1\times M}\) and \(Z\in\mathbb{R}^{N_2\times M}\) are the genotype matrices from the two GWAS, normalized to have columns with sample mean 0 and variance 1. \(Y_j\in\mathbb{R}^{N_1}\) and \(Z_j \in \mathbb{R}^{N_2}\) denote the \(j\)th columns of the two matrices.
\(\beta_1,\beta_2\in\mathbb{R}^M\) are the vectors of true per-variant genetic effect sizes for the two traits.
\(Y\beta_1\in\mathbb{R}^{N_1}\) and \(Z\beta_2\in\mathbb{R}^{N_2}\) are thus the vectors of per-individual genetic effects in the two GWAS.
\(\epsilon\in\mathbb{R}^{N_1}\) and \(\delta\in\mathbb{R}^{N_2}\) are the vectors of non-genetic effects in the two GWAS.

We model \(Y,Z,\beta_1,\beta_2, \delta,\epsilon\) as random variables with the following properties:

\[ \begin{align} \mathbb{Var(\delta)} &= (1-h_1^2)I & \text{ for some }h_1>0\\ \mathbb{Var(\epsilon)} &= (1-h_2^2)I& \text{ for some }h_2>0\\ \mathbb{Cov}(\delta_j,\epsilon_k) &= \begin{cases} \rho_e &\text{ if } j=k \le N_s \\ 0 & \text{ else } \end{cases} & \text{ for some }\rho_e>0 \label{cov_delta_epsilon}\\ \mathbb{Var}(\beta_1) &= \frac{1}{M} h_1^2 I \label{var_beta} \\ \mathbb{Var}(\beta_2) &= \frac{1}{M} h_2^2 I \label{var_gamma} \\ \mathbb{Cov}(\beta_1, \beta_2)&= \frac{1}{M}\rho_gI & \text{ for some }\rho_g>0 \label{cov_beta_gamma}\\ \mathbb{E}(\delta)&=0\label{e_delta}\\ \mathbb{E}(\beta_1)&=0 \label{e_beta} \\ \mathbb{E} Y &=0 \label{e_Y} \\ \mathbb{E}(\epsilon)&=0 \label{e_epsilon} \\ \mathbb{E}(\beta_2)&=0 \label{e_gamma} \\ \mathbb{E} (Z) &=0 \label{e_Z} \\ \mathbb{E}(Y_i^2) &=\mathbb{E}(Z_i^2)=1 &\text{ for all }i \label{y_z_var} \end{align} \]

Assumption (\(\ref{cov_delta_epsilon}\)) specifies that error is uncorrelated between GWAS, except for individuals who are participants both studies.
Assumptions (\(\ref{var_beta}\)), (\(\ref{var_gamma}\)), and (\(\ref{cov_beta_gamma}\)) specify an isotropically polygenic architecture, analogous to the one used in LDSC.
Assumptions (\(\ref{e_delta}\)), (\(\ref{e_beta}\)), (\(\ref{e_Y}\)), (\(\ref{e_Z}\)), and (\(\ref{y_z_var}\)) specify that the data has been pre-normalized. Here (\(\ref{e_Y}\)), (\(\ref{e_Z}\)), and (\(\ref{y_z_var}\)) elide the distinction between sample and population level normalization, which should only introduce minor error given the large sample size of a typical GWAS.

Furthermore, we assume the following relationships between the random variables, which are fairly standard in linear data-generating models:

The rows of \(Y\) and \(Z\) are mutually independent, except in the case that a row of \(Y\) and \(Z\) refer to the same individual (one of the \(N_s\) individuals present in both GWAS).
The rows of \(Y\) and \(Z\) are identically distributed. That is, the two GWAS samples are drawn independently from the sample population.
\(Y, \beta_1, \delta\) are all mutually independent, as are \(Z,\beta_2, \epsilon\).
\(\beta_1\) is independent of \(\epsilon\) and \(Z\), and \(\beta_2\) is independent of \(\delta\) and \(Y\).

Note that this implies that the true genetic effects have mean zero:

\[ \begin{align} \mathbb{E}( Y\beta_1) &= \mathbb{E}(Y) \mathbb{E}(\beta_1) &\text{ by independence}\\ &=0 \end{align} \]

and similarly, \(\mathbb{E}(Z\beta_2)=0\).

Define the following quantities related to Linkage Disequilibrium (LD):

The LD between SNP \(j\) and SNP \(k\) is denoted by \(r_{jk}:=\mathbb{E}(Y_{i,j}Y_{i,k})=\mathbb{E}(Z_{q,j}Z_{qk})\), (which does not depend on the individuals \(i\) and \(q\) by our assumption that the rows of \(Y\) and \(Z\) are iid).
The LD score of a SNP \(j\) is defined to be \(l_j:= \sum_k r_{jk}^2\). It quantifies the magnitude of the dependence between \(j\) and other SNPs.

Denote by \(\hat\beta_1,\hat\beta_2\in\mathbb{R}^M\) the marginal GWAS regression coefficients from GWAS 1 and 2. Note the critical distinction between \(\beta_1,\beta_2\), which are underlying causal, joint effects of genetic variants, and \(\hat\beta_1,\hat\beta_2\), which are the outputs of per-variant univariate marginal regressions.

Genetic Covariance

First, let us compute the genetic covariances between the two phenotypes. Let \(X\in\mathbb{R}^M\) denote the genotype of an arbitrary individual. By definition, genetic covariance is:

\[ \begin{align} &\mathbb{Cov}\left( \sum_j X_j \beta_{j,1}, \sum_k X_k \beta_{k,2} \right)\\ &= \mathbb{E}\left( \left(\sum_j X_j \beta_{j,1}\right)\left(\sum_k X_k \beta_{k,2} \right) \right)& \text{mean 0}\\ &=\sum_{j,k} \mathbb{E} \left( X_j X_k \beta_{j,1} \beta_{k,2} \right) & \text{linearity}\\ &=\sum_{j,k}\mathbb{E} (X_jX_k) \mathbb{E}(\beta_{j,1} \beta_{k,2}) & \text{independence}\\ &=\sum_{j}\mathbb{E} (X_j^2) \mathbb{E}(\beta_{j,1} \beta_{j,2}) & \text{by (\ref{cov_beta_gamma})}\\ &=\sum_j \mathbb{E} (\beta_{j,1} \beta_{j,2}) & \text{ by (\ref{y_z_var})}\\ &=\rho_g. & \text{by (\ref{cov_beta_gamma})} \end{align} \]

Genetic Correlation

Note that the model (\(\ref{dg1},\ref{dg2}\)) is an extension of the model used in LDSC. By the derivation of LDSC, we have \(\mathbb{Var}(\sum_j X_j \beta_{j,1})=h_1^2\) and \(\mathbb{Var}(\sum_j X_j \beta_{j,2})=h_2^2\).

Since we have assumed phenotypic variance is standardized to 1 (see \((\ref{e_Y})\), \((\ref{e_Z})\), and \((\ref{y_z_var})\)), it follows that the genetic correlation of the two traits can be computed as their genetic covariance divided by the square root of the product of their heritabilities:

\[ \frac{\rho_g}{h_1h_2}. \]

Regression equation

Using the same logic as in derivation of LDSC, we approximate the Wald \(\chi^2\) statistics of a variant \(j\) in the two GWAS:

\[ \begin{align} \chi_{j,1}^2 &\approx N \hat{\beta_1}_{j,1}^2=\frac{(y_1^T Y_{j})^2}{N_1}\\ \chi_{j,2}^2 &\approx N \hat{\beta_1}_{j,1}^2=\frac{(y_1^T Z_{j})^2}{N_2}. \end{align} \]

It follows that the corresponding z-statistics are:

\[ \begin{align} z_{j,1} &\approx \frac{(y_1^T Y_{j})}{\sqrt{N_1}} \label{d1}\\ z_{j,2} &\approx \frac{y_1^T Z_{j}}{\sqrt{N_2}} \label{d2}. \end{align} \]

Recall that in LDSC, the regression dependent variable was the Wald \(\chi^2\) statistic. In CT-LDSC, the regression dependent variable is the product of the \(z\)-statistics from the two GWAS. To derive the regression equation, we must estimate \(\mathbb{E}(z_{j,1}z_{j,2})\). By the Tower Law of Expectation (see Grimmett and Stirzaker pg.336³),

\[ \begin{align} \mathbb{E}(z_{j,1}z_{j,2})= \mathbb{E}\left( \mathbb{E}(z_{j,1}z_{j,2}|Z,Y) \right). \\ \end{align} \]

Evaluating the inner expectation,

\[ \begin{align} &\mathbb{E}(z_{j,1}z_{j,2}|Z,Y)\\ &\approx\frac{1}{\sqrt{N_1N_2}} \mathbb{E}\left(Y_j^T y_1 y_1^T Z_j|Z,Y \right) & \text{ by (\ref{d1}, \ref{d2})}\\ &=\frac{1}{\sqrt{N_1N_2}}Y^T_j \mathbb{E}\left( (Y\beta_1+\delta)(Z\beta_2+\epsilon)^T |Z,Y\right)Z_j & \text{ by (\ref{dg1},\ref{dg2})}\\ &=\frac{1}{\sqrt{N_1N_2}}Y^T_j \left(Y \mathbb{E}(\beta_1 \beta_2^T)Z^T + \mathbb{E}(\delta \beta_2^T) Z^T + Y \mathbb{E}(\beta_1 \epsilon^T) + \mathbb{E}(\epsilon \delta^T) \right)Z_j & \text{linearity, independence}\\ &=\frac{1}{\sqrt{N_1N_2}}Y^T_j\left(\frac{\rho_g}{M}YZ^T + \rho_e Q\right)Z_j & \text{ by (\ref{cov_beta_gamma},\ref{cov_delta_epsilon})}\\ &=\frac{1}{\sqrt{N_1N_2}}\left(\frac{\rho_g}{M}Y_j^TYZ^TZ_j + \rho_e Y_j^T QZ_j\right) \label{cond_exp} \end{align} \]

where \(Q\in\mathbb{R}^{N_1\times N_2}\) and

\[ \begin{align} Q_{i,j}&= \begin{cases} 1 & \text{ if }i=j \le N_s\\ 0 & \text{else} \end{cases}. \end{align} \]

Our goal is to derive an expression for the unconditional expectation of \(z_{j,1}z_{j,2}\). To achieve this, we take the expectation of both terms in \((\ref{cond_exp})\). First, we have

\[ \begin{align} &\mathbb{E}(Y_j^TYZ^TZ_j)\\ &=\mathbb{E} Y_j^T \left( \sum_{k=1}^M Y_{i,k} Z_{q,k} \right)_{(i,q)} Z_j \\ &=\mathbb{E} \sum_{i=1}^{N_1}\sum_{q=1}^{N_2}\sum_{k=1}^M\left( Y_{i,j}Y_{i,k} Z_{q,j} Z_{q,k} \right) & \text{ def of matrix product} \end{align} \]

There are \(N_1N_2\) possible pairs of the indexes \(i\) and \(q\). For each such pair, there are two cases:

Case 1: \(i\) and \(q\) refer to the same individual. There are \(N_s\) such index pairs. In this case we have

\[ \begin{align} &\sum_{k=1}^M \mathbb{E}(Y_{i,j} Y_{i,k} Z_{q,j} Z_{q,k} )\\ &=\sum_{k=1}^M \mathbb{E}(Y_{i,j}^2 Y_{i,k}^2)\\ &\approx \sum_{k=1}^M (1 + 2r_{j,k}^2) \label{isserlis_approx_line}\\ &=M + 2\sum_{k=1}^M r_{j,k}^2\\ &= M + 2l_j \end{align} \]

Where we have approximated the random variables as having a normal distribution, and used Isserlis's Theorem. Recall that we also used Isserlis's theorem at a similar point in the derivation of LDSC.

Case 2: \(i\) and \(q\) refer to different individuals. There are \(N_1N_2-N_s\) such pairs.

\[ \begin{align} &\sum_{k=1}^M \mathbb{E}(Y_{i,j} Y_{i,k} Z_{q,j} Z_{q,k} )\\ &=\sum_{k=1}^M \mathbb{E}(Y_{i,j} Y_{i,k}) \mathbb{E}(Z_{q,j} Z_{q,k} ) & &\text{ independence} \\ &=\sum_{k=1}^M r_{jk}^2\\ &= l_{j} \end{align} \]

Combining the two cases yields

\[ \begin{align} &\mathbb{E}(Y_j^TYZ^TZ_j)\\ &=N_s(M+2l_j) + (N_1N_2-N_s)l_j\\ &=MN_s + (N_1N_2+N_s)l_j\\ &\approx MN_s + N_1N_2 l_j. \label{exp_prod_1} \end{align} \]

Returning to the second term in \((\ref{cond_exp})\), we have that

\[ \begin{align} \mathbb{E}(Y_j^TQZ_j)=N_s \label{exp_prod_2} \end{align} \]

by the assumption of independence of the genotypes of distinct individuals.

Substituting \((\ref{exp_prod_1})\) and \((\ref{exp_prod_2})\) into \((\ref{cond_exp})\) yields:

\[ \begin{align} &\mathbb{E}(z_{j,1}z_{j,2})\\ &\approx \frac{1}{\sqrt{N_1N_2}}\left( \frac{\rho_g}{M} (MN_s + N_1N_2 l_j) +\rho_e N_s \right)\\ &= \frac{\rho_g\sqrt{N_1N_2}}{M}l_j + \frac{N_s}{\sqrt{N_1N_2} } (\rho_g + \rho_e ).\label{ct_ldsc_eqn} \end{align} \]

\((\ref{ct_ldsc_eqn})\) is the key regression equation in CT-LDSC.

Derivation for Binary Traits

Above, we assumed quantitative traits, but most disease traits are binary. Here, we modify the CT-LDSC derivation to cover the case of binary traits and ascertained case-control GWAS.

Liability-Scale Data Generating Model

We model binary traits via a liability threshold .

We have

\[ \begin{align} \psi_1&= Y\beta_1 + \delta \label{bin_dg1}\\ \psi_2&= Z\beta_2 + \epsilon \label{bin_dg2}\\ y_1 &= 1_{\psi_1>\tau_1}\\ y_2 &= 1_{\psi_2>\tau_2}\\ \end{align} \]

where

There are \(M \gg 0\) genetic variants.
There are \(N_1 \gg 0\) individuals in the trait 1 GWAS and \(N_2 \gg 0\) individuals in the trait 2 GWAS.
There are \(N_s\) individuals included in both GWAS. Without loss of generality, assume these \(N_s\) individuals are listed first in the lists of participants in both studies.
\(y_1\in \{0,1\}^{N_1}\) and \(y_2 \in \{0,1\}^{N_2}\) are the binary phenotype vectors for the trait 1 and trait 2 GWAS, respectively.
\(G_1\in\{0,1\}^{N_1\times M}\) and \(G_2\in\{0,1\}^{N_2\times M}\) are the raw genotype matrices of the two GWAS. Like the original paper, we work with haploid genomes for simplicity. The authors argue that results can easily be extended to diploid genomes.
\(p_j\in\mathbb{R}\) is the population-level probability that an arbitrary individual has variant \(j\).
\(Y\in\mathbb{R}^{N_1\times M}\) and \(Z\in\mathbb{R}^{N_2\times M}\) are the population-standardized genotype matrices from the two GWAS. Thus \(Y_{i,j}=\frac{G_{1,i,j}-p_j}{\sqrt{p_j(1-p_j)}}\) and \(Z_{i,j}=\frac{G_{2,i,j}-p_j}{\sqrt{p_j(1-p_j)}}\). Note that because of case-control ascertainment, the columns of \(Y\) and \(Z\) may have sample means far from 0 and variances far from 1.
\(Y_j\in\mathbb{R}^{N_1}\) and \(Z_j \in \mathbb{R}^{N_2}\) denote the \(j\)th columns of the two matrices.
\(\beta_1,\beta_2\in\mathbb{R}^M\) are the vectors of true per-variant genetic effect sizes for the two traits.
\(\epsilon\in\mathbb{R}^{N_1}\) and \(\delta\in\mathbb{R}^{N_2}\) are the vectors of non-genetic effects in the two GWAS.
\(\tau_1,\tau_2\in\mathbb{R}\) are the liability thresholds for traits 1 and 2 respectively.
\(\psi_1\in\mathbb{R}^{N_1}\) is the latent liability score for trait 1 for participants in GWAS 1. \(\psi_2\in\mathbb{R}^{N_2}\) is the latent liability score for trait 2 for participants in GWAS 2.

Let \(K_1 \in (0,1)\) denote the proportion of individuals with trait 1 in the general population, and let \(K_2 \in (0,1)\) denote the proportion of individuals with trait 2 in the general population. Because of case-control ascertainment these proportionals may not hold in the GWAS. We let \(P_1\in (0,1)\) denote the proportion of GWAS 1 participants with trait 1, and \(P_2\) denote the proportion of GWAS 2 participants with trait 2.

We model \(G_1,G_2,\beta_1,\beta_2,\delta,\epsilon\) as random variables with the following properties

\[ \begin{align} \mathbb{Var(\delta)} &= (1-h_1^2)I & \text{ for some }h_1>0\\ \mathbb{Var(\epsilon)} &= (1-h_2^2)I& \text{ for some }h_2>0\\ \mathbb{Cov}(\delta_j,\epsilon_k) &= \begin{cases} \rho_e &\text{ if } j=k \le N_s \\ 0 & \text{ else } \end{cases} & \text{ for some }\rho_e>0 \label{bin_cov_delta_epsilon}\\ \mathbb{Var}(\beta_1) &= \frac{1}{M} h_1^2 I \label{bin_var_beta} \\ \mathbb{Var}(\beta_2) &= \frac{1}{M} h_2^2 I \label{bin_var_gamma} \\ \mathbb{Cov}(\beta_1, \beta_2)&= \frac{1}{M}\rho_gI & \text{ for some }\rho_g>0 \label{bin_cov_beta_gamma}\\ \mathbb{E}(\delta)&=0\label{bin_e_delta}\\ \mathbb{E}(\beta_1)&=0 \label{bin_e_beta} \\ \mathbb{E}(\epsilon)&=0 \label{bin_e_epsilon} \\ \mathbb{E}(\beta_2)&=0 \label{bin_e_gamma} \\ \end{align} \]

\(h_1,h_2\) are the liability-scale heritabilities of trait 1 and trait 2, while \(\rho_g\) is the liability-scale genetic covariance.

We assume the following relationships between the random variables:

The rows of \(G_1\) and \(G_2\) are independent, except when they refer to the same individual (one of the \(N_s\) individuals present in both GWAS).
\(Y, \beta_1, \delta\) are all mutually independent, as are \(Z,\beta_2, \epsilon\).
\(\beta_1\) is independent of \(\epsilon\) and \(Z\), and \(\beta_2\) is independent of \(\delta\) and \(Y\).

For convenience, also define \(g_p\in\{0,1\}^{M}\) to be the genotype of an individual selected at random from the general population. Define \(x_p\in\mathbb{R}^{M}\) by \(x_{p,j}=\frac{g_{p,j}-p_j}{\sqrt{p_j(1-p_j)}}\) so that \(x_p\) is the corresponding standardized genotype vector. Note that we now have three genotype distributions:

The distribution of genotypes in the general population (\(g_p,x_p\)).
The distribution of genotypes in GWAS 1 (\(G_1,Y\)).
The distribution of genotypes in GWAS 2 (\(G_2,Z\)).

Define the following Linkage Disequilibrium (LD) quantities:

The LD between SNP \(j\) and SNP \(k\) is denoted by \(r_{jk}:=\mathbb{E}(x_{p,j}x_{p,k})\). Note that we use the population-level distribution here.
The LD score of a SNP \(j\) is defined to be \(l_j:= \sum_k r_{jk}^2\). It quantifies the magnitude of the dependence between \(j\) and other SNPs.

Marginal Regression Coefficients

Analogously to the above, let \(\hat\beta_1,\hat\beta_2\in\mathbb{R}^M\) denote the marginal sample GWAS regression coefficients from GWAS 1 and 2. Thus \(\hat\beta_{1,j}\) is the outcome of regressing \(y_1\) on \(Y_i\).
We will assume that ordinary linear regression is used.
We also introduce population level marginal GWAS regression coefficients \(\overline\beta_1,\overline\beta_2\in\mathbb{R}^{M}\). \(\hat\beta_1\to\overline\beta_1\) as \(N_1\to\infty\) and \(\hat\beta_2\to\overline\beta_2\) as \(N_2\to\infty\).
Recall that the various

Note on derivation

The derivation in the supplementary material to the original paper¹ contains several typos and implicit approximations. In the version above, I have corrected typos and make approximations explicit.

Brendan Bulik-Sullivan, Hilary K Finucane, Verneri Anttila, Alexander Gusev, Felix R Day, Po-Ru Loh, ReproGen Consortium, Psychiatric Genomics Consortium, Genetic Consortium for Anorexia Nervosa of the Wellcome Trust Case Control Consortium 3, Laramie Duncan, and others. An atlas of genetic correlations across human diseases and traits. Nature Genetics, 47(11):1236–1241, 2015. URL: https://pmc.ncbi.nlm.nih.gov/articles/PMC4797329/. ↩↩
Brendan K Bulik-Sullivan, Po-Ru Loh, Hilary K Finucane, Stephan Ripke, Jian Yang, Nick Patterson, Mark J Daly, Alkes L Price, and Benjamin M Neale. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nature Genetics, 47(3):291–295, 2015. URL: https://pmc.ncbi.nlm.nih.gov/articles/PMC4495769/. ↩
Geoffrey Grimmett and David Stirzaker. Probability and Random Processes: Third Edition. Oxford university press, 2001. URL: https://www.amazon.ca/Probability-Random-Processes-3th-third/dp/B006QYQVTS. ↩