Analyzing genome-wide association studies with an FDR controlling modification of the Bayesian Information Criterion

A neutral comparison of algorithms to minimize L0 penalties for high-dimensional variable selection

Variable selection methods based on L0 penalties have excellent theoretical properties to select sparse models in a high-dimensional setting. There exist modifications of the Bayesian Information Criterion (BIC) which either control the familywise error rate (mBIC) or the false discovery rate (mBIC2) in terms of which regressors are selected to enter a model. However, the minimization of L0 penalties comprises a mixed-integer problem which is known to be NP-hard and therefore becomes computationally challenging with increasing numbers of regressor variables. This is one reason why alternatives like the LASSO have become so popular, which involve convex optimization problems that are easier to solve. The last few years have seen some real progress in developing new algorithms to minimize L0 penalties. The aim of this article is to compare the performance of these algorithms in terms of minimizing L0 -based selection criteria. Simulation studies covering a wide range of scenarios that are inspired by genetic association studies are used to compare the values of selection criteria obtained with different algorithms. In addition, some statistical characteristics of the selected models and the runtime of algorithms are compared. Finally, the performance of the algorithms is illustrated in a real data example concerned with expression quantitative trait loci (eQTL) mapping.

Identification of Driver Epistatic Gene Pairs Combining Germline and Somatic Mutations in Cancer

Cancer arises from the complex interplay of various factors. Traditionally, the identification of driver genes focuses primarily on the analysis of somatic mutations. We describe a new method for the detection of driver gene pairs based on an epistasis analysis that considers both germline and somatic variations. Specifically, the identification of significantly mutated gene pairs entails the calculation of a contingency table, wherein one of the co-mutated genes can exhibit a germline variant. By adopting this approach, it is possible to select gene pairs in which the individual genes do not exhibit significant associations with cancer. Finally, a survival analysis is used to select clinically relevant gene pairs. To test the efficacy of the new algorithm, we analyzed the colon adenocarcinoma (COAD) and lung adenocarcinoma (LUAD) samples available at The Cancer Genome Atlas (TCGA). In the analysis of the COAD and LUAD samples, we identify epistatic gene pairs significantly mutated in tumor tissue with respect to normal tissue. We believe that further analysis of the gene pairs detected by our method will unveil new biological insights, enhancing a better description of the cancer mechanism.

BICOSS: Bayesian iterative conditional stochastic search for GWAS

Background

Single marker analysis (SMA) with linear mixed models for genome wide association studies has uncovered the contribution of genetic variants to many observed phenotypes. However, SMA has weak false discovery control. In addition, when a few variants have large effect sizes, SMA has low statistical power to detect small and medium effect sizes, leading to low recall of true causal single nucleotide polymorphisms (SNPs).

Results

We present the Bayesian Iterative Conditional Stochastic Search (BICOSS) method that controls false discovery rate and increases recall of variants with small and medium effect sizes. BICOSS iterates between a screening step and a Bayesian model selection step. A simulation study shows that, when compared to SMA, BICOSS dramatically reduces false discovery rate and allows for smaller effect sizes to be discovered. Finally, two real world applications show the utility and flexibility of BICOSS.

Conclusions

When compared to widely used SMA, BICOSS provides higher recall of true SNPs while dramatically reducing false discovery rate.

A genome-wide association study of plasma concentrations of warfarin enantiomers and metabolites in sub-Saharan black-African patients

Diversity in pharmacogenomic studies is poor, especially in relation to the inclusion of black African patients. Lack of funding and difficulties in recruitment, together with the requirement for large sample sizes because of the extensive genetic diversity in Africa, are amongst the factors which have hampered pharmacogenomic studies in Africa. Warfarin is widely used in sub-Saharan Africa, but as in other populations, dosing is highly variable due to genetic and non-genetic factors. In order to identify genetic factors determining warfarin response variability, we have conducted a genome-wide association study (GWAS) of plasma concentrations of warfarin enantiomers/metabolites in sub-Saharan black-Africans. This overcomes the issue of non-adherence and may have greater sensitivity at genome-wide level, to identify pharmacokinetic gene variants than focusing on mean weekly dose, the usual end-point used in previous studies. Participants recruited at 12 outpatient sites in Uganda and South Africa on stable warfarin dose were genotyped using the Illumina Infinium H3Africa Consortium Array v2. Imputation was conducted using the 1,000 Genomes Project phase III reference panel. Warfarin/metabolite plasma concentrations were determined by high-performance liquid chromatography with tandem mass spectrometry. Multivariable linear regression was undertaken, with adjustment made for five non-genetic covariates and ten principal components of genetic ancestry. After quality control procedures, 548 participants and 17,268,054 SNPs were retained. CYP2C9*8, CYP2C9*9, CYP2C9*11, and the CYP2C cluster SNP rs12777823 passed the Bonferroni-adjusted replication significance threshold (p < 3.21E-04) for warfarin/metabolite ratios. In an exploratory GWAS analysis, 373 unique SNPs in 13 genes, including CYP2C9*8, passed the Bonferroni-adjusted genome-wide significance threshold (p < 3.846E-9), with 325 (87%, all located on chromosome 10) SNPs being associated with the S-warfarin/R-warfarin outcome (top SNP rs11188082, CYP2C19 intron variant, p = 1.55E-17). Approximately 69% of these SNPs were in linkage disequilibrium (r² > 0.8) with CYP2C9*8 (n = 216) and rs12777823 (n = 8). Using a pharmacokinetic approach, we have shown that variants other than CYP2C9*2 and CYP2C9*3 are more important in sub-Saharan black-Africans, mainly due to the allele frequencies. In exploratory work, we conducted the first warfarin pharmacokinetics-related GWAS in sub-Saharan Africans and identified novel SNPs that will require external replication and functional characterization before they can be considered for inclusion in warfarin dosing algorithms.

Selecting predictive biomarkers from genomic data

Recently there have been tremendous efforts to develop statistical procedures which allow to determine subgroups of patients for which certain treatments are effective. This article focuses on the selection of prognostic and predictive genetic biomarkers based on a relatively large number of candidate Single Nucleotide Polymorphisms (SNPs). We consider models which include prognostic markers as main effects and predictive markers as interaction effects with treatment. We compare different high-dimensional selection approaches including adaptive lasso, a Bayesian adaptive version of the Sorted L-One Penalized Estimator (SLOBE) and a modified version of the Bayesian Information Criterion (mBIC2). These are compared with classical multiple testing procedures for individual markers. Having identified predictive markers we consider several different approaches how to specify subgroups susceptible to treatment. Our main conclusion is that selection based on mBIC2 and SLOBE has similar predictive performance as the adaptive lasso while including substantially fewer biomarkers.

High-LD SNP markers exhibiting pleiotropic effects on salt tolerance at germination and seedlings stages in spring wheat

Salt tolerance at germination and seedling growth stages was investigated. GWAS revealed nine genomic regions with pleiotropic effects on salt tolerance. Salt tolerant genotypes were identified for future breeding program. With 20% of the irrigated land worldwide affected by it, salinity is a serious threat to plant development and crop production. While wheat is the most stable food source worldwide, it has been classified as moderately tolerant to salinity. In several crop plants; such as barley, maize and rice, it has been shown that salinity tolerance at seed germination and seedling establishment is under polygenic control. As yield was the ultimate goal of breeders and geneticists, less attention has been paid to understanding the genetic architecture of salt tolerance at early stages. Thus, the genetic control of salt tolerance at these stages is poorly understood relative to the late stages. In the current study, 176 genotypes of spring wheat were tested for salinity tolerance at seed germination and seedling establishment. Genome-Wide Association Study (GWAS) has been used to identify the genomic regions/genes conferring salt tolerance at seed germination and seedling establishment. Salinity stress negatively impacted all germination and seedling development parameters. A set of 137 SNPs showed significant association with the traits of interest. Across the whole genome, 33 regions showed high linkage disequilibrium (LD). These high LD regions harbored 15 SNPs with pleiotropic effect (i.e. SNPs that control more than one trait). Nine genes belonging to different functional groups were found to be associated with the pleiotropic SNPs. Noteworthy, chromosome 2B harbored the gene TraesCS2B02G135900 that acts as a potassium transporter. Remarkably, one SNP marker, reported in an early study, associated with salt tolerance was validated in this study. Our findings represent potential targets of genetic manipulation to understand and improve salinity tolerance in wheat.

What has GWAS done for HLA and disease associations?

The major histocompatibility complex (MHC) is located in chromosome 6p21 and contains crucial regulators of immune response, including human leucocyte antigen (HLA) genes, alongside other genes with nonimmunological roles. More recently, a repertoire of noncoding RNA genes, including expressed pseudogenes, has also been identified. The MHC is the most gene dense and most polymorphic part of the human genome. The region exhibits haplotype-specific linkage disequilibrium patterns, contains the strongest cis- and trans-eQTLs/meQTLs in the genome and is known as a hot spot for disease associations. Another layer of complexity is provided to the region by the extreme structural variation and copy number variations. While the HLA-B gene has the highest number of alleles, the HLA-DR/DQ subregion is structurally most variable and shows the highest number of disease associations. Reliance on a single reference sequence has complicated the design, execution and analysis of GWAS for the MHC region and not infrequently, the MHC region has even been excluded from the analysis of GWAS data. Here, we contrast features of the MHC region with the rest of the genome and highlight its complexities, including its functional polymorphisms beyond those determined by single nucleotide polymorphisms or single amino acid residues. One of the several issues with customary GWAS analysis is that it does not address this additional layer of polymorphisms unique to the MHC region. We highlight alternative approaches that may assist with the analysis of GWAS data from the MHC region and unravel associations with all functional polymorphisms beyond single SNPs. We suggest that despite already showing the highest number of disease associations, the true extent of the involvement of the MHC region in disease genetics may not have been uncovered.

Bayesian and frequentist analysis of an Austrian genome-wide association study of colorectal cancer and advanced adenomas

Most genome-wide association studies (GWAS) were analyzed using single marker tests in combination with stringent correction procedures for multiple testing. Thus, a substantial proportion of associated single nucleotide polymorphisms (SNPs) remained undetected and may account for missing heritability in complex traits. Model selection procedures present a powerful alternative to identify associated SNPs in high-dimensional settings. In this GWAS including 1060 colorectal cancer cases, 689 cases of advanced colorectal adenomas and 4367 controls we pursued a dual approach to investigate genome-wide associations with disease risk applying both, single marker analysis and model selection based on the modified Bayesian information criterion, mBIC2, implemented in the software package MOSGWA. For different case-control comparisons, we report models including between 1-14 candidate SNPs. A genome-wide significant association of rs17659990 (P=5.43×10^-9, DOCK3, chromosome 3p21.2) with colorectal cancer risk was observed. Furthermore, 56 SNPs known to influence susceptibility to colorectal cancer and advanced adenoma were tested in a hypothesis-driven approach and several of them were found to be relevant in our Austrian cohort. After correction for multiple testing (α=8.9×10^-4), the most significant associations were observed for SNPs rs10505477 (P=6.08×10^-4) and rs6983267 (P=7.35×10^-4) of CASC8, rs3802842 (P=8.98×10^-5, COLCA1,2), and rs12953717 (P=4.64×10^-4, SMAD7). All previously unreported SNPs demand replication in additional samples. Reanalysis of existing GWAS datasets using model selection as tool to detect SNPs associated with a complex trait may present a promising resource to identify further genetic risk variants not only for colorectal cancer.

Joint genotype- and ancestry-based genome-wide association studies in admixed populations

In genome-wide association studies (GWAS) genetic loci that influence complex traits are localized by inspecting associations between genotypes of genetic markers and the values of the trait of interest. On the other hand, admixture mapping, which is performed in case of populations consisting of a recent mix of two ancestral groups, relies on the ancestry information at each locus (locus-specific ancestry). Recently it has been proposed to jointly model genotype and locus-specific ancestry within the framework of single marker tests. Here, we extend this approach for population-based GWAS in the direction of multimarker models. A modified version of the Bayesian information criterion is developed for building a multilocus model that accounts for the differential correlation structure due to linkage disequilibrium (LD) and admixture LD. Simulation studies and a real data example illustrate the advantages of this new approach compared to single-marker analysis or modern model selection strategies based on separately analyzing genotype and ancestry data, as well as to single-marker analysis combining genotypic and ancestry information. Depending on the signal strength, our procedure automatically chooses whether genotypic or locus-specific ancestry markers are added to the model. This results in a good compromise between the power to detect causal mutations and the precision of their localization. The proposed method has been implemented in R and is available at http://www.math.uni.wroc.pl/~mbogdan/admixtures/.

Controlling the Rate of GWAS False Discoveries

With the rise of both the number and the complexity of traits of interest, control of the false discovery rate (FDR) in genetic association studies has become an increasingly appealing and accepted target for multiple comparison adjustment. While a number of robust FDR-controlling strategies exist, the nature of this error rate is intimately tied to the precise way in which discoveries are counted, and the performance of FDR-controlling procedures is satisfactory only if there is a one-to-one correspondence between what scientists describe as unique discoveries and the number of rejected hypotheses. The presence of linkage disequilibrium between markers in genome-wide association studies (GWAS) often leads researchers to consider the signal associated to multiple neighboring SNPs as indicating the existence of a single genomic locus with possible influence on the phenotype. This a posteriori aggregation of rejected hypotheses results in inflation of the relevant FDR. We propose a novel approach to FDR control that is based on prescreening to identify the level of resolution of distinct hypotheses. We show how FDR-controlling strategies can be adapted to account for this initial selection both with theoretical results and simulations that mimic the dependence structure to be expected in GWAS. We demonstrate that our approach is versatile and useful when the data are analyzed using both tests based on single markers and multiple regression. We provide an R package that allows practitioners to apply our procedure on standard GWAS format data, and illustrate its performance on lipid traits in the North Finland Birth Cohort 66 cohort study.

An Adaptive Ridge Procedure for L0 Regularization

Penalized selection criteria like AIC or BIC are among the most popular methods for variable selection. Their theoretical properties have been studied intensively and are well understood, but making use of them in case of high-dimensional data is difficult due to the non-convex optimization problem induced by L0 penalties. In this paper we introduce an adaptive ridge procedure (AR), where iteratively weighted ridge problems are solved whose weights are updated in such a way that the procedure converges towards selection with L0 penalties. After introducing AR its specific shrinkage properties are studied in the particular case of orthogonal linear regression. Based on extensive simulations for the non-orthogonal case as well as for Poisson regression the performance of AR is studied and compared with SCAD and adaptive LASSO. Furthermore an efficient implementation of AR in the context of least-squares segmentation is presented. The paper ends with an illustrative example of applying AR to analyze GWAS data.

ICN: a normalization method for gene expression data considering the over-expression of informative genes

The global increase of gene expression has been frequently established in cancer microarray studies.

Further improvements to linear mixed models for genome-wide association studies

We examine improvements to the linear mixed model (LMM) that better correct for population structure and family relatedness in genome-wide association studies (GWAS). LMMs rely on the estimation of a genetic similarity matrix (GSM), which encodes the pairwise similarity between every two individuals in a cohort. These similarities are estimated from single nucleotide polymorphisms (SNPs) or other genetic variants. Traditionally, all available SNPs are used to estimate the GSM. In empirical studies across a wide range of synthetic and real data, we find that modifications to this approach improve GWAS performance as measured by type I error control and power. Specifically, when only population structure is present, a GSM constructed from SNPs that well predict the phenotype in combination with principal components as covariates controls type I error and yields more power than the traditional LMM. In any setting, with or without population structure or family relatedness, a GSM consisting of a mixture of two component GSMs, one constructed from all SNPs and another constructed from SNPs that well predict the phenotype again controls type I error and yields more power than the traditional LMM. Software implementing these improvements and the experimental comparisons are available at http://microsoft.com/science.