Score statistics and their variance-covariance matrix.

We assume we have data from L GWA studies, with each consisting of n^(l) subjects, l = 1,⋯,L. Each gene in that pathway can contain one or multiple SNP(s), while any two genes may have some overlapped SNPs. For simplicity, we use superscript l to represent an individual study. For subject i in study l, i = 1,⋯,n^(l), let be the dichotomous outcome (e.g., disease condition, case/control status) taking values from {0,1}, and let be the vector of covariates to be adjusted for. The centralized genotypes of q SNPs within a pathway are presented as a vector for subject i. We assume the following logistic regression model as the risk model

Under the self-contained null hypothesis H₀: γ = 0, we denote the maximum likelihood estimate of α^(l) as . Let and . The Rao’s score statistic vector on γ, which is the sum of score vectors from L participating studies, follows the asymptotic multivariate normal distribution N(0,V), where (1) and (2)

For study l, let be the number of subjects having their genotypes measured as (or imputed) at SNP t, where . As pointed out by Hu, Berndt (24) if the covariates and genotypes are uncorrelated or weakly correlated, the covariance between scores at SNPs t and S can be approximated as (3) where is the number of samples that have their genotypes available at both SNPs in study l, , and . Here, we assume that the Pearson’s correlation coefficient ρ_ts between two SNPs is the same among all participating studies. This assumption is valid as long as subjects from all studies are sampled from the same source population, or the population under study is relatively homogeneous, such as a study of subjects with European ancestry in the United States.

When only the summary statistics, i.e., the estimated marginal log odds ratios and their standard errors are available for each of the L studies, the score statistic at SNP t, defined by (Eq 1) can be approximated as (4)

Note that , thus according to (Eq 3), we have (5)

Assume that ρ_ts can be estimated from a public dataset (e.g., 1000 Genomes Project) and the sample sizes and are known, we can approximately recover the variance-covariance matrix V = (V_ts)_q×q of score statistics S = (S_t)_q×1. In cases when we only have the SNP p-value p and its marginal log odds ratio , we can compute its standard error as , where is the quantile satisfying , with representing a 1-df chi-squared random variable.

Combining score statistics for pathway analysis.

With recovered score statistics vector S and its variance-covariance matrix V, we can conduct a pathway association test using the framework of the ARTP method. The ARTP method first combines p-values of individual SNPs within a gene to form a gene-based association statistic (i.e., the gene-level p-value), and then combines the gene-level p-values into a final testing statistic for the pathway-outcome association. In the original ARTP method, [9] proposed the use of a resampling-based method to evaluate the significance level of the pathway association test. Here we integrate the SNP-level score statistics into the ARTP framework and use DSA [28] to evaluate the significance level, which is much faster than the original ARTP algorithm [31]. Below is a brief summary of the improved ARTP algorithm.

First we obtain the p-values of q_j distinct SNPs in gene j as . Let be their order statistics such that . For any predefined integer K and SNP-level cut points c₁ < ⋯ < c_k, we define the observed negative log product statistics for that gene at cut point c_k as

We sample M copies of vectors of the score statistic from the null distribution N(0,V) and convert each of them to be the tail probability of as , m = 1,⋯,M, which are then used to calculate , m = 1,⋯,M. The significance of can be estimated as

The ARTP statistic for testing association between gene j and the outcome is defined as . Note that for any , the set forms its empirical null distribution. The significance of therefore can be estimated as

This idea, which was given by [32], can be used to avoid the computationally challenging nested two-layer resampling procedure for evaluating p-values. The p-value of can be readily calculated as where . is the estimated gene-level p-value for the association between the outcome and the jth gene. To obtain the pathway p-value, a similar procedure as above can be applied to combine already established gene-level p-values , j = 1,⋯,J, through a set of K′ gene-level cut points d₁ < ⋯ <d_K′. For simplicity, let be the significance (p-value) of negative log product statistics defined on , j = 1,⋯,J at a specific cut point d_k, k = 1,⋯,K′. The ARTP statistic for the pathway association is defined as . The top genes, at which , can be regarded as the set of selected candidate genes that collectively convey the strongest pathway association signal.

In the following discussion, we will use the term sARTP to represent the proposed pathway analysis procedure using the SNP-level summary statistics as input, and reserve the term ARTP to represent the original ARTP procedure that requires the individual-level genetic data. Both procedures adopt the DSA algorithm to accelerate evaluating the significance level. When performing the pathway analysis in this paper, we set SNP-level cut points as (c₁,c₂) = (1,2), i.e., gene-level association is summarized by one or two most significant SNPs within each gene, and gene-level cut points as , k = 1,⋯,10, where J is the number of genes in a pathway, and is the largest integer that is less or equal to J / 20. We used M = 10⁵ DSA steps to assess the significance level of each pathway in the initial screening. For pathways with estimated p-values <10⁻⁴, we further refined their p-value estimates with M = 10⁷ or 10⁸ DSA steps.

Applying sARTP to meta-analysis result.

Many GWAS consortia usually publish their meta-analysis results by providing only the combined results from the fixed effects model, rather than the summary statistics from each participating study. We can apply sARTP to this meta-analysis result directly, with some modifications. First, since the reported marginal log odds ratios for each SNP by using the fixed effects inverse-variance weighting method is given by with its standard error given by (6)

Based on (Eq 4), we can see . By assuming large sample sizes and certain conditions (see S1 Text), we can also approximate the covariance between S_t and S_s, which is given by (Eq 5), as (7) where , and . Thus, using just the meta-analysis result, without knowing summary statistics from each participating study, we can still obtain S_t exactly, and approximately recover V_ts. As a result, we can carry out the pathway-based meta-analysis based on the SNP-level meta-analysis result as if it were summary data from a single study. We call this approach the Meta-analysis based sARTP (MsARTP).

However, to apply the MsARTP, we need additional sample size information n_t and n_st in order to properly estimate the variance-covariance matrix defined by (Eq 7). If the same set of SNPs are studied by all participating studies, we have n_t = n_s = n_st, and the approximation (Eq 7) becomes , i.e., we can obtain the estimated variance-covariance matrix without knowing n_t and n_st. But in most applications, not all GWAS choose the same SNP genotyping array, even after the imputation using the same reference genomes. As a result, the SNP coverage, i.e., the set of SNPs evaluated in each participating study can be quite different. In those situations, we need to know the SNP coverage information in each participating study in order to obtain n_s and n_st. We will show in the Results Section that using MsARTP with an inappropriate uniform coverage assumption (i.e., n_ts = n_t = n_s), which is commonly made by many multi-locus approaches, can lead to inflated type I error.

Given SNP-level summary statistics from each participating study, we can either apply sARTP directly, or first conduct a SNP-level meta-analysis, and then apply MsARTP to the meta-analysis result. These two approaches use the same score statistics, and different but consistent estimates for the variance-covariance matrix. Numeric experiments in the Results Section suggest that these two approaches generate vary similar pathway p-values.

Pathway and gene definition.

We downloaded definitions for 4,716 human and murine (mammalian) pathways (gene sets) from the MSigDB v5.0 (C2: curated gene sets). Genomic definitions for genes were downloaded from Homo sapiens genes NCBI36 and reference genome GRCh37.p13 using the Ensemble BioMart tool.

DIAGRAM study.

The DIAGRAM (DIAbetes Genetics Replication And Meta-analysis) consortium conducted a large-scale GWAS meta-analysis to characterize the genetic architecture of T2D [13]. We downloaded the summary statistics generated by the DIAGRAMv3 (Stage 1) GWAS meta-analysis from www.diagram-consortium.org [13]. The meta-analysis studied 12 GWAS with European ancestry consisting of 12,171 cases and 56,862 controls. Up to 2.5 million autosomal SNPs with minor allele frequencies (MAFs) larger than 1% were imputed using CEU samples from Phase II of the International HapMap Project. Study-specific covariates were adjusted in testing T2D-SNP association under an additive logistic regression model [13]. SNP-level summary statistics from each GWAS were first adjusted for residual population structure using the genomic control (GC) method [33], and then combined in the fixed effects meta-analysis.

We sorted 2.5 million autosomal SNPs by their corresponding meta-analysis sample sizes in S1 Fig, which shows that there are two major groups of SNPs with equal sample sizes. One group of 469,985 SNPs (19.0%) had 12,171 cases and 56,862 controls, which included all the available samples in the meta-analysis; another group of 1,431,361 SNPs (57.9%) had 9,580 cases and 53,810 controls. Since the calculation of covariance V_ts in (Eq 7) relies on n_ts, the number of samples having genotypes available at both SNP s and SNP t, in order to obtain an accurate estimate of n_ts, we focused on these two groups of SNPs, which in combination had a total of 1,901,346 SNPs. For any two SNPs in this reduced set, it is certain n_ts = min(n_t, n_s). The Pearson’s correlation coefficients ρ_ts were estimated using an external reference panel consisting of genotypes on 503 European subjects (CEU, TSI, FIN, GBR, and IBS) from the 1000 Genomes Project (Phase 3, v5, 2013/05/02).

GERA study.

We assembled a GWAS on T2D from the Genetic Epidemiology Research on Adult Health and Aging (GERA, dbGaP Study Accession: phs000674.v1.p1). The GERA project includes a cohort of over 100,000 adults who are members of the Kaiser Permanente Medical Care Plan, Northern California Region, and participating in the Kaiser Permanente Research Program on Genes, Environment, and Health (RPGEH). From the GERA data, we compiled a GWAS with 7,638 T2D cases and 54,319 controls (subjects without T2D) who self-reported to be non-Hispanic White Europeans in the RPGEH survey. We performed the genotype imputation with IMPUTE2 [34] using CEU reference samples from Phase II of the International HapMap Project. After removing SNPs with low imputation quality (r² < 0.3), we ended up with 2.4 million SNPs for further analysis. In the single-locus analysis, we adjusted for the categorized body mass index (BMI) provided in the downloaded dataset (adding a category for missing BMI), gender, year of birth (in five-year categories), a binary indicator on whether or not a participant was diagnosed with cancer (includes malignant tumors, neoplasms, lymphoma and sarcoma), and the top five eigenvectors for the adjustment of population stratification. In the following discussion, we refer this assembled T2D GWAS as the GERA study.

When analyzing the SNP-level summary data from the GERA study, the Pearson’s correlation coefficients ρ_ts were estimated using an external reference panel consisting of genotypes on 503 European subjects from the 1000 Genomes Project.

AGEN-T2D study.

The Asian Genetic Epidemiology Network (AGEN) consortium carried out a meta-analysis by combining eight GWAS of T2D with a total of 6,952 cases and 11,865 controls from eastern Asian populations [10]. The meta-analysis was conducted with the fixed effect model. We obtained SNP-level summary statistics on 2.6 million imputed and genotyped autosomal SNPs from AGEN, and used this summary data to evaluate whether pathway associations identified in European populations remain to be present in Asians. We adopted an external reference panel consisting of 312 eastern Asian subjects (103 from CHB, 105 from CHS, and 104 from JPT) from the 1000 Genomes Project for the variance-covariance matrix estimation in the pathway analysis.

Empirical size of sARTP and MsARTP.

To evaluate the empirical size of sARTP and MsARTP, we conducted a simulation study by using individual-level GWAS data of the pathway PUJANA_BREAST_CANCER_WITH_BRCA1_MUTATED_UP (including 728 SNPs in 50 genes) from the GERA study. We picked 12,000 samples randomly for this experiment. By keeping their genotypes unchanged, we randomly assigned 6,000 subjects as cases and the remaining as controls to generate 500,000 datasets. We split each dataset into three case-control studies, each with 2,000 cases and 2,000 controls. To mimic the scenario when not all studies have their genotypes measured on the same set of SNPs (such as the one occurred in the DIAGRAM and AGEN data), we assumed that each case-control study had genotypes measured on only half of SNPs in the pathway. For each generated dataset that consisted of three case-control studies, we applied sARTP to the SNP-level summary data obtained from each case-control study, and MsARTP to the meta-analysis result based on the three case-control studies, with the variance-covariance matrix estimated by an external reference panel (with 503 European reference samples from the 1000 Genomes Project), or an internal reference panel (with 500 samples randomly selected from the GERA data).

Based on results from the 500,000 generated datasets, this simulation study showed that both sARTP and MsARTP, using the internal or external reference samples, can well control their empirical sizes (Table 1). Given the same reference panel, the p-values estimated from sARTP and MsARTP are highly consistent (Pearson’s correlation coefficient > 0.99). Furthermore, the p-values of sARTP (or MsARTP) estimated with an external or internal reference panel are also very consistent (Pearson’s correlation coefficient > 0.99). More numeric experiments demonstrating the validity of sARTP under the null are described in S2 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Empirical sizes of the sARTP, MsARTP, and MsARTP-u procedures.

https://doi.org/10.1371/journal.pgen.1006122.t001

To demonstrate the importance of knowing n_t and n_ts when applying MsARTP to the meta-analysis result, we analyzed each simulated dataset using MsARTP assuming the uniform coverage (n_ts = n_t = n_s). We called this approach MsARTP-u. It is clear from Table 1 that MsARTP-u assuming the uniform coverage suffers from inflated type I errors with either the internal or external reference panel.

Empirical power of sARTP and MsARTP for pathway-based meta-analysis.

We conducted a set of simulation studies to compare the power of different strategies to carry out pathway analysis when SNP-level summary statistics were available from multiple studies. We considered a hypothetical pathway consisting of 50 genes randomly selected from chromosome 17, each with 20 randomly chosen SNPs. The joint genotype distribution at the 20 SNPs within each gene was defined by the observed genotypes in the GERA study. We further assumed that all genes in that pathway are independent. This assumption is unnecessary for sARTP and MsARTP, but it was introduced for simplifying the simulation. For the risk model, we assumed the first () genes were associated with the outcome. Within each outcome-associated gene, we picked the SNP with its MAF closest to the median MAF level within the gene to be functional. We considered the following risk model (8) where is the genotype (encoded as 0, 1, or 2 according to counts of minor alleles) at the functional SNP within gene l. Under this model, is also the marginal log odds ratio for the lth functional SNP [9]. Given the sample sizes of cases and controls, and the MAF of the lth functional SNP, was chosen such that the theoretical power of the trend test to detect the lth functional SNP is equal to (), with 0.05 as the targeted type I error rate. For every pair of , we generated 1,000 datasets, each consisting of three case-control studies, with the same sample size and SNP coverage configurations used for evaluating the empirical size. Given the genotype distribution in the general population, individual-level genotype data for a case-control study can be generated according to the assumed risk (model 8).

We assumed that only SNP-level summary statistics from each of the three studies were available. For each simulated dataset, we applied sARTP and MsARTP, using either an internal or external reference panel to estimate the variance-covariance matrix. The sARTP and MsARTP approaches integrate association evidence across SNP-level summary statistics, which are obtained by pooling information from all participating studies on individual SNPs. As a comparison, we also considered a naïve approach, in which we first applied sARTP to analyze the summary statistics from each study separately, and then combined the three pathway p-values with Fisher’s method. This naïve approach could be useful when the researchers do not have access to the SNP-level summary data but the pathway p-values from individual studies. The empirical powers are compared at the type I error level of 0.05, and are summarized in Table 2. It is obvious that the pathway-based meta-analysis using sARTP, with either the internal or external reference panel, have almost the same level of power as the MsARTP method. It is also evident that both sARTP and MsARTP are more powerful than the naïve approach, which suggests that it is always be beneficial to have the SNP-level summary statistics from each participating study, or SNP-level meta-analysis result when conducting a pathway analysis.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Power comparisons under the type I error rate of 0.05 when analyzing data from three studies.

https://doi.org/10.1371/journal.pgen.1006122.t002

Given the SNP coverage information, the MsARTP method is a valid pathway association test that has well controlled type I error and similar power to the sARTP method. In the following analysis, either sARTP or MsARTP is chosen depending on the type of available data. For the sake of simplicity, we always label the chosen procedure as sARTP.

Power comparison between sARTP and aSPUsPath.

Since the aSPUsPath method in the current aSPU package cannot handle summary data from multiple studies, or meta-analysis results from studies with varied SNP coverage, we focused on the scenario with just one study, and adopted the similar simulation strategy as the one used by [8] to compare the power between sARTP and aSPUsPath. We simulated haplotypes on a set of SNPs within a gene in the general population using the algorithm of Wang and Elston [35]. Then the joint genotypes on a subject can be formed by randomly pairing two haplotypes. In brief, we first chose the MAF for each SNP by randomly sampling a value from the uniform distribution U(0.1,0.4). Then for the set of SNPs in a gene we sampled a latent vector from a multivariate normal distribution with a covariance matrix Cov(z_i, z_j) = ρ^|i−j|,1≤i,j≤q, where ρ was sampled from the uniform distribution U(0,0.8) for a given gene. We randomly picked 50% of the SNPs and converted their simulated z_i into minor and major alleles (coded as 0, 1), with the cuts chosen for each z_i such that the resultant minor allele has its frequency defined by the specified MAF. For the remaining SNPs, we used the same algorithm to dichotomize −z_i into minor and major alleles. This created a more realistic haplotype structure such that a haplotype can consist of a mixture of minor and major alleles. Genotypes on SNPs from different genes were generated independently.

Given the number of genes (20, 50, or 80) in a pathway, the proportion of genes (5%, 10%, 20%, and 30%) associated with the outcome, and a chosen common value for all log odds ratios (γ*) in the risk (model 8), we repeated the following steps to generate 1,000 case-control studies, with each consisting of 1,000 cases and 1,000 controls. First, the number of SNPs within each gene was randomly chosen from 10 to 100. Second, for each randomly selected outcome-associated gene, we randomly picked a functional SNP. Third, we use the aforementioned algorithm of Wang and Elston [35] to generate the individual-level genotype data for a case-control study according to the specified risk model. We also considered the situation where all γ* in the risk (model 8) had the same magnitude but different directions. More precisely, when generating a case-control study at the third step, we defined the risk (model 8) by randomly choosing the direction of each log odds ratio to be positive or negative with equal probability. Furthermore, we considered a more complex scenario where each outcome-associated gene had one or two functional SNPs, each with equal probability.

All simulation results are given in S1 and S2 Tables. It is clear that sARTP are generally more powerful than aSPUsPath, especially when the signal-to-noise ratio (the proportion of genes including a functional SNP) is relatively low. The two types of tests tend to have comparable performance when the signal-to-noise ratio increases to 30%, although it is uncommon for a candidate pathway to have such a high signal-to-noise ratio in real applications. For example, among the 4,713 candidate pathways analyzed in the next section, only 4.2% and 0.9% of the pathways have over 20% and 30% of their genes that are likely to contain association signals (i.e., with gene-level p-values < 0.05).

From S1 and S2 Tables, we also notice that the advantage of sARTP over aSPUsPath is more evident if not all minor alleles of the functional SNPs are deleterious (or protective) variants (i.e., γ* in the risk (model 8) are not all positive). This is expected, as the sARTP approach does not take the effect direction of the minor allele at each SNP into consideration, while aSPUsPath integrates a set of candidate statistics, including the one similar to the burden test that assumes all minor alleles are either deleterious or protective. When this assumption is not valid, the inclusion of the burden test statistic in aSPUsPath is unlikely to enhance the power, but certainly would increase the multiple-testing penalty.

Findings from the European populations.

Since our goal was to identify new susceptibility loci for T2D through the pathway analysis, we excluded 170 high evidence T2D associated SNPs that were either listed in [13] or found from the GWAS Catalog satisfying the following three conditions simultaneously: (1) were investigated by GWAS of samples with European ancestry; (2) had reported p-values <10⁻⁷ on the initial study; and (3) were replicated on independent studies. We excluded 195 SNPs that has their single-locus testing p-values less than 10⁻⁷ in either DIAGRAM or GERA data to ensure that the pathway analysis result was not driven by a single SNP. In addition, we further excluded genes within a ±500kb region from each of the removed SNPs to eliminate potential association signals that could be caused by linkage disequilibrium (LD) with the index SNPs.

We conducted three types of pathway-based meta-analyses using sARTP, including the one using the DIAGRAM SNP-level summary statistics, the one using the GERA SNP-level summary statistics, and the pathway meta-analysis combining SNP-level summary statistics from both DIAGRAM and GERA studies. When applying the pathway-based meta-analysis to a single gene, we refer to this as the gene-level meta-analysis. We used the external reference panel of 503 Europeans from the 1000 Genomes Project to estimate the variance-covariance matrix.

Before performing a pathway analysis, we applied LD filtering to remove redundant SNPs. For any two SNPs with their pairwise squared Pearson’s correlation coefficient > 0.9 estimated from the external reference panel from the 1000 Genomes Project, we removed the one with a smaller value defined as, 2f(1−f)n₀n₁n⁻¹, where n₀ and n₁ are numbers of controls and cases, n = n₀ + n₁, and f is the MAF based on the reference panel. This value is proportional to the non-centrality parameter of the trend test statistic at a given SNP. We also excluded SNPs with MAF < 1%. After all SNP filtering steps, we had a total of 4,713 pathways for the analysis. The summary of the number of genes and SNPs used in each pathway analysis is given in S2 Fig.

The DIAGRAM study had a genomic control inflation factor λ_GC = 1.10 based on the published meta-analysis result. The assembled GERA T2D GWAS had λ_GC = 1.08. When conducting the pathway analysis on each of two studies, we adjusted the inflation by using the corresponding to rescale the standard error of estimated log odds ratio at each SNP. The single-locus meta-analysis combining results from DIAGRAM and GERA datasets had an inflation factor λ_GC = 1.067 after each study had adjusted for its own inflation factor. We further adjusted this inflation in the pathway and gene-level meta-analysis when combining SNP-level summary statistics from both studies using formulas (4) and (5).

The Q-Q plots of gene-level and pathway p-values are given in Fig 2. Gene-level p-value Q-Q plots based on the three analyses show no sign of inflation with their λ_GC close to 1.0, but suggest that there are enriched gene-level association signals at the tail end. The pathway p-value Q-Q plots, on the other hand, shift away from the diagonal identify line and have much higher λ_GC, which suggests that T2D associated genes are preferably included in pathways under study. In fact, it can be seen from S3 Fig that a gene with a smaller gene-level meta-analysis p-value tends to be included in more pathways, even though the 4,713 pathways collected from MSigDB v5.0 are not specifically defined for the study of T2D.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Q-Q plots of gene-level and pathway p-values based on the sARTP procedure on the DIAGRAM study, the GERA study, and the two studies combined.

(Left) Q-Q plots of gene-level p-values on 15,946 genes based on the sARTP gene-based analysis of the DIAGRAM study (DIAGRAM), the GERA study (GERA), and the two studies combined (META). (Right) Q-Q plots of pathway p-values on 4,713 pathways based on the sARTP pathway analysis of the DIAGRAM study (DIAGRAM), the GERA study (GERA), and the two studies combined (META).

https://doi.org/10.1371/journal.pgen.1006122.g002

Fig 2 illustrates that the gene and pathway level signal from the GERA study tends to be slightly stronger than that from the DIAGRAM study. The main reason is that the DIAGRAM summary result had gone through two rounds of inflation adjustments, with the first round done at each participating study, and the second round on the meta-analysis result. Also, its second round adjustment (λ_GC = 1.10) is larger than the one applied to the GERA study (λ_GC = 1.08). Adjusting for λ_GC in the pathway analysis could be too conservative, since some proportion of the inflation can be caused by the real polygenic effect. A less conservative adjustment could be possible, but it might not be adequate. More discussions on this issue are given in the Discussion Section.

Based on the pathway meta-analysis on a total of 4,713 pathways, we identified 43 significant pathways with p-values less than 1.06×10⁻⁵, the family-wise significant threshold based on the Bonferroni correction. Their pathway meta-analysis results as well as results from individual studies are summarized in Table 3. More detailed results on each of 43 significant pathways are given in the S6–S48 Figs and S6 Table. There are a total of 15,946 unique genes in all 4,713 pathways. The top 50 genes with smallest gene-level p-values based on the gene meta-analysis are listed in S3 Table. Because of the LD filtering, a gene belonging to two pathways might end up with slightly different sets of SNPs. To remove this ambiguity, we obtained the gene-level p-values by conducting a gene-level meta-analysis on each gene separately.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Summary of 43 significant pathways detected by the pathway meta-analysis based on the DIAGRAM and GERA studies.

https://doi.org/10.1371/journal.pgen.1006122.t003

From Table 3, we can notice that some identified pathways have relatively weak association signals from each of the two studies, but have very significant p-values based on the pathway meta-analysis on the two studies combined. For example, the pathway RIZ_ERYTHROID_DIFFERENTIATION has p-values of 0.0233 and 0.0231 based on DIAGRAM and GERA studies, respectively. Combining these two p-values using Fisher’s method yields a p-value of 0.0046. On the other hand, the pathway meta-analysis produces a much more significant result (p = 6.15×10⁻⁷). This demonstrates the power advantage of the pathway meta-analysis over the approach that simply combines the pathways p-values from individual studies. The aforementioned simulation studies also confirmed this observation (Table 2).

In Fig 3, we illustrate the connection between the 43 significant pathways and a group of genes showing association evidence. For the purpose of illustration, in the figure we only focus on 46 genes that are covered by the 43 pathways and have their gene-level meta-analysis p-values less than 0.001. It is evident from Fig 3 that a cluster of 4 genes, UBE2Z, SNF8, GIP, and ATP5G1, has the most significant gene-level p-values (S3 Table), and contribute association signals to 20 out of 43 significant pathways (S6–S25 Figs). These 4 genes overlap each other at chromosome 17q21. This region contains a previously unidentified genome-wide significant synonymous SNP rs1058018 (meta-analysis p = 3.06×10⁻⁸) after two rounds of inflation adjustments. More detailed information on SNP rs1058018 and SNPs in that region are given in S4 Table, S4 and S5 Figs. By conditioning on rs1058018, none of the other SNPs in this region are significant based on the conditional association analysis using the GERA individual-level GWAS data. Based on GTEx data v6, rs1058018 is a cis eQTL for UBE2Z in blood (p = 7.9×10⁻¹⁵). UBE2Z is involved in Class I MHC antigen processing and presentation (GeneCards). The region at 17q21 was previously implicated to be associated with T2D through a candidate gene/loci approach [36]. Although genes at the 17q21 region carry the strongest association signal, 11 out of those 20 pathways remain to be globally significant (p < 1.06×10⁻⁵) after excluding those genes from the pathway definition (Table 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Heat map of gene-level p-values on selected genes within 43 significant pathways based on the DIAGRAM and GERA studies.

There are 46 unique genes in the 43 significant pathways that have their gene-level meta-analysis p-values less than 0.001. Each row in the plot represents one of 43 significant pathways. Each column represents one of the 46 unique genes. The chromosome IDs of 46 unique genes are given in parentheses. The color of each cell represents the gene-level p-value (in the −log₁₀ scale). A cell for a gene that is not included in a pathway is colored gray in the corresponding entry. The orders of genes (x-axis) and pathways (y-axis) are arranged according to their gene and pathway meta-analysis p-values.

https://doi.org/10.1371/journal.pgen.1006122.g003

The majority of 43 identified pathways are enriched with signals from multiple chromosomal regions as demonstrated by the Q-Q plots of their SNP-level and gene-level p-values (S6–S48 Figs). For example, the strongest T2D-associated pathway, SCHLOSSER_SERUM_RESPONSE_UP, consists of 103 genes, which includes two genes with p-values < 0.001, and has 20 genes with p-values between 0.001 and 0.05 (S26 Fig, and Supplemental data). We conducted the ingenuity pathway analysis on those 22 genes with p-values less than 0.05, and found enrichment of these genes in caveolae-mediated cytosis (important for removal of low/high density lipoproteins), and lipid metabolism pathways, and in functions/diseases related to differentiation of phagocytes and transport of proteins.

It is assuring that our pathway analysis detected several pathways that are natural candidates underlying the development of T2D, including the pathways KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG (S31 Fig), KEGG_TYPE_II_DIABETES_MELLITUS (S41 Fig), KEGG_INSULIN_SIGNALING_PATHWAY (S43 Fig), and REACTOME_REGULATION_OF_BETA_CELL_DEVELOPMENT (S33 Fig). It is worth emphasizing that these pathways were analyzed after excluding genes in the neighborhood of 170 GWAS established T2D loci and 195 SNPs with p-values <10⁻⁷ on either DIAGRAM or GERA data, which suggests that these well-defined T2D-related pathways are enriched with additional unidentified and contributory T2D-associated genes.

Among the 43 globally significant pathways, there are multiple ones that are defined according to specific gene expression patterns on various tumor types, including pancreatic adenocarcinoma (S21 Fig), hepatocellular carcinoma (HCC) (S27 and S32 Figs), bladder carcinoma (S42 Fig), nasopharyngeal carcinoma (S24 Fig), and familial breast cancer (S34 and S37 Figs). It is well recognized that T2D patients have elevated risk of cancer at multiple cancer sites, such as the liver and pancreas [37, 38]. These findings can provide valuable insights into the genetic basis underlying the connection between T2D and a host of different cancers.

In the above analysis, we used the sARTP method with the gene-level association evidence summarized by one or two most significant SNPs within each gene, under the assumption that there are at most two independent association signals within a given gene. We also applied sARTP by using 3 SNP-level cut points (i.e., (c₁,c₂,c₃) = (1,2,3)) to reanalyze the 4,713 pathways based on the combined data of DIAGRAM and GERA. It appears that results obtained by sARTP with 3 SNP-level cut points are very consistent with those with 2 cut points (S49 Fig).

Findings from Eastern Asian populations.

We reanalyzed the 43 significant pathways identified from the European populations using summary-level data generated by the AGEN-T2D study. An inflation factor λ_GC = 1.03 calculated from the AGEN-T2D meta-analysis was adjusted in the pathway meta-analysis. The genetic regions excluded from analyzing the DIAGRAM and GERA studies were also excluded from the AGEN study. The results were summarized in Table 4. There are 10 out of 43 pathways with the unadjusted p-value less than 0.05, suggesting that many pathways identified from the European populations were also enriched with T2D-associated genes in the eastern Asian populations (S7 Table). Among the 43 pathways, we were able to identify 4 significant T2D-associated pathways at the false discovery rate (FDR [39]) of 0.05 (S27, S12, S32 and S36 Figs), and 3 additional T2D-associated pathways at the FDR of 0.1 (S47, S44, and S25 Figs). All the pathway p-values remain basically the same level if we further excluded genes within ±500kb regions surrounding the GWAS T2D loci established in eastern Asian populations. These results support the presence of trans-ethnic pathway effect on T2D in European and eastern Asian populations [11, 12].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Pathway p-values and FDR adjusted p-values based on the AGEN-T2D study.

https://doi.org/10.1371/journal.pgen.1006122.t004

Given the existing epidemiologic evidence on the close connection between T2D and the liver cancer, it is noteworthy that the two HCC related pathways (S27 and S32 Figs) identified in European populations remain to be significant in eastern Asian populations at the FDR of 0.05 (Table 4). The pathway PATIL_LIVER_CANCER consists of 653 genes (after data preprocessing) that are highly expressed in HCC and are enriched with genes having functions related to cell growth, cell cycle, metabolism, and cell proliferation [40]. The other pathway, HOSHIDA_LIVER_CANCER_SUBCLASS_S3 consists of 240 genes that show similar gene expression variation patterns and together define a HCC subtype with its unique histologic, molecular and clinical characteristics [41]. These two pathways have only 6 genes in common, and none of the 6 genes has a gene-level p-value < 0.05 in either European or eastern Asian data. More in depth investigations of these two complementary pathways could lead to further understanding the connection between T2D and the liver cancer.

The genome-wide significant SNP rs1058018 at the 17q21 region identified through the combined analysis of DIAGRAM and GERA studies turned out to be null in the AGEN-T2D study (p = 0.29). This could be due to the relatively small sample size of the AGEN-T2D study, or the genetic risk heterogeneity at the 17q21 locus among different ethnic populations. Nevertheless, 2 out of the 20 pathways (S12 and S25 Figs) that contain genes within the 17q21 region are still significant at the FDR of 0.1. Among the 23 pathways that do not contain any gene within the 17q21 region, 5 pathways remain significant at the FDR of 0.1 (S27, S32, S36, S47 and S44 Figs).