Abstract

In addition to methods that can identify common variants associated with susceptibility to common diseases, there has been increasing interest in approaches that can identify rare genetic variants. We use the simulated data provided to the participants of Genetic Analysis Workshop 17 (GAW17) to identify both rare and common single-nucleotide polymorphisms and pathways associated with disease status. We apply a rare variant collapsing approach and the usual association tests for common variants to identify candidates for further analysis using pathway-based and tree-based ensemble approaches. We use the mean log p-value approach to identify a top set of pathways and compare it to those used in simulation of GAW17 dataset. We conclude that the mean log p-value approach is able to identify those pathways in the top list and also related pathways. We also use the stochastic gradient boosting approach for the selected subset of single-nucleotide polymorphisms. When compared the result of this tree-based method with the list of single-nucleotide polymorphisms used in dataset simulation, in addition to correct SNPs we observe number of false positives.

Background

Many genome-wide association studies (GWAS) have been conducted in the search for specific genetic variants associated with common diseases. In testing for association with common polymorphisms, those variants that were identified were able to explain a modest proportion of disease heritability. This led to the hypothesis that multiple rare variants may play a role in complex disease etiology [1][2][3]. The multiple rare variants or common disease/rare variant hypothesis states that multiple rare variants with moderate to high penetrances underlie the susceptibility to a common disease. It is likely that both common and rare genetic variants contribute to disease risk.

Approaches targeted at uncovering associations between common polymorphisms and disease are generally underpowered for detecting the influence of rare variants. To identify disease-associated rare variants, investigators have proposed several methods based on the collapsing of low-frequency single-nucleotide polymorphisms (SNPs) [4–7].

In this analysis we use the methods proposed by Li and Leal [4] and Morris and Zeggini [8] to identify rare variants, and we use association analysis to identify common variants that confer liability to disease. The rationale behind this collapsing approach is that although the probability that an individual carries more than one rare allele may be low, in aggregate rare alleles may be common enough to account for variation in common traits. The goal is then to test for an association of an accumulation of rare minor alleles with the disease trait, by combining information across multiple variant sites.

We begin our analyses with the collapsing methods and extend the analyses in two ways. First, we use the mean log p-value (MLP) [9], which is a method that takes into account information about SNP function and ontologic pathway. The MLP can be thought of as a way to group together SNPs by their functional implication. It was originally developed for the analysis of gene expression data for a better understanding of the underlying mechanisms. Thus, by further analyzing the results of the rare variant collapsing approach using MLP analysis, we exploit both the spatial and the functional associations of SNPs implicated in a disease. Second, we use an empirical approach, stochastic gradient boosting (SGB), to discern groups of SNPs conferring liability to disease. SGB is an ensemble tree-based method [10] that is useful for empirically detecting sets of genes associated with a disease.

Methods

Data

Our analyses focus on the case-control data provided to the participants of Genetic Analysis Workshop 17 [11]. We selected the first of 200 simulated replicates for analysis. We conducted the Hardy-Weinberg equilibrium test using PLINK [12] and excluded from further analysis all markers with deviations (p-value less than 0.0001 in control subjects). We conducted population stratification analysis by first excluding correlated markers and then using the multidimensional scaling methods in PLINK.

Significant markers

We use the collapsing method proposed by Li and Leal [4] and Morris and Zeggini [8] to identify possible variants among rare SNPs. For this analysis we use the CCRaVAT (Case-Control Rare Variant Analysis Tool) software package [13]. The collapsing method is as follows: first, we divide the markers into groups on the basis of predefined criteria (either genes or sliding windows of defined sequence length); next, we collapse marker data based on an indicator variable that shows whether a subject carries any rare variants; and finally, using a Pearson chi-square test, we test the significance of the difference in proportions between case subjects and control subjects who carry rare variant minor alleles. When cell counts are small, we use a Fisher exact test instead.

We consider several approaches for the collapsing criterion, including gene-based collapsing and sliding windows of five different sizes (1 kb, 5 kb, 25 kb, 50 kb, and 100 kb). The resulting p-values are recorded for further analysis. In addition, we test the common variant SNPs using the Pearson chi-square test. Again, the resulting p-values are retained for further analysis.

MLP approach

We use the MLP approach [9] to incorporate functional and pathway information about genes into our analysis. The MLP analysis was developed in the context of gene expression analysis. The idea is to first assign a statistic (e.g., a p-value) to each gene. The genes are mapped onto gene sets or pathways by utilizing gene annotation databases, such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) [14], the Ingenuity Pathway Analysis (IPA) [15], and the Gene Ontology (GO) Biological Processes databases [16]. Permutation tests are used to determine a p-value for a gene set and to identify the top set of gene sets.

In our analysis, we explore both rare and common variants. To assign a p-value to a gene, we use the results of collapsing and association tests. Thus a gene can have multiple p-values associated with it, especially in the rare variant analysis, because the windows overlap. We examine several ways to assign the p-value to a gene in order to explore the utility of each, three for rare variants and two for common variants. For rare variants, we use (1) the p-value from the association test based on gene-wise collapsing or (2) the minimum p-value among association tests based on the collapsing within 5-kb sliding windows located in a gene. This window size is based on the preliminary explorations of varying window sizes ranging from 1 kb to 100 kb. In addition, we use (3) the mean log p-value among association tests based on the collapsing within 5-kb sliding windows located in a gene. For common variants, we use either (1) the minimum p-value among SNPs in a gene or (2) the mean log p-value among SNPs in a gene.

We create gene sets, consisting of groups of genes, on the basis of one of the databases (KEGG, IPA, or GO). The gene set statistic is subsequently calculated as the MLP of the gene statistics for each gene set. The permutation procedure described by Raghavan et al. [9] is used to obtain the gene set p-value. The gene sets are rank-ordered by p-value. We examine the top 20 sets and present the top 6 sets in this report.

Stochastic gradient boosting

SGB [10] is an ensemble tree-based method that uses an independently drawn random sample of individuals and SNPs to construct a small tree, typically containing 2 to 12 terminal nodes. The tree is grown as a result of recursively partitioning a node and contributes a small portion to the overall model. Each consecutive tree is built for the prediction residuals (from all preceding trees) of an independently drawn random sample. The final SGB model and its prediction perform by combining weighted individual tree contributions, with weight being a shrinkage parameter appropriately selected to reduce overfitting.

The SGB method produces a variable importance measure that can be used to identify top predictors. For tree methods variable importance scores show the relative contribution of each of the variables to predicting the outcome. For ensemble methods, such as SGB, the variable importance scores are averaged across all trees.

We apply the SGB method, using TreeNet, developed by Salford Systems [17], to the data consisting of the first replicate of the affected phenotype, multidimensional scaling components, environmental predictors (Age, Smoke, and Sex), and SNPs. We start with a set of the top common and rare SNPs passed from the collapsing approach and association tests. We then use the set of all SNPs provided in the GAW17 dataset.

Results

The initial data analysis of minor allele frequencies (MAFs) in the case-control data showed 3,224 rare variants (MAF between 1% and 5%) and 18,131 very rare variants (MAF less than 1%). There were 209 (30%) affected case subjects and 488 (70%) unaffected control subjects. Twenty-nine percent of males and thirty-one percent of females were affected. Smoking differed with case status. Smoking was prevalent in 35.9% of the affected subjects, in contrast to 21.7% of the unaffected subjects. After filtering, the data set included 22,615 SNPs. The population stratification results (the first two components) are shown in Figure 1. Three clusters were identified, corresponding to three populations. The resulting dimensions were carried forward for stratification.

Figure 1

Plot of the first two components of multidimensional scaling

We performed collapsing for various window sizes (1 kb, 5 kb, 25 kb, 50 kb, and 100 kb) as well as gene-wise collapsing. The Manhattan plot of p-values produced by the collapsing approach for 5-kb sliding windows is shown in Figure 2. The results from the 5-kb analysis were carried into further analyses.

We used the MLP approach to identify the top gene sets based on the KEGG, IPA, and GO databases. We examined the results from the MLP analysis using IPA pathways in greater detail, because this was the database used to simulate the GAW17 data. The top 20 gene sets were examined. Results for the top six sets are summarized in Table 1. The three gene statistics for rare variants and the two gene statistics for common variants described in the “MLP Approach” subsection of the Methods section are presented. The top gene sets based on the minimum of window p-values shows the Notch, Hypoxia, Nitric Oxide, and vascular endothelial growth factor (VEGF) signaling pathways. Both the VEGF and the Notch signaling pathways control initiation and differentiation in angiogenesis, a process leading to blood vessel formation or remodeling. The VEGF pathway is also among the top 15 pathways identified using the KEGG database and the same statistic (results not shown).

Table 1

Major IPA pathways identified by the MLP approach using five statistics for rare and common variants and their corresponding p-values

Statistic

Top 6 pathways selected

p-value

Rare variants

Gene-wise collapsing

1. Androgen and estrogen metabolism

0.0006

2. Sphingolipid metabolism

0.0006

3. Phenylalanine metabolism

0.0015

4. Death receptor signaling

0.002

5. Stilbene, coumarine, and lignin biosynthesis

0.0036

6. TWEAK signaling

0.0043

Minimum p-value of 5-kb sliding window collapsing within a gene

1. Notch signaling

0.0248

2. Hypoxia signaling in the cardiovascular system

0.0403

3. Nitric oxide signaling in the cardiovascular system

0.0409

4. VEGF signaling

0.0585

5. Glutamate receptor signaling

0.0642

6. Glutamate metabolism

0.0741

Mean log p-value of 5-kb sliding window collapsing within a gene

1. Cyanoamino acid metabolism

0.0032

2. Ubiquinone biosynthesis

0.0148

3. Nitrogen metabolism

0.0267

4. Alanine and aspartate metabolism

0.0392

5. GABA receptor signaling

0.0423

6. FXR/RXR activation

0.0438

Common variants

Minimum p-value among SNPs

1. Apoptosis signaling

0.0083

2. Pyrimidine metabolism

0.0232

3. CNTF signaling

0.0429

4. FLT3 signaling in hematopoietic progenitor cells

0.0601

5. Role of NANOG in mammalian embryonic stem cell pluripotency

0.0847

6. EGF signaling

0.0908

Mean log p-value for SNPs in a gene

1. Pyrimidine metabolism

0.0021

2. CNTF signaling

0.0127

3. Melanocyte development and pigmentation signaling

0.0221

4. JAK/Stat signaling

0.0327

5. IL-15 signaling

0.0356

6. FLT3 signaling in hematopoietic progenitor cells

0.045

We applied the SGB approach to the data containing preselected SNPs, population stratification results, and environmental variables. We used TreeNet to perform the SGB. We built 5,000 trees (iterations) with a maximum of 8 nodes. We chose a shrinkage parameter of 0.01 as appropriate for a data set of this dimension [18, 19]. The top set of SNPs was selected using a variable importance measure of 7.00 as a cutoff threshold. The corresponding genes were also recorded. The results of the SGB approach are shown in Table 2. The table contains the top SNPs selected and their corresponding genes. The results that match the simulated model are shown in bold. The SGB analyses using the complete set of SNPs did not show an improvement over prior runs (results not shown).

Table 2

Top SNPs and corresponding genes identified using the SGB approach

Top SNPs identified (from highest to lowest variable importance)

SNPs

C13S523

C9S3621

C6S6142

C5S237

C9S4860

C22S1374

C2S5630

C11S60

C13S905

C1S6542

C7S2893

C19S282

C6S1097

C10S2632

C2S955

C14S784

C1S5779

C7S2446

C9S1469

C12S4668

C2S1087

C2S2148

C1S10506

C6S2129

C15S3343

C12S4188

C14S1863

C2S4601

C6S2469

C10S2533

C19S609

C10S5515

C1S5530

C17S3017

C9S1225

C13S522

C12S3028

C5S3461

C19S1762

C1S9584

C19S1849

C9S4013

C22S1405

C12S622

C12S7056

C2S7558

C16S3421

C12S552

C2S4407

C1S996

C22S1351

C20S2310

C22S1158

C15S4060

C17S1262

C3S1305

C7S158

C10S387

C17S2377

C7S1877

C1S9718

C10S4422

C4S2872

C7S3971

C2S689

C8S3322

C10S6566

C14S20

C7S1076

C11S3224

C1S7413

C22S146

C8S4238

C8S4028

C18S2322

C6S6040

C12S5220

C6S6177

C19S3382

C19S2528

C1S9506

C4S4283

C12S3528

C11S2585

C17S2376

C12S5446

C17S4841

C1S10200

C4S2239

C7S3613

C5S4072

C11S6503

C11S4881

C1S10800

C9S123

C2S7414

C2S1139

C3S3962

C7S3490

C10S5783

C11S1683

C9S2613

C11S2532

C7S4111

C18S2310

C2S4079

C6S2366

C8S627

C2S6985

C1S7941

C11S5292

C4S4339

C3S3938

C6S5380

C22S875

C1S7092

C7S2590

C11S2871

C6S2216

C6S5677

C7S4646

C8S850

C8S271

C4S2296

C10S386

C9S5111

C15S3138

C1S7427

C17S3510

C3S96

C22S385

C1S3900

C3S4638

C21S672

C1S1388

C10S2683

C13S1168

C7S3697

C2S4909

C11S1280

C2S2154

C12S4591

C3S1176

C22S2039

C11S3320

C2S873

C9S3100

C2S7390

C12S5526

C11S1599

C6S4552

C1S10256

C10S3243

C12S5510

C4S2678

C4S2970

C2S8207

C16S560

C6S7138

C17S321

C20S1844

C12S5445

C10S2670

C1S4009

C17S2026

C9S3554

C13S1660

C14S590

C10S6432

C9S759

C19S4625

C1S9511

C8S934

C6S4242

C18S1560

C4S97

C15S3744

C7S397

C19S4658

C12S4534

C9S2083

C19S5271

C7S3898

C1S4838

C9S1607

C4S3834

C11S5644

C15S2848

C10S3777

C3S3657

C14S122

C14S3426

C16S1482

C4S3076

C9S2542

C2S6995

C21S778

C9S1835

C15S3559

C8S3416

C5S2032

C1S10813

C10S5690

C1S3676

C6S4400

C13S163

C22S645

C12S3039

C1S10164

C6S7164

C22S1222

C4S649

C19S277

C1S7408

C1S1542

C4S186

Genes

TNFRSF25

AHSA2

ADH1B

GPR85

OR13C5

SYTL2

PSME2

PTPRS

ARHGEF10L

RGPD3

C4ORF33

PTPRZ1

CDK5RAP2

PDGFD

VTI1B

OR10H3

KIF17

RGPD4

TKTL2

PAX4

STOM

EXPH5

BEGAIN

CYP4F2

PDE4B

ACVR1C

ANP32C

SMO

GOLGA1

OR8D4

PAQR5

ZNF486

PTGER3

LY75

PLEKHG4B

ZC3HC1

BRD3

CLEC2D

TSPAN3

NPHS1

MSH4

PPIG

ZNF474

AKR1B1

CARD9

KLRK1

ADAMTS7

ZNF576

STXBP3

WDR75

ABLIM3

SLC37A3

ECHDC3

OR6C1

ALPK3

LYPD5TMC4

HIPK1

LOC729332

FOXI1

GATA4

ERCC6

OR6C65

SLC28A1

C20ORF32

VTCN1

UGT1A10

PGBD1

TNFRSF10D

PGBD3

SRGAP1

AKAP13

PRIC285

ARNT

COL6A3

BAT2

CDCA2

HKDC1

PLXNC1

ZNF213

PIGP

ADAM15

MTERFD2

SLC44A4

PTK2B

MAT1A

C12ORF63

USP31

ETS2

OR10J1

CRELD1

PSMB8

EXT1

CYP2C8

CHPT1

LOC100132786

HIRA

OR10J5

LOC100130135

MDN1

SAMD12

SORCS1

TRAFD1

TRPV3

ARVCF

UAP1

SETD2

VNN1

MLZE

CASP7

CAMKK2

KCNJ12

TOP3B

LRRN2

GOLGB1

FUCA2

TG

DCLRE1A

ANAPC5

SLFN13

SUSD2

NUAK2

TMCC1

UTRN

ANKRD15

TACC2

ZNF26

PIP4K2B

SEC14L3

IKBKE

TRPC1

AGPAT4

PDCD1LG2

ATHL1

TNFRSF19

BRCA1

SMTN

LAMB3

CRIPAK

PMS2

AQP3

GALNTL4

FLT1

FAM117A

LIMK2

DUSP10

GRK4

TSPAN13

C9ORF131

HPS5

TRPC4

RECQL5

LOC100132621

KIAA0133

AFAP1

NPC1L1

NPR2

NELL1

FREM2

HRH4

ARL6IP2

PF4V1

NCF1

POLR1E

OR8H1

PIBF1

B4GALT6

OXER1

STBD1

FBXO24

SMC5

OR9G4

RNASE6

MCART2

FSHR

FAM13A1

FBXL13

C9ORF79

OR4D9

FLJ10357

ZNF57

BCL11A

PDLIM5

RELN

ROR2

AHNAK

ACIN1

ZNF77

Boldface indicates results that match the simulated model.

The MLP approach correctly placed the VEGF signaling pathway and the two pathways related to the cardiovascular system (Hypoxia and Nitric Oxide) among the top six pathways. There are only 5 (out of 216) SNPs correctly identified using the SGB approach; their MAFs range from 0.2% to 17%. Most of the SNPs (211) placed in the top list are false positives.

Discussion and conclusions

Traditionally, GWAS test for association of disease with common polymorphisms. Polymorphisms with population frequencies of 5% or more could be tested directly or indirectly for association with disease risk or quantitative traits, and GWAS have identified many genetic variants associated with disease traits. Replication of these results has not been consistently successful. More recently, methods have been proposed to identify multiple rare variants with small individual effect sizes that may be implicated in complex multigenic diseases. These methods are based on grouping rare variants by their physical proximity in order to combine information across them. According to Hirschhorn [20], one of the primary goals of GWAS is to discover the biologic pathways underlying polygenic diseases.

The MLP method is an effort in this direction, where both common and rare variants are considered on the basis of their functional implication in disease etiology. Our goal here was to exploit both the spatial and the functional associations of SNPs implicated in a disease to identify the underlying biologic pathways. Pathways used to simulate the affection status in the GAW17 data set were among the top four pathways identified by the MLP approach based on statistic 2 for rare variants. More specifically, the three pathways (Hypoxia, Nitric Oxide, and VEGF) were used to simulate the data. In addition, as explained in what follows, the top four pathways, including the Notch pathway, may be part of a cascade of interrelated pathways. Enriched signaling pathways in our analysis may overlap functionally and indicate processes leading to angiogenesis. Hypoxia signaling can trigger the VEGF cascade in cancer tissue angiogenesis, and the Notch processes downstream from the Hypoxia and VEGF pathways lead to a differentiation of newly formed vessels. Notch signaling can also down-regulate VEGF expression in a feedback loop. Thus the MLP approach based on statistic 2 for rare variants appears to be able to also identify related pathways and may be promising for the discovery of biological pathways implicated in disease etiology by rare variants.

One of the goals of this analysis was to compare results from a variety of functional and pathway databases and from a number of gene statistics. Our results indicate that using the IPA database, the gene statistic that identified the most relevant pathways was the minimum p-value derived from collapsing rare variants within 5-kb sliding windows residing in a gene. These results highlight the importance of using the most appropriate pathway database for the analysis, an aspect not explicitly discussed in the literature. Our analysis also indicates that (1) the results can be influenced by the density of coverage of rare variant SNPs in a gene; (2) gene-based collapsing may be too broad and may dilute the underlying information; and (3) using the mean of the window p-values may mute the signal considerably. In future work, we will evaluate alternative approaches to mapping SNP-level p-values to gene-level p-values as well as methods for combining the rare and common variants analyses.

The top pathways identified using the MLP method intersect with the pathways that contain genes from the results of the SGB approach. There are 5 correct SNPs out of 216 residing in 5 correct genes out of 188 corresponding to the top selected SNPs. The large number of false positives may be due to correlation of the SNPs identified by the SGB approach with the true SNPs used in the simulation model. Our current work includes studying methods to bridge the two approaches utilizing the functional information and the statistical correlation, respectively.

Currently, many analytic strategies rely on GWAS with one-SNP-at-a-time analyses. Although this approach has certainly yielded many promising candidates, it requires large samples to mitigate type I errors. One-SNP-at-a-time analyses generally do not take advantage of all the information present in the data, and failure of replication is commonplace. Recent attempts have been made to incorporate information from rare variants into the analysis by aggregating across SNPs that are in close proximity to each other. We have extended this further by leveraging information from SNPs that are either functionally related (MLP approach) or statistically correlated (SGB approach) with the hope of obtaining results that are credible and logically interpretable. These methods would, of course, need to be further evaluated in other data sets and other settings.

Declarations

Acknowledgments

The Genetic Analysis Workshop is supported by National Institutes of Health grant R01 GM031575. We are grateful to Jesse Berlin and Paul Stang for their help with access to the data and resources. We would also like to thank Qingqin Li and Dai Wang for several discussions about the design and methods we used.

Competing interests

The authors declare that there are no competing interests.

Authors’ contributions

YC carried out the collapsing approach, pathway and tree-based analyses and participated in writing the manuscript. MW, NR, SF designed the analysis plan and participated in writhing the manuscript. FD prepared the data for analysis and edited the manuscript. All authors read and approved the final manuscript.

Copyright

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.