. 2020 Apr 16;106(5):646–658. doi: 10.1016/j.ajhg.2020.03.011

Pleiotropy-Based Decomposition of Genetic Risk Scores: Association and Interaction Analysis for Type 2 Diabetes and CAD

Daniel I Chasman ^1,^2,^3,^4,^5,^∗, Franco Giulianini ¹, Olga V Demler ^1,³, Miriam S Udler ^3,^4,⁶

PMCID: PMC7212269 PMID: 32302534

Abstract

Genetic risk for a disease in the population may be represented as a genetic risk score (GRS) constructed as the sum of inherited risk alleles, weighted by allelic effects established in an independent population. While this formulation captures overall genetic risk, it typically does not address risk due to specific biological mechanisms or pathways that may nevertheless be important for interpretation or treatment response. Here, a GRS for disease is resolved into independent or nearly independent components pertaining to biological mechanisms inferred from pleiotropic relationships. The component GRSs’ weights are derived from the singular value decomposition (SVD) of the matrix of appropriately scaled genetic effects, i.e., beta coefficients, of the disease variants across a panel of the disease-related phenotypes. The SVD-based formalism also associates combinations of disease-related phenotypes with inferred disease pathways. Applied to incident type 2 diabetes (T2D) in the Women’s Genome Health Study (N = 23,294), component GRSs discriminate glycemic control and lipid-based genetic risk, while revealing significant interactions between specific components and BMI or physical activity, the latter not observed with a GRS for overall T2D genetic liability. Applied to coronary artery disease (CAD) in both the WGHS and in JUPITER (N = 8,749), a randomized trial of rosuvastatin for primary prevention of CVD, component GRSs discriminate genetic risk associated with LDL-C from risk associated with reciprocal genetic effects on triglycerides and HDL-C. They also inform the pharmacogenetics of statin treatment by demonstrating that benefit from rosuvastatin is as strongly related to genetic risk from triglycerides and HDL-C as from LDL-C.

Keywords: genetic risk score, polygenic risk score, type 2 diabetes, coronary artery disease, singular value decomposition, interaction, statin

Introduction

The effects of pre-specified genetic variants on a phenotype, e.g., incident T2D or CAD, may be added to yield a single number, the genetic risk score (GRS), which approximates the inherited genetic liability for an individual. The contribution of each variant in the GRS derives from the product of a weight corresponding to the effect estimate, typically a regression beta coefficient from an independent sample, together with a suitable representation of the inherited genotype. The accuracy of the weights is critical to the utility of the GRS, and it has been shown that as sample sizes for GWASs continue to increase, the effect estimates asymptotically reach their true effects such that a GRS may ultimately represent a substantially complete manifestation of inherited common genetic risk.¹ It is now common to anticipate that such GRS representations of inherited genetic risk may join conventional clinical measures in applications for personalized medicine, especially for diseases with appreciable heritability.

By combining effects of all relevant genetic variants under the assumption of additivity, the GRS approach implicitly also amalgamates all genetically influenced biologicial functions contributing to the relevant phenotype. The GRS does not distinguish among separate genetic liabilities invested in each contributing biological pathway, although such distinctions may be important. For example, since different therapeutic strategies for the same condition may target different pathways, knowing an individual’s pathway-specific genetic risk burden may inform the choice of treatment option. Indeed, testing for the interaction between a clinical exposure, e.g., treatment or a lifestyle choice, and a GRS implicitly assumes that the ratio of effects in the exposed compared to unexposed strata will be the same for all of the constituent SNPs.² However, if the effects of an exposure are mediated through a specific pathway, this ratio may differ for SNPs according to their biological roles. As a consequence, the ability to detect statistical interactions of an expsoure with a single GRS will be degraded compared with a GRS that represents only a relevant pathway.

While pathways may be broadly defined³^,⁴^,⁵^,⁶^,⁷ as combinations of genes assigned to shared function, an alternative approach might leverage pleiotropy in SNP associations across phenotypes. Such an approach would depend solely on relative SNP effects on phenotypes and obviate the need to assign SNPs to genes or to infer knowledge of gene function.⁸^,⁹^,¹⁰ For common variants with small effects such as those from GWASs, pleiotropy may be viewed as the consequence of modest pertubations of various underlying pathways. Therefore, variants with similar profiles of phenotypic effects may be assigned to the same pathway. This perspective recently suggested five or six groupings of SNPs from GWASs leading to etiologically distinct forms of T2D.¹¹^,¹²

Here, we present an approach that leverages pleiotropy to decompose GRSs for disease via the singular value decomposition (SVD) of a matrix of appropriately scaled SNP effects, i.e., GWAS beta coefficients, for disease-related phenotypes. We apply the approach to incident T2D and CAD. We show that the weights of the SVD decomposition delineate combinations of the disease-related phenotypes, revealing recognizable disease etiologies and their associated gene targets. Moreover, the GRS components reveal or strengthen statistical evidence for interactions with clinical exposures for incident type 2 diabetes (T2D) and with randomized statin allocation for incident cardiovascular disease (CVD), compared to interactions assessed using the overall disease GRSs.

Subjects and Methods

Cohorts

The Women’s Genome Health Study (WGHS)

The WGHS is a prospective cohort for genome-wide genetic analysis derived from the Women’s Health Study (WHS) that was a balanced 2×2 factorial placebo-controled trial of aspirin (100 mg every other day) and vitamin E (600 IU every other day) over 10 years for primary prevention of CVD and cancer among female health care professionals, 45 years old and free of known prevalent endpoint conditions at baseline.¹³ Since the end of the trial period in 2004, follow-up has continued in observational mode. Approximately 71% WHS participants provided a baseline blood sample within which lipid fractions were measured and DNA was extracted for genotyping. Demographic and lifestyle data, including physical activity measures, were collected by questionnaire. Follow-up was by annual questionnaire. Incident coronary heart disease (CHD) was verified by physician review of medical records according to WHO standards and defined as a composite of myocardial infarction, coronary revascularization, and coronary heart disease death. This definition is essentially equivalent to the definition of coronary artery disease (CAD) that was used in the GWAS to discovery genetic variants incorporated into the risk score analysis, although the latter may have also included cases of angina for some cohorts contributing to the GWAS.¹⁴^,¹⁵ Incident type 2 diabetes was ascertained by self-report, followed by physician phone interview or supplemental questionnaire.¹⁶ For analysis of incident diabetes, 2.4% of study participants with diabetes at baseline were excluded.

JUPITER Study

JUPITER (justification for the use of statins in prevention, an intervention trial evaluating rosuvastatin) was a randomized, placebo-controlled trial of rosuvastatin (20 mg/d) for primary prevention of incident CVD among individuals with elevated C-reactive protein (≥2 mg/L) and moderate to low LDL-C (<130 mg/dl) over a median follow-up of 1.9 years.¹⁷ Within the baseline and 12 month follow-up blood samples, plasma lipids were determined by standard clinical assays. The primary endpoint of the trial was a composite nonfatal myocardial infarction, nonfatal stroke, hospitalization for unstable angina, arterial revascularization, or confirmed death from cardiovascular cause. The decision to include stroke in the primary endpoint of the trial recognizes the sharing of risk factors and risk reduction with statin treatment by stroke and coronary events. For these same reasons, we used the primary endpoint as the target for associations of the CAD GRS even though it is not as well matched to the endpoint of CAD that was used in the discovery GWAS¹⁵ as is the CHD endpoint in the WHS/WGHS (above).

Genotyping

Genotype data within the WGHS has been described previously.¹³ Briefly, genome-wide genetic information was collected using the Illumina HumanHap300 duo “+” array. The subset for the current analysis includes 23,294 women of verified European ancestry. Genotype data within the approximately 71% of JUPITER participants who provided consent for genetic analysis have been described previously.¹⁸ The subset for the current analysis includes 8,749 men and women of verified European ancestry, roughly equally allocated to statin treatment or placebo. Briefly, genome-wide genetic information was collected using the Illumina Omni 1M array. Imputation of genotypes was performed with Minimac using the University of Michigan server and the Haplotype Reference Consortium reference panel for the WGHS or the 1000 Genomes (phase 3, v5) cosmopolitan reference panel for JUPITER.¹⁹ Table S1 shows demographic information for the WGHS and JUPITER samples.

Pleitoropy-Based Decomposition of a Disease GRS

SNP Selection and Coding

Genome-wide significant SNPs for T2D and CAD have been identified in recent GWASs.¹⁵^,²⁰ We considered only the independent, primary index SNPs from these studies as published, excluding SNPs arising from secondary conditional analysis. For T2D and CAD separately, the method focuses on the matrix of effects (i.e., beta coefficients) of these primary SNPs on disease-related phenotypes (see Tables S2, S3, S4, and S5 for these phenotypes, sources of the GWAS summary statistics, and a list of included primary index SNPs). The signs of the effects were determined by enforcing coding consistent with the risk allele for the related disease, except for strand-ambiguous SNPs that could not be resolved unambiguously by allele frequency and were excluded. Because the GWAS for T2D was conducted using the 1000 Genomes reference panel but GWAS summary statistics for many of the T2D-related phenotypes were available only for SNPs in the HapMap reference panel, only some of the T2D index SNPs could be included in the analysis. After additional exclusions such as coding ambiguity above and other features of the method described below, 81 of the 243 T2D index SNPs (Table S3) and 136 of the 148 CAD index SNPs (Table S5) were available for further analysis. For some applications, SNP effects for binary phenotypes from logistic regression were transformed to the liability scale as described²¹ assuming prevalence (K) of 10% for T2D and 5% for CAD.

Scaling, Imputation, and Zeroing of Beta Coefficients

For each phenotype j, beta coefficients were scaled by dividing by the phenotype’s standard deviation, i.e., the square root of the variance estimated by slope from regression of the standard errors (SE) of the phenotype’s betas across SNPs on the reciprocal of the GWAS sample size $N_{j}$ times the SNP variance $2 m a f_{i} (1 - m a f_{i})$ for SNP i with minor allele frequency $m a f_{i}$ :²²

s . e ._{i, j}^{2} \sim v a r (p h e n o_{j}) * \frac{1}{N_{j} 2 m a f_{i} (1 - m a f_{i})},

(Equation 1)

where $m a f_{i}$ was taken from GWAS summary statistics if available and from the 1000 Genomes (phase 1, v5) European subsample if not (Tables S2 and S4). Due to the nearly complete independence of the genome-wide significant index SNPs (i.e., no linkage disequilibrium [LD]), the small number of beta coefficients that were not available for some phenotypes (<10% at most) could be imputed by multiple regression using beta coefficients of other phenotypes as independent variables, although we note that more sophisticated approaches would be needed for imputation in settings with appreciable LD, e.g., Rüeger et al.²³ For secondary analysis using SNP genotypes that were standardized, each beta coefficient was further scaled by multiplying by the SNP standard deviation, i.e., $2 m a f_{i} (1 - m a f_{i})$ . A significance threshold for the beta coefficients was an adjustable parameter of the method. For each p value threshold, p_t, beta coefficients with p_beta > p_t were set to 0. However, the wide range of total number samples included in GWAS summary statistics for various phenotypes will affect power, resulting in fewer true association signals from smaller samples at each p value threshold. As an alternative to p value thresholding, we also considered significance thresholds based on the Benjamin-Hochberg false discovery rate (FDR).²⁴ The FDR will also tend to be more stringent for GWASs conducted with smaller sample sizes where power for true but weak associations will be limiting. Subsequent to thresholding, phenotypes were retained if there was at least one beta coefficient across SNPs that was not zero. Similarly, SNPs were retained if there was at least one non-zero beta coefficient across phenotypes. Since some values of p_t (or a FDR threshold) may result in assignment of beta coefficients to 0 for all phenotypes (and SNPs), this parameter influences the number of SNPs (and phenotypes) included in each decomposition.

GRS Decomposition

A GRS, $S_{disease}$ , for N_SNP disease SNPs across N_ind individuals may be computed by multiplication of the 1 × N_SNP matrix $B_{disease}$ of SNP beta coefficients for a reference trait (e.g., a disease) and the N_SNP × N_ind individuals matrix $G$ of additive encoded and centered genotypes (i.e., mean = 0):

S_{disease} = B_{disease} G .

(Equation 2)

Consider now the matrix, $B_{phenos}$ , of N_pheno phenotypes × N_SNP SNPs scaled beta coefficients related to a particular disease or reference trait (previous section). The matrix $S_{phenos}$ of GRSs for phenotypes in rows and individuals in columns may be caculated as:

S_{phenos} = B_{phenos} G .

(Equation 3)

The standard SVD of $B_{phenos}$ induces two matrices U and V, such that

B_{phenos} = U Σ V^{T}

(Equation 4)

where the left singular vectors constitute the columns of $U$ , an N_pheno × N_pheno matrix representing orthonormal linear combination of the phenotypes, while the right singular vectors constitute the columns of $V$ , an N_SNP × N_SNP matrix representing orthonormal linear combinations of the SNPs. $Σ$ is an N_pheno × N_SNP matrix with the non-zero singular values only for the p = min(N_SNP, N_pheno) diagonal elements in descending order. Usually, but not necessarily, there will be more SNPs (N_SNP) than phenotypes (N_pheno). The signs of the singular vectors are arbitrary. If the mean value of the columns of $B_{phenos}$ , corresponding to the scaled and oriented beta coefficients for each SNP across phenotypes, is zero or approximately zero (Results), then the columns of $V$ will represent (approximately) eigenvectors (“eigenSNPs” after Alter et al.²⁵) of the covariance matrix of the SNP beta coefficients times the number of SNPs, and the corresponding the squared elements of $Σ$ will represent (approximately) the eigenvalues, since

N_{S N P s} c o v (B_{phenos}) \approx B_{phenos}^{T} B_{phenos} = V Σ^{T} U^{T} U Σ V^{T} = V Σ^{2} V^{T} .

(Equation 5)

Similarly, the rows of $U$ will represent (approximately) the corresponding eigenvectors of the matrix of covariances between genetic effects on the phenotypes (“eigenphenos”). Substituting Equation 4 into Equation 3 gives $S_{phenos} = U Σ V^{T} G$ . However, we are interested in the GRSs derived from linear combinations of SNPs that represent the variance components of $B_{phenos}$ , i.e., the columns of $V$ . In particular, only the top r components (i.e., columns) with non-zero singular values explaining most, say 80%, of the variance may be relevant, corresponding to a matrix $V_{r}$ of rank r (≤ P). The r SVD-derived GRS components are therefore:

S_{SVD} = V_{r}^{T} G .

(Equation 6)

We suggest two types of sensitivity analysis for this decomposition. First, to assess the impact of the significance threshold, p_t, for setting the beta coefficients to zero on the biology represented by components of the decomposition, we identify the top phenotypes and SNPs contributing at least 80% of the loadings, i.e., norm of top loadings ≥ 0.8, on the “eigenphenos” and “eigenSNPs” across the various p value and FDR thresholds. These results are presented as tables representing Venn diagrams, i.e., “Venn tables,” that show phenotypes and SNPs shared for different thresholding strategies of the beta coefficents. Second, to assess the contribution of the correlation structure among the beta coefficients and validate the rule-of-thumb threshold of 80% cumulative variance for determining r, we implement a chi-square-like statistic and permutation procedures for evaluating its significance. Consider the distribution of eigenvalues for the SVD of matrices $B_{p h e n o s}^{p e r m}$ derived by random resampling without replacement of the scaled and zeroed beta coefficients $B_{phenos}$ across SNPs separately for each phenotype over N iterations (here N = 100,000). Once these null ordered eigenvalue distributions are scaled to their standard deviation, the nested cumulative sums of their squares will have the form of neseted sums of correlated approximate chi-square statistics. Similarly, the eigenvalues of the unpermuted betas $Σ$ from SVD of $B_{phenos}$ (Equation 4) can be scaled by the standard deviations of the empirical null distributions from the permutations and then squared and incorporated into cumulative sums. The significance of the cumulative sums is estimated by comparison to the empirical null distribution from the permultations.

Since the overall signs of the columns of $V$ are arbitrary, they were chosen such that the estimated association effects of the genetic risk scores $S_{SVD}$ were positive, i.e., odds ratio > 1. Regressing the coefficients for disease liability on these r eigenSNPs yields the residuals, $B_{resid}$ , such that

B_{disease} = \sum_{i = 1}^{r} γ_{i} V_{r_{i}}^{T} + B_{resid}

(Equation 7)

with $γ_{i}$ corresponding to the regression coefficient for the i^th row of $V_{r}^{T}$ , i.e., $V_{r_{i}}^{T},$ resulting in the decomposition of $S_{disease}$ as:

S_{disease} = B_{disease} G = (\sum_{i = 1}^{r} γ_{i} V_{i}^{T} + B_{resid}) G = \sum_{i = 1}^{r} γ_{i} V_{r_{i}}^{T} G + B_{resid} G .

(Equation 8)

SNPs from the disease or reference trait that do not strongly contribute to the reference-related phenotypes will be highly loaded in $B_{resid}$ and therefore will be respresented in the corresponding GRS_resid = $B_{resid} G$ .

Additional Scaling of Beta Coefficients

The GRSs in the rows of $S_{SVD}$ and $S_{resid}$ are not fully independent due to differences in the SNP standard deviations that depend on $m a f_{i}$ . If the beta coefficients constituting $B_{phenos}$ are scaled by SNP standard deviations, the component GRSs resulting from the SVD are now independent. However, the scaled beta coefficients represent the square root of the proportion of phenotypic variance explained by the SNPs for the various phenotypes, and therefore the interpretation of the GRS associations with disease will be different from the interpretation of the unscaled beta coefficients. See Appendix A for details.

Other Statistical Analysis

The association of GRSs with outcomes with suitable covariates (including 10 principle components of sub-European population structure as computed by EIGENSTRAT²⁶) and their interactions was performed using Cox proportional hazards regression. GRSs and covariates were standardized for all analyses. All programing was performed using R.

Ethics Statement

All analysis in the WGHS and JUPITER cohorts was approved by the local institutional review board (IRB) of the Brigham and Women’s Hospital and Partners Healthcare.

Results

Overview

An individual’s genetic risk score (GRS) for a reference phenotype, e.g., a disease, is constructed by summing the inherited genetic contribution at relevant genetic variants, e.g., SNPs, each weighted by its effect taken from an independent source, e.g., the beta coefficient from GWAS. The overall goal of the current method is to decompose such a GRS into a series of separate, largely independent GRSs each reflecting a different biological pathway (Subjects and Methods). The approach assumes that the relevant pathways are implicit in pleiotropic relationships among SNP association effects, i.e., beta coefficients from GWAS, across phenotypes related to the reference phenotype. The signs of the beta coefficients are determined by requiring that allele coding is both (1) the same across the reference and all of the related phenotypes and (2) always implies the same direction of effect for the reference, e.g., disease risk. The matrix of the beta coefficients is transformed through the SVD into orthonormal components that are interpreted as contributions of the genetic variants to separate pathways. Due to the selection of the SNPs from a reference phenotype GWAS and the enforcement of consistent allele coding, such a definition of pathways is implicitly conditional on genetic effects on the reference phenotype. The contributions of the SNPs to each of the components provide weights for constructing component, i.e., pathway-specific, GRSs Here, we apply the pleiotropy-based SVD decomposition method to GRSs derived from recent GWASs for T2D and CAD.¹⁵^,²⁰

Decomposition of the GRS for T2D

We collected association beta coefficients for the primary T2D index SNPs, i.e., excluding secondary conditional signals, from a recent GWAS²⁰ as well as 54 additional T2D-related phenotypes for which there are publicly available GWAS summary statistics (Table S2, S3, S4, and S5). After suitable encoding and scaling (Subjects and Methods), beta coefficients were set to zero when their corresponding p values were greater than a threshold, p_t, which was set to p_t = 0.05 but also explored at p_t = 0.005 and 0.5. We also explored FDR thresholds of 0.05, 0.1, and 0.75 that yielded similar numbers of non-zero beta coefficients as the p value thresholds (Subjects and Methods; Tables S3 and S6). Although the reference T2D GWAS was conducted with the 1000 Genomes reference, many of the available GWAS summary statistics for the T2D related phenotypes included only the subset of 1000 Genomes genetic variation that is represented in the HapMap reference. Thus, for threshold p_t = 0.05, 81 primary index SNPs with unambiguous coding were available across all GWAS summary statistics and had non-zero beta coefficients for at least one of the T2D-related phenotypes (Tables S3 and S6). The mean of scaled and zeroed beta coefficients for each SNP across phenotypes was approximately 0 (Figure S1). Twelve components explained 80% of the variance of the SVD decomposition, i.e., of the sum of squared eigenvalues (Figure 1). The estimated number of components based on the significance of the cumulative chi-square-like statistic was also 12, and similar across the other p value and FDR thresholds for setting beta coefficients to zero (Subjects and Methods, Table S7).

SVD of Phenotypic Effects (i.e., Beta Coefficients) for T2D Index SNPs from GWAS²⁰ across T2D-Related Phenotypes

$U_{r}$ matrix of component weights for phenotypes is shown at left, $V_{r}^{T}$ matrix of component weights for SNPs is shown on top. The squared weights of each component are show in the inset in the center. Circle area is proportional to SVD coefficients. Red and blue circles, respectively, indicate positive and negative coeffficients. SNPs are indicated using the nearest or likely gene assignment (see Table S3). In the inset, “cum.” refers to the the cumulative proportion of the total sum of the squared eigenvalues.

The loadings on the components the decomposition are informative about the biology they represent (Subjects and Methods and Venn table Table S8). In SVD of the beta coefficients thresholded by p_t = 0.05 (threshold 2 in Table S8A1), the first component emphasized direct measures of beta cell function (CIR, fasting proinsulin, disposition index, incremental insulin at 30 min, insulin at 30 min), the second component emphasized lipid metabolism (HDL-C, TG, adiponectin), and the third component emphasized metabolic measures related to insulin resistance (fatty acids, 2 h glucose, fasting insulin, HOMA-IR). Among the loci, again with p_t = 0.05 (threshold 2 in Table S8B1), the first component was dominated by SNPs rs7903146 (TCF7L2), rs7756992 (CDKAL1), and rs3802177 (SLC30A8); the second component was dominated by rs10096633 (LPL) and rs1260326 (GCKR); and the third component was dominated by rs1800961 (HNF4A) and rs1260326 (GCKR) the latter with an opposite sign compared to its contribution in the second component (Figure 1). In sensitivity analyses using different values of p_t or the FDR method for setting non-significant betas to zero, the loadings were somewhat stable, such that CIR remained a top contributors to the first component, as did HDL-C and triglycerides to the second component, and HDL-C and fatty acids to the third component (category 1, 2, and 3 in Table S8A). Top SNP contributors across thresholds were the SNPs at CDKAL1, TCF7L2 to the first component, the SNPs at GCKR and LPL to the second component, and the SNP at GCKR to the third component (category 1, 2, and 3 in Table S8B). Some consistency also persisted for higher-order components in the sensitivity analysis, e.g., phenotypes and loci that appear for all significance thresholds as well as both the p value and FDR methods for significance thresholding.

Using linear regression, we quantified the extent to which the components approximate the beta coefficients for the index T2D SNPs (Subjects and Methods). For p_t = 0.05, the first component explained 36.5% on the per allele scale or 45.3% on the standardized genotype scale of the variance in the index T2D betas, while the top 10 components explained 54% and 58% of the variance, respectively (Figure S2). T2D beta coefficients for some loci, e.g., RAI1, CTTNBP2, ZBED3, were very closely approximated by regression on the components (i.e., residual < 2.5% of absolute beta coefficient value), while others were not, e.g., ATP1B2, 140% (Figures 2 and S2 for other values of p_t). The decomposition using beta coefficients corresponding to standardized genotypes had similar properties (Figures S3–S6).

Beta Coefficients for T2D Index SNPs and Their Residuals after Regression against the Top 12 12 SVD Components of Effects on T2D-Related Phenotypes

Associations and Interactions of the Decomposed T2D GRS

In the WGHS, the GRS constructed from the beta coefficients for T2D was strongly associated with incident T2D (T2D GRS, HR, HR [95%CI] = 1.33 [1.27–1.38], p = 1.6 × 10⁻³⁸), while the association of a GRS calculated from the residuals weights after accounting for the components was much less strong (T2D resids GRS, 1.08 (1.03–1.13), 4.9 × 10⁻⁴) (Table 1). As expected, the component GRSs were minimally correlated (median absolute = 0.052, range absolute < 0.001–0.235), and essentially uncorrelated for decomposition based on standardized genotypes (Table S9, Equations A2 and A3). For the decomposition using p_t = 0.05, the effect of the component 1 GRS representing the insulin secretion measures (1.22 (1.17–1.27), 4.1 × 10⁻²⁰) was only somewhat less that the effect of the overall T2D GRS, followed by smaller significant effects for components 4, 3, 2, 7 GRSs, ordered by signficiance (Table 1). For decompositions using the other values of p_t, the component 1 GRS also had the largest effect on T2D, although the subsequent ordering of components with lesser effects varied across the thresholds (Table S10).

Table 1.

Association of a T2D GRS and Its Components with Incident T2D in the WGHS and Interaction with Triglycerides, BMI, and Physical Activity

		BMI^a						Triglycerides^a
GRS	Overall	Low	High	p int.	Low	High	p int.	Low	High	p int.
T2D.grs	1.33 (1.27–1.38), 1.6e−38	1.46 (1.31–1.62), 3e−12	1.29 (1.23–1.35), 1.4e−26	0.038	1.36 (1.29–1.44), 2.1e−27	1.28 (1.19–1.37), 9.3e−13	0.1500	1.46 (1.34–1.59), 1.8e−17	1.27 (1.21–1.33), 1.8e−21	0.007
T2D.resids.grs	1.08 (1.03–1.13), 4.9e−04	1.06 (0.96–1.18), 0.25	1.08 (1.03–1.13), 0.0013	0.757	1.05 (0.99–1.11), 0.083	1.12 (1.05–1.20), 0.0011	0.1600	1.14 (1.05–1.25), 0.0024	1.07 (1.02–1.13), 0.0067	0.251
comp.1	1.22 (1.17–1.27), 4.1e−20	1.41 (1.27–1.56), 1.3e−10	1.20 (1.15–1.26), 1.8e−14	0.009	0.77 (0.73–0.81), 2.3e−20	1.16 (1.08–1.24), 1.8e−05	0.0500	1.29 (1.18–1.41), 8.1e−09	1.18 (1.12–1.24), 3.2e−11	0.099
comp.2	1.07 (1.02–1.11), 0.0029	1.18 (1.06–1.32), 0.0021	1.06 (1.02–1.12), 0.01	0.068	1.08 (1.02–1.14), 0.0075	1.06 (0.99–1.14), 0.073	0.7900	0.99 (0.91–1.08), 0.77	1.02 (0.97–1.07), 0.40	0.458
comp.3	1.07 (1.03–1.12), 0.0016	0.92 (0.82–1.02), 0.10	1.06 (1.01–1.11), 0.024	0.019	1.06 (1.01–1.13), 0.029	1.09 (1.02–1.17), 0.017	0.7100	1.16 (1.06–1.26), 0.0012	1.07 (1.01–1.12), 0.011	0.114
comp.4	1.13 (1.08–1.18), 3.2e−08	1.09 (0.98–1.21), 0.12	1.14 (1.08–1.19), 1.4e−07	0.489	1.15 (1.09–1.21), 1.1e−06	1.09 (1.02–1.17), 0.012	0.2500	1.17 (1.07–1.27), 5.2e−04	1.13 (1.08–1.19), 5.7e−07	0.548
comp.5	1.01 (0.97–1.06), 0.56	1.04 (0.93–1.15), 0.52	0.98 (0.93–1.03), 0.36	0.293	1.05 (0.99–1.11), 0.098	0.96 (0.90–1.03), 0.28	0.0640	1.02 (0.94–1.12), 0.58	1.00 (0.95–1.05), 0.95	0.613
comp.6	1.05 (1.01–1.10), 0.024	1.06 (0.95–1.17), 0.33	1.04 (0.99–1.09), 0.084	0.788	1.05 (1.00–1.11), 0.069	1.06 (0.99–1.13), 0.11	0.8900	1.06 (0.97–1.16), 0.20	1.04 (0.99–1.09), 0.12	0.797
comp.7	1.06 (1.02–1.11), 0.0044	0.99 (0.89–1.10), 0.86	1.04 (0.99–1.09), 0.091	0.500	1.07 (1.01–1.13), 0.018	1.07 (1.00–1.15), 0.038	0.8300	1.03 (0.94–1.12), 0.57	1.06 (1.01–1.11), 0.023	0.567
comp.8	1.02 (0.97–1.06), 0.46	1.06 (0.95–1.17), 0.32	1.01 (0.96–1.06), 0.78	0.418	1.02 (0.96–1.08), 0.49	1.02 (0.95–1.09), 0.56	0.9800	1.02 (0.94–1.11), 0.64	1.01 (0.96–1.06), 0.84	0.813
comp.9	1.04 (0.99–1.08), 0.093	0.98 (0.88–1.09), 0.73	1.03 (0.99–1.08), 0.18	0.373	1.02 (0.97–1.08), 0.42	1.05 (0.98–1.13), 0.13	0.5000	1.02 (0.94–1.12), 0.58	1.04 (0.99–1.09), 0.16	0.831
comp.10	1.01 (0.96–1.05), 0.74	1.05 (0.94–1.17), 0.37	1.00 (0.96–1.05), 0.88	0.446	1.07 (1.01–1.13), 0.02	0.92 (0.86–0.98), 0.013	0.0009	0.97 (0.89–1.06), 0.51	1.02 (0.97–1.07), 0.43	0.261
comp.11	1.07 (1.02–1.12), 0.0022	1.13 (1.02–1.26), 0.02	1.06 (1.01–1.11), 0.018	0.207	1.12 (1.06–1.19), 4e−05	0.98 (0.92–1.05), 0.57	0.0029	1.13 (1.04–1.23), 0.0059	1.06 (1.01–1.11), 0.023	0.203
comp.12	1.02 (0.98–1.07), 0.29	1.01 (0.91–1.13), 0.82	1.02 (0.98–1.07), 0.33	0.859	1.01 (0.96–1.07), 0.64	1.04 (0.97–1.11), 0.28	0.5600	1.00 (0.92–1.09), 0.96	1.01 (0.97–1.07), 0.58	0.811

Open in a new tab

Hazard ratios, 95% confidence intervals, and p values from separate Cox-proportional hazards models for each GRS adjusted for baseline age and 10 principal components of sub-European population structure. “Overall” refers to the entire population. High and low refer to analysis stratified by each of the interaction variables. P int. refers to the significance of the interaction between each specified component and the interaction variable in analysis on the entire population. Components derived from SVD of phenotype beta coefficients using p_t = 0.05 (Subjects and Methods).

BMI and triglyceride strata were defined by median values.

Physical activity strata defined by inactive (<3.75 mets) v. active (≥3.75 mets).

The components revealed biological functions that may be particularly relevant to interactions between genetics and three clinical variables that are very strongly associated with incident T2D in the WGHS: BMI, triglyceride levels, and physical activity (Table S11). There were significant interactions of the overall T2D GRS with BMI (p_int = 0.038) and triglycerides (p_int = 0.007) but not with physical activity (p_int = 0.15) (Table 1). For p_t = 0.05, the interaction with BMI was recapitulated by components 1 (p_int = 0.009) and 3 (p_int = 0.019) GRSs, emphasizing glucose-, insulin-, and HDL-C-related traits. No significant interaction was observed in the GRS from the residuals. GRS components 1 and 3 had more significant interactions with BMI than the overall GRS with BMI: higher compared to lower BMI attenuated the absolute effect of the overall T2D GRS (1.29 [1.23–1.35] v. 1.46 [1.31–1.62]) and the component 1 GRS (1.20 [1.15–1.26] v. 1.41 [1.27–1.56]), but there were opposite effects for the component 3 GRS (1.06 [1.01–1.11] v. 0.92 [0.82–1.02]). When the SVD decomposition was applied to beta coefficients zeroed by an FDR threshold of 0.1 resulting in a comparable number non-zero beta coefficients as a p value threshold (p_t) of 0.05 (Tables S6 and S12), the overall effects of components 1 and 3 were similar while the interaction term for component 1 was an order of magnitutde more significant but now attenuated for component 3, p_int = 0.00075 and 0.13, respectively.

For physical activity, there was no significant interaction with the overall T2D GRS but the component 10 GRS (which emphasizes glucose-related and fatty acid measures) interacted strongly (p_int = 0.001), showing opposite effects in the two strata (0.92 [0.86–0.98] v. 1.07 [1.01–1.13]). The component 11 GRS (which shares phenotypes with component 10 but often differs in sign) also showed a significant interaction due to exclusive effects in the low physical activity stratum (1.12 [1.06–1.19] v. 0.98 [0.92–1.05], p_int = 0.003), while the differences between strata for the component GRS 1 were just significant (1.16 [1.08–1.24] v. 0.77 [0.73–0.81], p_int = 0.05). However, a significant interaction between triglyceride levels and the overall T2D GRS was not recapitulated by the GRS components. There were both significant overall effects and significant interactions using component GRSs arising from beta coefficients derived from the standardized genotypes (Table S13). The component 1 GRS showed the strongest overall effect (1.22 [1.17–1.28], p = 1.8 × 10⁻²⁰), as well as quantitative interactions with BMI (p_int = 0.008) and physical activity (p_int = 0.05). With the standardized genotypes, the component 2 GRS also interacted with BMI (p_int = 0.03), while component GRSs 8 and 12 interacted with physical activity, but there were no interactions with triglycerides (Table S13).

Decomposition of the GRS for CAD

We also decomposed a GRS including 130 qualifying SNPs from the recent CAD GWAS¹⁵ through pleiotropy among LDL-C, HDL-C, triglycerides, C-reactive protein (CRP), systolic and diastolic blood pressure, eGFRcrea, and type 2 diabetes with and without adjustment for BMI. Most (i.e., >80%) of the variance among the beta coefficients for the CAD-related phenotypes was explained by the first two or three components, although all components are depicted in the figure (Figures 3 and S7, Tables S14 and S15). The components explain less of the variance in the beta coefficients for CAD than for the decomposition of T2D but still reached 34% with all 9 (Figure S8). For p_t = 0.05, the top component emphasized LDL-C, the second component emphasized reciprocal effects on triglycerides and HDL-C, and the third component emphasized blood pressure (Figure 3). These loadings were stable across significance threholds (of p_t or the FDR) for setting beta coefficients to zero. Additional components included effects also for CRP and T2D, as well as contrasting signs for some phenotypes, e.g., component 3 v. 8, component 6 v. 9. The top SNP contributions were rs7412 (APOE) and rs11591147 (PCSK9) for component 1, rs116843064 (ANGPTL4) and rs17091891 (LPL) for component 2, and eight SNPs for component 3 (Table S16). The extent to which individual SNP beta coefficients for CAD are reconstituded from the components is less than for T2D, but is nevertheless substantial for some (Figure 4). When based on beta coefficients scaled for standardized genotypes, proportion of the CAD beta coefficients explained by the components was only 12% but other aspects were similar to the decomposition based on per allele betas (Figures S9–S12).

SVD of Phenotypic Effects (i.e., Beta Coefficients) for CAD Index SNPs from GWAS¹⁵ across CAD-Related Phenotypes

Presentation as in Figure 1. SNPs are indicated using the nearest or likely gene assignment (see Table S5).

Beta Coefficients for CAD Index SNPs and Their Residuals after rRegression against All (i.e., 9) SVD Components of Effects on CAD-Related Phenotypes

Associations and Interactions of the Decomposed CAD GRS

In the WGHS, GRSs derived from the pleiotropy-based SVD components for the CAD SNPs were minimally correlated (median absolute = 0.042, range absolute = 0.002–0.282) and not correlated for beta coefficients derived from beta coefficients scaled for standardized genotypes (Table S17). While the GRS for CAD overall was strongly associated with incident CHD (1.40 [1.32–1.49], p = 3.3 × 10⁻³⁰), the decomposition for p_t = 0.05 indicated that much of the effect was due to effects captured by the component 2 GRS (triglycerides and HDL-C, 1.16 [1.10–1.23], p = 7.5 × 10⁻⁷) and the component 3 GRS (blood pressure, trigycerides, 1.15 [1.09–1.22], p = 2.7 × 10⁻⁶) rather than the component 1 GRS (1.05 [0.99–1.11], p = 0.12), which is largely related to LDL-C (Table 2). The GRS constructed from weights in the residuals also remained associated with incident CHD (1.12 [1.05–1.18], p = 1.9 × 10⁻⁴), suggesting that the CAD-related phenotypes do not fully explain the effects of the CAD GWAS SNPs, consistent with the moderate variance in the CAD betas explained by the SVD components (Figure S8). Results were similar for other values of p_t (Table S18) and for analysis corresponding to standardized genotypes (Table S19).

Table 2.

Association of a CAD GRS and Its Components with the CHD in the WGHS

GRS	Overall
CAD.grs	1.40 (1.32–1.49), 3.3e−30
CAD.resids.grs	1.12 (1.05–1.18), 1.9e−04
comp.1	1.05 (0.99–1.11), 0.12
comp.2	1.16 (1.10–1.23), 7.5e−07
comp.3	1.15 (1.09–1.22), 2.7e−06
comp.4	1.04 (0.99–1.11), 0.14
comp.5	1.02 (0.96–1.08), 0.52
comp.6	1.08 (1.02–1.14), 0.012
comp.7	1.04 (0.98–1.10), 0.23
comp.8	1.05 (0.99–1.11), 0.12
comp.9	1.00 (0.95–1.07), 0.87

Open in a new tab

Hazard ratios, 95% confidence intervals, and p values from Cox-proportional hazards models adjusted for baselne age and 10 principal components of sub-European population structure. Components derived from SVD of phenotype beta coefficients using p_t = 0.05 (Subjects and Methods).

In JUPITER, a trial of rosuvastatin v. placebo for primary prevention, the correlations of the CAD-related GRSs were minimal and similar to what was observed in the WGHS (Table S20). There was no significant association of the CAD GRS on the incident primary event, a composite of MI, stroke, angina, revascularization, and cardiovascular death (1.12 [0.98–1.27], p = 0.096) in models adjusted for drug allocation (Table 3). However, as in the WGHS, a genetic effect was concentrated in the component 2 GRS (1.19 [1.04–1.36], p = 0.01) and to a lesser extent with marginal significance in the component 3 GRS (1.12 [0.99–1.28], p = 0.08). Overall effects were similarly concentrated in the component 2 GRS for other values of p_t (Table S21) and in analysis corresponding to standardized genotypes (Table S22). While there is an overwhelming effect of rosuvastatin in this population in reducing the incidence of the primary event (Table S11 and Ridker et al.¹⁷), the overall CAD GRS only showed a marginal interaction rosuvastatin v. placebo allocation (p_int = 0.06) (Table 3). However, there were significant interactions for the component 1 (p = 0.018, LDL-C) and component 2 (p = 0.042, triglycerides and HDL-C) GRSs but not the CAD residual GRS. The interactions were due to strongly significant associations in the placebo arm, which was also observed for the overall CAD GRS, and essentially null associations in the rosuvastatin arm. The interactions persisted in analysis adjusted for LDL-C, HDL-C, and triglyceride levels at baseline (Table 3, see also Table S22).

Table 3.

Association of a CAD GRS and Its Components with the Primary Endpoint of JUPITER and Interactions with Randomized Statin Allocation

		Stratified by Statin Allocation			Stratified by Statin Allocation, Baseline Lipid^aAdjusted
GRS	Overall^b	Placebo	Statin	p int.	Placebo	Statin	p int.
CAD.grs	1.12 (0.98–1.27), 0.096	1.22 (1.03–1.44), 0.018	0.95 (0.77–1.18), 0.65	0.064	1.21 (1.02–1.42), 0.026	0.96 (0.78–1.18), 0.69	0.079
CAD.resids.grs	0.92 (0.80–1.05), 0.20	1.19 (1.01–1.40), 0.04	0.98 (0.80–1.21), 0.89	0.149	1.18 (1.00–1.38), 0.054	0.99 (0.80–1.22), 0.90	0.172
comp.1	1.07 (0.94–1.22), 0.30	1.21 (1.02–1.44), 0.028	0.87 (0.72–1.07), 0.19	0.018	1.20 (1.00–1.43), 0.047	0.87 (0.71–1.07), 0.20	0.019
comp.2	1.19 (1.04–1.36), 0.01	1.32 (1.12–1.57), 0.0013	1.00 (0.81–1.23), 0.98	0.042	1.29 (1.09–1.54), 0.0038	0.95 (0.77–1.18), 0.66	0.042
comp.3	1.12 (0.99–1.28), 0.075	1.21 (1.02–1.42), 0.026	1.00 (0.81–1.23), 0.99	0.166	1.21 (1.03–1.43), 0.022	1.00 (0.81–1.24), 0.99	0.167
comp.4	1.03 (0.91–1.18), 0.60	1.01 (0.86–1.18), 0.94	1.09 (0.88–1.34), 0.45	0.631	1.00 (0.84–1.17), 0.96	1.11 (0.90–1.38), 0.33	0.571
comp.5	1.04 (0.91–1.19), 0.54	1.06 (0.90–1.25), 0.51	1.01 (0.82–1.24), 0.95	0.666	1.06 (0.90–1.25), 0.50	0.99 (0.80–1.23), 0.96	0.614
comp.6	1.04 (0.91–1.18), 0.57	1.00 (0.85–1.18), 0.99	1.10 (0.90–1.36), 0.35	0.456	1.00 (0.85–1.18), 0.98	1.12 (0.90–1.38), 0.31	0.427
comp.7	1.03 (0.90–1.17), 0.71	0.99 (0.84–1.17), 0.91	1.10 (0.89–1.36), 0.38	0.48	0.99 (0.84–1.17), 0.92	1.10 (0.89–1.36), 0.40	0.451
comp.8	1.10 (0.97–1.25), 0.15	1.11 (0.94–1.31), 0.22	1.10 (0.89–1.37), 0.38	0.831	1.12 (0.95–1.31), 0.19	1.11 (0.90–1.38), 0.34	0.824
comp.9	1.14 (1.00–1.30), 0.047	1.07 (0.91–1.26), 0.41	1.26 (1.02–1.57), 0.033	0.202	1.07 (0.91–1.26), 0.40	1.27 (1.02–1.57), 0.031	0.195

Open in a new tab

Hazard ratios, 95% confidence intervals, and p values from Cox-proportional hazards models minimally adjusted for age, sex, and 10 principal components of sub-European population structure or additional adjustment for baseline lipids as indicated. P int. refers to the p value for the interaction between each specified component and randomized rosuvastatin allocation. Components derived from SVD of phenotype beta coefficients using p_t = 0.05 (Subjects and Methods).

Baseline LDL-C, HDL-C, triglycerides.

Additional adjustment for randomized statin allocation.

Discussion

We describe a decomposition of SNP weights used to construct a GRS for a reference trait, e.g., a disease, through SVD applied to the scaled effects of SNPs for related phenotypes, i.e., betas, combined with residual effects. A critical element of the procedure is enforcement of a consistent orientation of SNP coding across the reference trait and the related phenotypes such that the components and therefore the pathways they may represent are defined with respect to biological process resulting in the reference trait, e.g., disease risk. Thus, the choice of SNPs and their pleiotropic relationships can be viewed as conditional on the effects on the disease outcome. By application to T2D and CAD, we demonstrate the potential of the components to reveal biological pathways through the relative weights and signs of both SNPs and phenotypes for each component; and we provide semi-formal approaches for assessing the sensitivity of the decomposition to significance thresholds for the beta coefficients. Applying the component weights in constructing GRSs, we assess pathway-specific contributions to an individual’s genetic risk of T2D or CAD. The decomposition may be particularly useful for revealing gene-by-covariate interactions that are mediated primarily by specific pathways and therefore not ideally manifest when tested using a GRS that aggregates all genetic etiologies. While essentially all of the overall effects in the WGHS would have survived significance standards consistent with multiple testing, this was not uniformly true for the associations in JUPITER or for the interactions. We note, however, there are strong prior hypotheses for the interactions, which may mitigate the need to account for multiple testing. Replication in future studies would strengthen the findings.

This work complements an application of non-negative matrix factorization of T2D-related genetic associations.¹² Compared to that work, we updated some GWASs (e.g., for T2D itself) while also excluding phenotypes for which public GWAS summary statistics did not include standard errors or sample size, both necessary for scaling. We also excluded GWASs corresponding to analysis stratified by sex to simplify the interpretation. Methodologically, the other work formulated the issue of decomposition in terms of soft clustering of SNPs in a Bayesian framework rather than the SVD and identified non-exclusive assignments of disease variants to components, which were interpreted as disease etiologies. That method also estimated an optimal number of such components, compared to the choice of “r” in the current approach via the proportion of variance explained or the significance of the cumulative chi-square-like statistic . Although the other work was primarily focused on clustering loci rather than formally decomposing a risk score, it derived GRSs through empirical thresholds for assignment of specific variants using cluster membership weights as opposed to inclusion of all SNPs through decomposition of weights here. For T2D, we do not find an exact correspondence of the components from the two approaches. However, our component 1 is similar to their “beta-cell” component, and collectively, our components 2 and 3 are reminiscent of their “lipodystrophy” and “liver/lipid” components.

In both the WGHS and JUPITER, there was risk associated with a GRS component related to LDL-C levels and a second component with reciprocal effects on triglycerides and HDL-C, consistent with previous work on lipid risk in atherosclerosis.²⁷ Perhaps surprisingly, this second component was at least as important as the first for the pharmacogenetic effect on outcomes with rosuvastatin treatment. The prominent loci for the first component, likely affecting PCSK9 (rs11591147) and APOE (rs7412), are known determinants of LDL-C response to rosuvastatin.¹⁸ The results here suggest that genetic effects at LPL (rs17091891) and ANGPTL4 (rs116843064) may be equally important for rosuvastatin-based risk reduction. Effects of both components were minimally attenuated by adjustment for baseline LDL-C, HDL-C, and triglycerides.

Similarly, interactions between the component GRSs for incident T2D in the WGHS may reveal biology. For BMI, the quantitative interaction effect with GRS component 1 implies that elevated BMI attenuates genetic effects on T2D related to beta-cell function, perhaps via a dysregulation phenomenon. By contrast, for component GRS component 3, the qualitative nature of the interaction effect may suggest persistent but opposite genetic effects in the two BMI strata, and focuses attention on biology shared by HOMA-IR and HOMA-B that is also related, in a reciprocal manner, to HDL-C but does not involve triglycerides. For physical activity, the strong qualitative interaction with the component 10 GRS implicates biology that is distributed over fatty acid and glucose metabolism, driven by a top genetic effect implicating PNPLA3 (rs738408), which is known for involvement with liver function.²⁸

The approach here represents an example of the general utility of matrix factorization, including SVD, for understanding complex biological processes and leveraging this understanding for clinical translation.²⁵^,²⁹^,³⁰ The methods have typically been used to resolve gene expression patterns conditional on cell/tissue type or to delineate somatic determinants of oncogenesis. Application to GRSs for complex trait genetics through pleiotropic relationships is relatively less developed formally, although it has intuitive appeal.³¹^,³²^,³³ In general, whether an orthogonal variance decomposition approach, as is done here, or a factorization method as was done previously for T2D¹² is superior will depend on exact relationships among the disease and the disease-related phenotypes. For example, clustering may be superior if the system can be modeled by just a few dominant pathways or non-linear relationships, whereas SVD may be superior when the pathways are essentially orthogonal (as for CAD) or for capturing finer structure involving minor but still biologically important pathways that do not contribute sufficiently to justify augmenting the the number of clusters in the clustering approach. Moreover, the SVD method has a more direct and principled connection to the construction of GRSs from the decomposed beta coefficient matrix, while also enforcing the property of orthogonality that aids interpretation. It is also more tractable computationally, although this aspect may have only minimal practical consequences. Regardless, any method that models genetic associations with a GRS will have the potential limitation in assuming that linear combinations of genetics without interactions are suitable representations for underlying biological functions. This assumption may be reasonable, although possibly incomplete, for the common genetic variation discovered from GWAS of complex traits. A highly polygenetic genetic model of small additive, i.e., linear, effects without interactions of common GWAS SNPs tends to explain much of complex trait heritability.³⁴

The availability of disease-related phenotypes and the precision of pleiotropic genetic effects are critical to the implementation of the current approach, and yet both remain sub-optimal in the GWAS literature. For example, by comparison to T2D, CAD SNP effects were less well explained by the disease-related phenotypes, likely in part due to a relative dearth of choices for pleiotropic phenotypes and corresponding GWASs. Nevertheless, recent progress with GWAS across all phenotypic information in the UKBiobank demonstrates the potential for genetic analysis to provide more complete accounting of pleiotropy in large cohorts with comprehensive phenotyping; and the potential of complete clinical profiles and genetics on entire populations is expected to continue to grow as genetic analysis is broadly deployed in health care systems as proposed (see Web Resources).³⁵^,³⁶^,³⁷^,³⁸^,³⁹^,⁴⁰^,⁴¹^,⁴² As GRSs become more refined through more accurate and precise estimation of SNP effects in large, multi-ethnic samples, they may begin to complement conventional clinical measures in standard clinical practice.¹^,⁴³ By isolating the biological pathways implicitly reflected in pleiotropic relationships, the method described here or elaborations thereof may contribute to how GRSs are implemented in optimizing healthcare with genomics.

Acknowledgments

The authors gratefully acknowledge Professors Julie E. Buring and Paul M Ridker for access to data in the WGHS and JUPITER. This work was funded in part by grants supporting the Women’s Health Study from the NHLBI/NIH (HL043851 and HL080467) and the NCI/NIH (CA047988 and UM1CA182913), Julie Buring and I-Min Lee, PIs, and the Gene-by-Lifestyle Consortium (HL118305 from the NHLBI/NIH, Rao PI, D.I.C. PI of subcontract). Funding for genotyping in the WGHS was provided by Amgen (Ridker, D.I.C., PIs). Funding for genotyping in JUPITER was provided by AstraZeneca (Ridker, D.I.C., PIs). O.V.D. is supported by grant 5K01HL135342 from the NHLBI/NIH. M.S.U. is supported by K23DK114551 from the NIDDK/NIH.

Published: April 16, 2020

Footnotes

Supplemental Data can be found online at https://doi.org/10.1016/j.ajhg.2020.03.011.

Appendix A: Correlation of Component GRSs and Standardization of Beta Coefficients

While the columns of $V$ are orthonormal, the component risk scores $S$ are not independent. Consider the covariance of between GRSs for components i and j over k individuals:

c o v (S_{i} S_{j}) = (\frac{1}{N}) \sum_{k} (S_{i, k} - \overline{S_{i}}) (S_{j, k} - \overline{S_{j}})

= (\frac{1}{N}) \sum_{k} S_{i, k} S_{j, k} - S_{i, k} \overline{S_{j}} - S_{j, k} \overline{S_{i}} + \overline{S_{i}} {\overline{S}}_{j}

(Equation A1)

= (\frac{1}{N}) \sum_{k} S_{i, k} S_{j, k}

Since

\sum_{k} S_{i, k} = \sum_{k} \sum_{l} V_{i, l} G_{l, k} = \sum_{l} V_{i, l} \sum_{k} G_{l, k} = 0

for centered encoded genotypes (i.e., mean = 0). Thus,

c o v (S) = (\frac{1}{N}) V G {(V G)}^{T} = (\frac{1}{N}) V G G^{T} V^{T}

= V diag (v a r (g_{1}), v a r (g_{2}), \dots, v a r (g_{N S N P s})) V^{T}

(Equation A2)

when the SNPs are independent, i.e., no LD, with $v a r (g_{1})$ equalling the variance in the encoded genotypes for SNP i. $V a r (g_{i})$ will differ for each SNP i according to minor allele frequency, giving rise to non-zero off-diagonal terms of $c o v (S)$ , i.e., non-zero covariation between component risk scores.

If the genotypes are standardized such that $v a r (g_{i}) = 1$ for all SNPs, i and the beta coefficients in matrix $B_{phenos}$ from Equation 3 in the Subjects and Methods are scaled accordingly (such that the effect per scaled allele is the same as the effects with betas and alleles without standardization), then the SVD of the standardized coefficients gives rise to new right singular vectors $V_{std}$ , and Equation A2 reduces to

= V_{std} V_{std}^{T} = I

(Equation A3)

indicating no correlation between the risk scores for beta coefficients scaled to standardized genotypes.

However, such standardized beta coefficients will no longer represent the mean change in phenotype per allele. Instead, for phenotypes with unit variance (as enforced by the method), the scaled betas will be the square root of the proportion variance phenotype explained by each (standardized) SNP, which is a population-based measure whose value will have been influenced by population dynamics, selection, etc. The SVD of these standardized beta coefficients will be non-trivially different from the SVD of the beta coefficients without standardization, since the relative scale factors, i.e., the standard deviation in SNP genotypes, varies from 0.31 to 0.71 for SNPs with MAF 0.05–0.5 (Subjects and Methods). Therefore, interpretations of the decompositions based on the different representations of the beta coefficients are different.

Web Resources

Neale lab summary statistics for associations in the UK Biobank, http://www.nealelab.is/blog/2017/7/19/rapid-gwas-of-thousands-of-phenotypes-for-337000-samples-in-the-uk-biobank
R statistical software, https://www.r-project.org/

Declaration of Interests

The authors declare no competing interests.

Supplemental Data

Document S1. Figures S1–S12 and Tables S1, S3, and S5–S22

mmc1.pdf^{(1.9MB, pdf)}

Data S1. Tables S2 and S4

mmc2.xlsx^{(62.3KB, xlsx)}

Document S2. Article plus Supplemental Information

mmc3.pdf^{(3MB, pdf)}

References

1.Chatterjee N., Shi J., García-Closas M. Developing and evaluating polygenic risk prediction models for stratified disease prevention. Nat. Rev. Genet. 2016;17:392–406. doi: 10.1038/nrg.2016.27. [DOI] [PMC free article] [PubMed] [Google Scholar]
2.Aschard H. A perspective on interaction effects in genetic association studies. Genet. Epidemiol. 2016;40:678–688. doi: 10.1002/gepi.21989. [DOI] [PMC free article] [PubMed] [Google Scholar]
3.Subramanian A., Tamayo P., Mootha V.K., Mukherjee S., Ebert B.L., Gillette M.A., Paulovich A., Pomeroy S.L., Golub T.R., Lander E.S., Mesirov J.P. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. USA. 2005;102:15545–15550. doi: 10.1073/pnas.0506580102. [DOI] [PMC free article] [PubMed] [Google Scholar]
4.Kanehisa M., Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30. doi: 10.1093/nar/28.1.27. [DOI] [PMC free article] [PubMed] [Google Scholar]
5.Pers T.H., Karjalainen J.M., Chan Y., Westra H.J., Wood A.R., Yang J., Lui J.C., Vedantam S., Gustafsson S., Esko T., Genetic Investigation of ANthropometric Traits (GIANT) Consortium Biological interpretation of genome-wide association studies using predicted gene functions. Nat. Commun. 2015;6:5890. doi: 10.1038/ncomms6890. [DOI] [PMC free article] [PubMed] [Google Scholar]
6.The Gene Ontology Consortium The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res. 2019;47(D1):D330–D338. doi: 10.1093/nar/gky1055. [DOI] [PMC free article] [PubMed] [Google Scholar]
7.Ashburner M., Ball C.A., Blake J.A., Botstein D., Butler H., Cherry J.M., Davis A.P., Dolinski K., Dwight S.S., Eppig J.T., The Gene Ontology Consortium Gene ontology: tool for the unification of biology. Nat. Genet. 2000;25:25–29. doi: 10.1038/75556. [DOI] [PMC free article] [PubMed] [Google Scholar]
8.Liu J.Z., McRae A.F., Nyholt D.R., Medland S.E., Wray N.R., Brown K.M., Hayward N.K., Montgomery G.W., Visscher P.M., Martin N.G., Macgregor S., AMFS Investigators A versatile gene-based test for genome-wide association studies. Am. J. Hum. Genet. 2010;87:139–145. doi: 10.1016/j.ajhg.2010.06.009. [DOI] [PMC free article] [PubMed] [Google Scholar]
9.de Leeuw C.A., Mooij J.M., Heskes T., Posthuma D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput. Biol. 2015;11:e1004219. doi: 10.1371/journal.pcbi.1004219. [DOI] [PMC free article] [PubMed] [Google Scholar]
10.Segrè A.V., Groop L., Mootha V.K., Daly M.J., Altshuler D., DIAGRAM Consortium. MAGIC investigators Common inherited variation in mitochondrial genes is not enriched for associations with type 2 diabetes or related glycemic traits. PLoS Genet. 2010;6:6. doi: 10.1371/journal.pgen.1001058. [DOI] [PMC free article] [PubMed] [Google Scholar]
11.Mahajan A., Wessel J., Willems S.M., Zhao W., Robertson N.R., Chu A.Y., Gan W., Kitajima H., Taliun D., Rayner N.W., ExomeBP Consortium. MAGIC Consortium. GIANT Consortium Refining the accuracy of validated target identification through coding variant fine-mapping in type 2 diabetes. Nat. Genet. 2018;50:559–571. doi: 10.1038/s41588-018-0084-1. [DOI] [PMC free article] [PubMed] [Google Scholar]
12.Udler M.S., Kim J., von Grotthuss M., Bonàs-Guarch S., Cole J.B., Chiou J., Boehnke M., Laakso M., Atzmon G., Glaser B., Christopher D. Anderson on behalf of METASTROKE and the ISGC Type 2 diabetes genetic loci informed by multi-trait associations point to disease mechanisms and subtypes: A soft clustering analysis. PLoS Med. 2018;15:e1002654. doi: 10.1371/journal.pmed.1002654. [DOI] [PMC free article] [PubMed] [Google Scholar]
13.Ridker P.M., Chasman D.I., Zee R.Y., Parker A., Rose L., Cook N.R., Buring J.E., Women’s Genome Health Study Working Group Rationale, design, and methodology of the Women’s Genome Health Study: a genome-wide association study of more than 25,000 initially healthy american women. Clin. Chem. 2008;54:249–255. doi: 10.1373/clinchem.2007.099366. [DOI] [PubMed] [Google Scholar]
14.Nikpay M., Goel A., Won H.H., Hall L.M., Willenborg C., Kanoni S., Saleheen D., Kyriakou T., Nelson C.P., Hopewell J.C. A comprehensive 1,000 Genomes-based genome-wide association meta-analysis of coronary artery disease. Nat. Genet. 2015;47:1121–1130. doi: 10.1038/ng.3396. [DOI] [PMC free article] [PubMed] [Google Scholar]
15.van der Harst P., Verweij N. Identification of 64 Novel Genetic Loci Provides an Expanded View on the Genetic Architecture of Coronary Artery Disease. Circ. Res. 2018;122:433–443. doi: 10.1161/CIRCRESAHA.117.312086. [DOI] [PMC free article] [PubMed] [Google Scholar]
16.Liu S., Lee I.M., Song Y., Van Denburgh M., Cook N.R., Manson J.E., Buring J.E. Vitamin E and risk of type 2 diabetes in the women’s health study randomized controlled trial. Diabetes. 2006;55:2856–2862. doi: 10.2337/db06-0456. [DOI] [PubMed] [Google Scholar]
17.Ridker P.M., Danielson E., Fonseca F.A., Genest J., Gotto A.M., Jr., Kastelein J.J., Koenig W., Libby P., Lorenzatti A.J., MacFadyen J.G., JUPITER Study Group Rosuvastatin to prevent vascular events in men and women with elevated C-reactive protein. N. Engl. J. Med. 2008;359:2195–2207. doi: 10.1056/NEJMoa0807646. [DOI] [PubMed] [Google Scholar]
18.Chasman D.I., Giulianini F., MacFadyen J., Barratt B.J., Nyberg F., Ridker P.M. Genetic determinants of statin-induced low-density lipoprotein cholesterol reduction: the Justification for the Use of Statins in Prevention: an Intervention Trial Evaluating Rosuvastatin (JUPITER) trial. Circ Cardiovasc Genet. 2012;5:257–264. doi: 10.1161/CIRCGENETICS.111.961144. [DOI] [PubMed] [Google Scholar]
19.Das S., Forer L., Schönherr S., Sidore C., Locke A.E., Kwong A., Vrieze S.I., Chew E.Y., Levy S., McGue M. Next-generation genotype imputation service and methods. Nat. Genet. 2016;48:1284–1287. doi: 10.1038/ng.3656. [DOI] [PMC free article] [PubMed] [Google Scholar]
20.Mahajan A., Taliun D., Thurner M., Robertson N.R., Torres J.M., Rayner N.W., Payne A.J., Steinthorsdottir V., Scott R.A., Grarup N. Fine-mapping type 2 diabetes loci to single-variant resolution using high-density imputation and islet-specific epigenome maps. Nat. Genet. 2018;50:1505–1513. doi: 10.1038/s41588-018-0241-6. [DOI] [PMC free article] [PubMed] [Google Scholar]
21.Gillett A.C., Vassos E., Lewis C.M. Transforming Summary Statistics from Logistic Regression to the Liability Scale: Application to Genetic and Environmental Risk Scores. Hum. Hered. 2018;83:210–224. doi: 10.1159/000495697. [DOI] [PubMed] [Google Scholar]
22.Winkler T.W., Day F.R., Croteau-Chonka D.C., Wood A.R., Locke A.E., Mägi R., Ferreira T., Fall T., Graff M., Justice A.E., Genetic Investigation of Anthropometric Traits (GIANT) Consortium Quality control and conduct of genome-wide association meta-analyses. Nat. Protoc. 2014;9:1192–1212. doi: 10.1038/nprot.2014.071. [DOI] [PMC free article] [PubMed] [Google Scholar]
23.Rüeger S., McDaid A., Kutalik Z. Evaluation and application of summary statistic imputation to discover new height-associated loci. PLoS Genet. 2018;14:e1007371. doi: 10.1371/journal.pgen.1007371. [DOI] [PMC free article] [PubMed] [Google Scholar]
24.Benjamini Y., Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. B. 1995;57:289–300. [Google Scholar]
25.Alter O., Brown P.O., Botstein D. Singular value decomposition for genome-wide expression data processing and modeling. Proc. Natl. Acad. Sci. USA. 2000;97:10101–10106. doi: 10.1073/pnas.97.18.10101. [DOI] [PMC free article] [PubMed] [Google Scholar]
26.Price A.L., Patterson N.J., Plenge R.M., Weinblatt M.E., Shadick N.A., Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 2006;38:904–909. doi: 10.1038/ng1847. [DOI] [PubMed] [Google Scholar]
27.Do R., Willer C.J., Schmidt E.M., Sengupta S., Gao C., Peloso G.M., Gustafsson S., Kanoni S., Ganna A., Chen J. Common variants associated with plasma triglycerides and risk for coronary artery disease. Nat. Genet. 2013;45:1345–1352. doi: 10.1038/ng.2795. [DOI] [PMC free article] [PubMed] [Google Scholar]
28.Speliotes E.K., Yerges-Armstrong L.M., Wu J., Hernaez R., Kim L.J., Palmer C.D., Gudnason V., Eiriksdottir G., Garcia M.E., Launer L.J., NASH CRN. GIANT Consortium. MAGIC Investigators. GOLD Consortium Genome-wide association analysis identifies variants associated with nonalcoholic fatty liver disease that have distinct effects on metabolic traits. PLoS Genet. 2011;7:e1001324. doi: 10.1371/journal.pgen.1001324. [DOI] [PMC free article] [PubMed] [Google Scholar]
29.Hore V., Viñuela A., Buil A., Knight J., McCarthy M.I., Small K., Marchini J. Tensor decomposition for multiple-tissue gene expression experiments. Nat. Genet. 2016;48:1094–1100. doi: 10.1038/ng.3624. [DOI] [PMC free article] [PubMed] [Google Scholar]
30.Stein-O’Brien G.L., Arora R., Culhane A.C., Favorov A.V., Garmire L.X., Greene C.S., Goff L.A., Li Y., Ngom A., Ochs M.F. Enter the Matrix: Factorization Uncovers Knowledge from Omics. Trends Genet. 2018;34:790–805. doi: 10.1016/j.tig.2018.07.003. [DOI] [PMC free article] [PubMed] [Google Scholar]
31.Wang K., Gaitsch H., Poon H., Cox N.J., Rzhetsky A. Classification of common human diseases derived from shared genetic and environmental determinants. Nat. Genet. 2017;49:1319–1325. doi: 10.1038/ng.3931. [DOI] [PMC free article] [PubMed] [Google Scholar]
32.Denny J.C., Ritchie M.D., Basford M.A., Pulley J.M., Bastarache L., Brown-Gentry K., Wang D., Masys D.R., Roden D.M., Crawford D.C. PheWAS: demonstrating the feasibility of a phenome-wide scan to discover gene-disease associations. Bioinformatics. 2010;26:1205–1210. doi: 10.1093/bioinformatics/btq126. [DOI] [PMC free article] [PubMed] [Google Scholar]
33.Gratten J., Visscher P.M. Genetic pleiotropy in complex traits and diseases: implications for genomic medicine. Genome Med. 2016;8:78. doi: 10.1186/s13073-016-0332-x. [DOI] [PMC free article] [PubMed] [Google Scholar]
34.Yang J., Bakshi A., Zhu Z., Hemani G., Vinkhuyzen A.A., Lee S.H., Robinson M.R., Perry J.R., Nolte I.M., van Vliet-Ostaptchouk J.V., LifeLines Cohort Study Genetic variance estimation with imputed variants finds negligible missing heritability for human height and body mass index. Nat. Genet. 2015;47:1114–1120. doi: 10.1038/ng.3390. [DOI] [PMC free article] [PubMed] [Google Scholar]
35.Denny J.C., Rutter J.L., Goldstein D.B., Philippakis A., Smoller J.W., Jenkins G., Dishman E., All of Us Research Program Investigators The “All of Us” Research Program. N. Engl. J. Med. 2019;381:668–676. doi: 10.1056/NEJMsr1809937. [DOI] [PMC free article] [PubMed] [Google Scholar]
36.Pickrell J.K., Berisa T., Liu J.Z., Ségurel L., Tung J.Y., Hinds D.A. Detection and interpretation of shared genetic influences on 42 human traits. Nat. Genet. 2016;48:709–717. doi: 10.1038/ng.3570. [DOI] [PMC free article] [PubMed] [Google Scholar]
37.Banda Y., Kvale M.N., Hoffmann T.J., Hesselson S.E., Ranatunga D., Tang H., Sabatti C., Croen L.A., Dispensa B.P., Henderson M. Characterizing Race/Ethnicity and Genetic Ancestry for 100,000 Subjects in the Genetic Epidemiology Research on Adult Health and Aging (GERA) Cohort. Genetics. 2015;200:1285–1295. doi: 10.1534/genetics.115.178616. [DOI] [PMC free article] [PubMed] [Google Scholar]
38.Carey D.J., Fetterolf S.N., Davis F.D., Faucett W.A., Kirchner H.L., Mirshahi U., Murray M.F., Smelser D.T., Gerhard G.S., Ledbetter D.H. The Geisinger MyCode community health initiative: an electronic health record-linked biobank for precision medicine research. Genet. Med. 2016;18:906–913. doi: 10.1038/gim.2015.187. [DOI] [PMC free article] [PubMed] [Google Scholar]
39.Bycroft C., Freeman C., Petkova D., Band G., Elliott L.T., Sharp K., Motyer A., Vukcevic D., Delaneau O., O’Connell J. The UK Biobank resource with deep phenotyping and genomic data. Nature. 2018;562:203–209. doi: 10.1038/s41586-018-0579-z. [DOI] [PMC free article] [PubMed] [Google Scholar]
40.Nagai A., Hirata M., Kamatani Y., Muto K., Matsuda K., Kiyohara Y., Ninomiya T., Tamakoshi A., Yamagata Z., Mushiroda T., BioBank Japan Cooperative Hospital Group Overview of the BioBank Japan Project: Study design and profile. J. Epidemiol. 2017;27(3S):S2–S8. doi: 10.1016/j.je.2016.12.005. [DOI] [PMC free article] [PubMed] [Google Scholar]
41.Ritchie M.D., Denny J.C., Crawford D.C., Ramirez A.H., Weiner J.B., Pulley J.M., Basford M.A., Brown-Gentry K., Balser J.R., Masys D.R. Robust replication of genotype-phenotype associations across multiple diseases in an electronic medical record. Am. J. Hum. Genet. 2010;86:560–572. doi: 10.1016/j.ajhg.2010.03.003. [DOI] [PMC free article] [PubMed] [Google Scholar]
42.Gaziano J.M., Concato J., Brophy M., Fiore L., Pyarajan S., Breeling J., Whitbourne S., Deen J., Shannon C., Humphries D. Million Veteran Program: A mega-biobank to study genetic influences on health and disease. J. Clin. Epidemiol. 2016;70:214–223. doi: 10.1016/j.jclinepi.2015.09.016. [DOI] [PubMed] [Google Scholar]
43.Martin A.R., Kanai M., Kamatani Y., Okada Y., Neale B.M., Daly M.J. Clinical use of current polygenic risk scores may exacerbate health disparities. Nat. Genet. 2019;51:584–591. doi: 10.1038/s41588-019-0379-x. [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Document S1. Figures S1–S12 and Tables S1, S3, and S5–S22

mmc1.pdf^{(1.9MB, pdf)}

Data S1. Tables S2 and S4

mmc2.xlsx^{(62.3KB, xlsx)}

Document S2. Article plus Supplemental Information

mmc3.pdf^{(3MB, pdf)}

PERMALINK

Pleiotropy-Based Decomposition of Genetic Risk Scores: Association and Interaction Analysis for Type 2 Diabetes and CAD

Daniel I Chasman

Franco Giulianini

Olga V Demler

Miriam S Udler

Abstract

Introduction

Subjects and Methods

Cohorts

The Women’s Genome Health Study (WGHS)

JUPITER Study

Genotyping

Pleitoropy-Based Decomposition of a Disease GRS

SNP Selection and Coding

Scaling, Imputation, and Zeroing of Beta Coefficients

GRS Decomposition

Additional Scaling of Beta Coefficients

Other Statistical Analysis

Ethics Statement

Results

Overview

Decomposition of the GRS for T2D

Figure 1.

Figure 2.

Associations and Interactions of the Decomposed T2D GRS

Table 1.

Decomposition of the GRS for CAD

Figure 3.

Figure 4.

Associations and Interactions of the Decomposed CAD GRS

Table 2.

Table 3.

Discussion

Acknowledgments

Footnotes

Appendix A: Correlation of Component GRSs and Standardization of Beta Coefficients

Web Resources

Declaration of Interests

Supplemental Data

References

Associated Data

Supplementary Materials

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases