Germline HLA landscape does not predict efficacy of pembrolizumab monotherapy across solid tumor types.

We evaluated the impact of class I and class II human leukocyte antigen (HLA) genotypes, heterozygosity, and diversity on the efficacy of pembrolizumab. Seventeen pembrolizumab clinical trials across eight tumor types and one basket trial in patients with advanced solid tumors were included (n > 3,500 analyzed). Germline DNA was genotyped using a custom genotyping array. HLA diversity (measured by heterozygosity and evolutionary divergence) across class I loci was not associated with improved response to pembrolizumab, either within each tumor type evaluated or across all patients. Similarly, HLA heterozygosity at each class I and class II gene was not associated with response to pembrolizumab after accounting for the number of tests conducted. No conclusive association between HLA genotype and response to pembrolizumab was identified in this dataset. Germline HLA genotype or diversity alone is not an important independent determinant of response to pembrolizumab and should not be used for clinical decision-making in patients treated with pembrolizumab.

GWAS of longitudinal trajectories at biobank scale.

Biobanks linked to massive, longitudinal electronic health record (EHR) data make numerous new genetic research questions feasible. One among these is the study of biomarker trajectories. For example, high blood pressure measurements over visits strongly predict stroke onset, and consistently high fasting glucose and Hb1Ac levels define diabetes. Recent research reveals that not only the mean level of biomarker trajectories but also their fluctuations, or within-subject (WS) variability, are risk factors for many diseases. Glycemic variation, for instance, is recently considered an important clinical metric in diabetes management. It is crucial to identify the genetic factors that shift the mean or alter the WS variability of a biomarker trajectory. Compared to traditional cross-sectional studies, trajectory analysis utilizes more data points and captures a complete picture of the impact of time-varying factors, including medication history and lifestyle. Currently, there are no efficient tools for genome-wide association studies (GWASs) of biomarker trajectories at the biobank scale, even for just mean effects. We propose TrajGWAS, a linear mixed effect model-based method for testing genetic effects that shift the mean or alter the WS variability of a biomarker trajectory. It is scalable to biobank data with 100,000 to 1,000,000 individuals and many longitudinal measurements and robust to distributional assumptions. Simulation studies corroborate that TrajGWAS controls the type I error rate and is powerful. Analysis of eleven biomarkers measured longitudinally and extracted from UK Biobank primary care data for more than 150,000 participants with 1,800,000 observations reveals loci that significantly alter the mean or WS variability.

Multi-trait analysis of gene-by-environment interactions in large-scale genetic studies.

Identifying genotype-by-environment interaction (GEI) is challenging because the GEI analysis generally has low power. Large-scale consortium-based studies are ultimately needed to achieve adequate power for identifying GEI. We introduce Multi-Trait Analysis of Gene-Environment Interactions (MTAGEI), a powerful, robust, and computationally efficient framework to test gene-environment interactions on multiple traits in large data sets, such as the UK Biobank (UKB). To facilitate the meta-analysis of GEI studies in a consortium, MTAGEI efficiently generates summary statistics of genetic associations for multiple traits under different environmental conditions and integrates the summary statistics for GEI analysis. MTAGEI enhances the power of GEI analysis by aggregating GEI signals across multiple traits and variants that would otherwise be difficult to detect individually. MTAGEI achieves robustness by combining complementary tests under a wide spectrum of genetic architectures. We demonstrate the advantages of MTAGEI over existing single-trait-based GEI tests through extensive simulation studies and the analysis of the whole exome sequencing data from the UKB.

Applying polygenic risk score methods to pharmacogenomics GWAS: challenges and opportunities.

Polygenic risk scores (PRSs) have emerged as promising tools for the prediction of human diseases and complex traits in disease genome-wide association studies (GWAS). Applying PRSs to pharmacogenomics (PGx) studies has begun to show great potential for improving patient stratification and drug response prediction. However, there are unique challenges that arise when applying PRSs to PGx GWAS beyond those typically encountered in disease GWAS (e.g. Eurocentric or trans-ethnic bias). These challenges include: (i) the lack of knowledge about whether PGx or disease GWAS/variants should be used in the base cohort (BC); (ii) the small sample sizes in PGx GWAS with corresponding low power and (iii) the more complex PRS statistical modeling required for handling both prognostic and predictive effects simultaneously. To gain insights in this landscape about the general trends, challenges and possible solutions, we first conduct a systematic review of both PRS applications and PRS method development in PGx GWAS. To further address the challenges, we propose (i) a novel PRS application strategy by leveraging both PGx and disease GWAS summary statistics in the BC for PRS construction and (ii) a new Bayesian method (PRS-PGx-Bayesx) to reduce Eurocentric or cross-population PRS prediction bias. Extensive simulations are conducted to demonstrate their advantages over existing PRS methods applied in PGx GWAS. Our systematic review and methodology research work not only highlights current gaps and key considerations while applying PRS methods to PGx GWAS, but also provides possible solutions for better PGx PRS applications and future research.

A fast and powerful linear mixed model approach for genotype-environment interaction tests in large-scale GWAS.

Genotype-by-environment interaction (GEI or GxE) plays an important role in understanding complex human traits. However, it is usually challenging to detect GEI signals efficiently and accurately while adjusting for population stratification and sample relatedness in large-scale genome-wide association studies (GWAS). Here we propose a fast and powerful linear mixed model-based approach, fastGWA-GE, to test for GEI effect and G + GxE joint effect. Our extensive simulations show that fastGWA-GE outperforms other existing GEI test methods by controlling genomic inflation better, providing larger power and running hundreds to thousands of times faster. We performed a fastGWA-GE analysis of ~7.27 million variants on 452 249 individuals of European ancestry for 13 quantitative traits and five environment variables in the UK Biobank GWAS data and identified 96 significant signals (72 variants across 57 loci) with GEI test P-values < 1 × 10-9, including 27 novel GEI associations, which highlights the effectiveness of fastGWA-GE in GEI signal discovery in large-scale GWAS.

Integrating multiple traits for improving polygenic risk prediction in disease and pharmacogenomics GWAS.

Polygenic risk score (PRS) has been recently developed for predicting complex traits and drug responses. It remains unknown whether multi-trait PRS (mtPRS) methods, by integrating information from multiple genetically correlated traits, can improve prediction accuracy and power for PRS analysis compared with single-trait PRS (stPRS) methods. In this paper, we first review commonly used mtPRS methods and find that they do not directly model the underlying genetic correlations among traits, which has been shown to be useful in guiding multi-trait association analysis in the literature. To overcome this limitation, we propose a mtPRS-PCA method to combine PRSs from multiple traits with weights obtained from performing principal component analysis (PCA) on the genetic correlation matrix. To accommodate various genetic architectures covering different effect directions, signal sparseness and across-trait correlation structures, we further propose an omnibus mtPRS method (mtPRS-O) by combining P values from mtPRS-PCA, mtPRS-ML (mtPRS based on machine learning) and stPRSs using Cauchy Combination Test. Our extensive simulation studies show that mtPRS-PCA outperforms other mtPRS methods in both disease and pharmacogenomics (PGx) genome-wide association studies (GWAS) contexts when traits are similarly correlated, with dense signal effects and in similar effect directions, and mtPRS-O is consistently superior to most other methods due to its robustness under various genetic architectures. We further apply mtPRS-PCA, mtPRS-O and other methods to PGx GWAS data from a randomized clinical trial in the cardiovascular domain and demonstrate performance improvement of mtPRS-PCA in both prediction accuracy and patient stratification as well as the robustness of mtPRS-O in PRS association test.

Robust genetic model-based SNP-set association test using CauchyGM.

MOTIVATION: Association testing on genome-wide association studies (GWAS) data is commonly performed under a single (mostly additive) genetic model framework. However, the underlying true genetic mechanisms are often unknown in practice for most complex traits. When the employed inheritance model deviates from the underlying model, statistical power may be reduced. To overcome this challenge, an integrative association test that directly infers the underlying genetic model from GWAS data has previously been proposed for single-SNP analysis. RESULTS: In this article, we propose a Cauchy combination Genetic Model-based association test (CauchyGM) under a generalized linear model framework for SNP-set level analysis. CauchyGM does not require prior knowledge on the underlying inheritance pattern of each SNP. It performs a score test that first estimates an individual P-value of each SNP in an SNP-set with both minor allele frequency (MAF) > 1% and three genotypes and further aggregates the rest SNPs using SKAT. CauchyGM then combines the correlated P-values across multiple SNPs and different genetic models within the set using Cauchy Combination Test. To further accommodate both sparse and dense signal patterns, we also propose an omnibus association test (CauchyGM-O) by combining CauchyGM with SKAT and the burden test. Our extensive simulations show that both CauchyGM and CauchyGM-O maintain the type I error well at the genome-wide significance level and provide substantial power improvement compared to existing methods. We apply our methods to a pharmacogenomic GWAS data from a large cardiovascular randomized clinical trial. Both CauchyGM and CauchyGM-O identify several novel genome-wide significant genes. AVAILABILITY AND IMPLEMENTATION: The R package CauchyGM is publicly available on github: https://github.com/ykim03517/CauchyGM. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

AWOT and CWOT for genotype and genotype-by-treatment interaction joint analysis in pharmacogenetics GWAS.

MOTIVATION: Pharmacogenomics (PGx) research holds the promise for detecting association between genetic variants and drug responses in randomized clinical trials, but it is limited by small populations and thus has low power to detect signals. It is critical to increase the power of PGx genome-wide association studies (GWAS) with small sample sizes so that variant-drug-response association discoveries are not limited to common variants with extremely large effect. RESULTS: In this article, we first discuss the challenges of PGx GWAS studies and then propose the adaptively weighted joint test (AWOT) and Cauchy Weighted jOint Test (CWOT), which are two flexible and robust joint tests of the single nucleotide polymorphism main effect and genotype-by-treatment interaction effect for continuous and binary endpoints. Two analytic procedures are proposed to accurately calculate the joint test P-value. We evaluate AWOT and CWOT through extensive simulations under various scenarios. The results show that the proposed AWOT and CWOT control type I error well and outperform existing methods in detecting the most interesting signal patterns in PGx settings (i.e. with strong genotype-by-treatment interaction effects, but weak genotype main effects). We demonstrate the value of AWOT and CWOT by applying them to the PGx GWAS from the Bezlotoxumab Clostridium difficile MODIFY I/II Phase 3 trials. AVAILABILITY AND IMPLEMENTATION: The R package COWT is publicly available on CRAN https://cran.r-project.org/web/packages/cwot/index.html. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Composite Kernel Association Test (CKAT) for SNP-set joint assessment of genotype and genotype-by-treatment interaction in Pharmacogenetics studies.

MOTIVATION: It is of substantial interest to discover novel genetic markers that influence drug response in order to develop personalized treatment strategies that maximize therapeutic efficacy and safety. To help enable such discoveries, we focus on testing the association between the cumulative effect of multiple single nucleotide polymorphisms (SNPs) in a particular genomic region and a drug response of interest. However, the currently existing methods are either computational inefficient or not able to control type I error and provide decent power for whole exome or genome analysis in Pharmacogenetics (PGx) studies with small sample sizes. RESULTS: In this article, we propose the Composite Kernel Association Test (CKAT), a flexible and robust kernel machine-based approach to jointly test the genetic main effect and SNP-treatment interaction effect for SNP-sets in Pharmacogenetics (PGx) assessments embedded within randomized clinical trials. An analytic procedure is developed to accurately calculate the P-value so that computationally extensive procedures (e.g. permutation or perturbation) can be avoided. We evaluate CKAT through extensive simulation studies and application to the gene-level association test of the reduction in Clostridium difficile infection recurrence in patients treated with bezlotoxumab. The results demonstrate that the proposed CKAT controls type I error well for PGx studies, is efficient for whole exome/genome association analysis and provides better power performance than existing methods across multiple scenarios. AVAILABILITY AND IMPLEMENTATION: The R package CKAT is publicly available on CRAN https://cran.r-project.org/web/packages/CKAT/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Effect of Endogenous Clostridioides difficile Toxin Antibodies on Recurrence of C. difficile Infection.

BACKGROUND: Endogenous antibodies (eAbs) against Clostridioides (Clostridium) difficile toxins may protect against recurrence of C. difficile infection (rCDI). This hypothesis was tested using placebo group data from MODIFY (Monoclonal Antibodies for C. difficile Therapy) I and II (NCT01241552 and NCT01513239, respectively), global, randomized phase 3 trials that assessed the efficacy and safety of the antitoxin monoclonal antibodies bezlotoxumab and actoxumab in participants receiving antibiotic therapy for CDI. METHODS: A placebo infusion (normal saline) was administered on study day 1. Serum samples were collected on day 1, week 4, and week 12, and eAb-A and eAb-B titers were measured by 2 validated electrochemiluminescence immunoassays. Rates of initial clinical cure and rCDI were summarized by eAb titer category (low, medium, high) at each time point. RESULTS: Serum eAb titers were available from a total of 773 participants. The proportion of participants with high eAb-A and eAb-B titers increased over time. Rates of initial clinical cure were similar across eAb titer categories. There was no correlation between eAb-A titers and rCDI rate at any time point. However, there was a negative correlation between rCDI and eAb-B titer on day 1 and week 4. rCDI occurred in 22% of participants with high eAb-B titers at baseline compared with 35% with low or medium titers (P = .015). CONCLUSIONS: Higher eAb titers against toxin B, but not toxin A, were associated with protection against rCDI. These data are consistent with the observed efficacy of bezlotoxumab, and lack of efficacy of actoxumab, in the MODIFY trials. CLINICAL TRIALS REGISTRATION: NCT01241552 and NCT01513239.

Bezlotoxumab for prevention of Clostridium difficile infection recurrence: Distinguishing relapse from reinfection with whole genome sequencing.

BACKGROUND: Bezlotoxumab has been shown to prevent Clostridium difficile infection recurrence (rCDI) in high-risk patients. METHODS: We used whole genome sequencing to estimate the impact of bezlotoxumab on same-strain relapse or new-strain reinfection in MODIFY I/II trials. Reinfection with a new strain and relapse with the same strain were differentiated by the comparison of ribotype (RT) and pair-wise single-nucleotide whole genome sequencing (WGS) variations (PWSNV). Relapse was assigned if the baseline RT and the RT isolated during rCDI were the same, and if PWSNVs were ≤ 2. Reinfection was assigned if the baseline RT and the RT isolated during rCDI were different, or if the RT was the same but PWSNVs were > 10. Unknown status was assigned if the RT was the same but PWSNVs were 3-10. RESULTS: 259 rCDI events were evaluable (50 [19.3%] reinfection; 198 [76.4%] relapse). The proportion of relapses was higher for ribotype 027 (84.5%) compared with other ribotypes (74.1%). Cumulative incidence of relapse was significantly lower for bezlotoxumab versus no bezlotoxumab (p < 0.0001), with a non-significant trend towards reduction for reinfection (p = 0.14). CONCLUSION: Bezlotoxumab treatment significantly reduced the rate of CDI relapse versus a regimen without bezlotoxumab. (NCT01241552/NCT01513239).

STOPGAP: a database for systematic target opportunity assessment by genetic association predictions.

SUMMARY: We developed the STOPGAP (Systematic Target OPportunity assessment by Genetic Association Predictions) database, an extensive catalog of human genetic associations mapped to effector gene candidates. STOPGAP draws on a variety of publicly available GWAS associations, linkage disequilibrium (LD) measures, functional genomic and variant annotation sources. Algorithms were developed to merge the association data, partition associations into non-overlapping LD clusters, map variants to genes and produce a variant-to-gene score used to rank the relative confidence among potential effector genes. This database can be used for a multitude of investigations into the genes and genetic mechanisms underlying inter-individual variation in human traits, as well as supporting drug discovery applications. AVAILABILITY AND IMPLEMENTATION: Shell, R, Perl and Python scripts and STOPGAP R data files (version 2.5.1 at publication) are available at https://github.com/StatGenPRD/STOPGAP . Some of the most useful STOPGAP fields can be queried through an R Shiny web application at http://stopgapwebapp.com . CONTACT: matthew.r.nelson@gsk.com. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

No evidence of large genetic effects on steroid response in asthma patients.

BACKGROUND: Inhaled corticosteroids (ICSs) are considered the most effective anti-inflammatory therapy for asthma control and management; however, there is substantial treatment response variability. OBJECTIVE: We sought to identify genetic markers of ICS response by conducting the largest pharmacogenetic investigation to date in 2672 ICS-treated patients with asthma. METHODS: Genotyping and imputation was performed in fluticasone furoate (FF) or fluticasone propionate-treated patients with asthma from 3 phase IIB and 4 phase IIIA randomized, double-blind, placebo-controlled, parallel group, multicenter studies. The primary end point analyzed was change in trough FEV1 (ΔFEV1) from baseline to 8 to 12 weeks of treatment. RESULTS: More than 9.8 million common genetic variants (minor allele frequency ≥ 1%) were analyzed to test for association with ΔFEV1. No genetic variant met the prespecified threshold for statistical significance. CONCLUSIONS: This study provides no evidence to confirm previously reported associations between candidate genetic variants and ICS response (ΔFEV1) in patients with asthma. In addition, no variant satisfied the criterion for genome-wide significance in our study. Common genetic variants are therefore unlikely to prove useful as predictive biomarkers of ICS response in patients with asthma.

Multi-population classical HLA type imputation.

Statistical imputation of classical HLA alleles in case-control studies has become established as a valuable tool for identifying and fine-mapping signals of disease association in the MHC. Imputation into diverse populations has, however, remained challenging, mainly because of the additional haplotypic heterogeneity introduced by combining reference panels of different sources. We present an HLA type imputation model, HLA*IMP:02, designed to operate on a multi-population reference panel. HLA*IMP:02 is based on a graphical representation of haplotype structure. We present a probabilistic algorithm to build such models for the HLA region, accommodating genotyping error, haplotypic heterogeneity and the need for maximum accuracy at the HLA loci, generalizing the work of Browning and Browning (2007) and Ron et al. (1998). HLA*IMP:02 achieves an average 4-digit imputation accuracy on diverse European panels of 97% (call rate 97%). On non-European samples, 2-digit performance is over 90% for most loci and ethnicities where data available. HLA*IMP:02 supports imputation of HLA-DPB1 and HLA-DRB3-5, is highly tolerant of missing data in the imputation panel and works on standard genotype data from popular genotyping chips. It is publicly available in source code and as a user-friendly web service framework.

Statistical assessment of biomarker replicability using MAJAR method.

In the era of precision medicine, many biomarkers have been discovered to be associated with drug efficacy and safety responses, which can be used for patient stratification and drug response prediction. Due to the small sample size and limited power of randomized clinical studies, meta-analysis is usually conducted to aggregate all available studies to maximize the power for identifying prognostic and predictive biomarkers. However, it is often challenging to find an independent study to replicate the discoveries from the meta-analysis (e.g. meta-analysis of pharmacogenomics genome-wide association studies (PGx GWAS)), which seriously limits the potential impacts of the discovered biomarkers. To overcome this challenge, we develop a novel statistical framework, MAJAR (meta-analysis of joint effect associations for biomarker replicability assessment), to jointly test prognostic and predictive effects and assess the replicability of identified biomarkers by implementing an enhanced expectation-maximization algorithm and calculating their posterior-probability-of-replicabilities and Bayesian false discovery rates (Fdr). Extensive simulation studies were conducted to compare the performance of MAJAR and existing methods in terms of Fdr, power, and computational efficiency. The simulation results showed improved statistical power with well-controlled Fdr of MAJAR over existing methods and robustness to outliers under different data generation processes. We further demonstrated the advantages of MAJAR over existing methods by applying MAJAR to the PGx GWAS summary statistics data from a large cardiovascular randomized clinical trial. Compared to testing main effects only, MAJAR identified 12 novel variants associated with the treatment-related low-density lipoprotein cholesterol reduction from baseline.

Deep sequencing of the LRRK2 gene in 14,002 individuals reveals evidence of purifying selection and independent origin of the p.Arg1628Pro mutation in Europe.

Genetic variation in LRRK2 predisposes to Parkinson disease (PD), which underpins its development as a therapeutic target. Here, we aimed to identify novel genotype-phenotype associations that might support developing LRRK2 therapies for other conditions. We sequenced the 51 exons of LRRK2 in cases comprising 12 common diseases (n = 9,582), and in 4,420 population controls. We identified 739 single-nucleotide variants, 62% of which were observed in only one person, including 316 novel exonic variants. We found evidence of purifying selection for the LRRK2 gene and a trend suggesting that this is more pronounced in the central (ROC-COR-kinase) core protein domains of LRRK2 than the flanking domains. Population genetic analyses revealed that LRRK2 is not especially polymorphic or differentiated in comparison to 201 other drug target genes. Among Europeans, we identified 17 carriers (0.13%) of pathogenic LRRK2 mutations that were not significantly enriched within any disease or in those reporting a family history of PD. Analysis of pathogenic mutations within Europe reveals that the p.Arg1628Pro (c4883G>C) mutation arose independently in Europe and Asia. Taken together, these findings demonstrate how targeted deep sequencing can help to reveal fundamental characteristics of clinically important loci.

Pharmacogenomics polygenic risk score for drug response prediction using PRS-PGx methods.

Polygenic risk scores (PRS) have been successfully developed for the prediction of human diseases and complex traits in the past years. For drug response prediction in randomized clinical trials, a common practice is to apply PRS built from a disease genome-wide association study (GWAS) directly to a corresponding pharmacogenomics (PGx) setting. Here, we show that such an approach relies on stringent assumptions about the prognostic and predictive effects of the selected genetic variants. We propose a shift from disease PRS to PGx PRS approaches by simultaneously modeling both the prognostic and predictive effects and further make this shift possible by developing a series of PRS-PGx methods, including a novel Bayesian regression approach (PRS-PGx-Bayes). Simulation studies show that PRS-PGx methods generally outperform the disease PRS methods and PRS-PGx-Bayes is superior to all other PRS-PGx methods. We further apply the PRS-PGx methods to PGx GWAS data from a large cardiovascular randomized clinical trial (IMPROVE-IT) to predict treatment related LDL cholesterol reduction. The results demonstrate substantial improvement of PRS-PGx-Bayes in both prediction accuracy and the capability of capturing the treatment-specific predictive effects while compared with the disease PRS approaches.

A statistical perspective on baseline adjustment in pharmacogenomic genome-wide association studies of quantitative change.

In pharmacogenetic (PGx) studies, drug response phenotypes are often measured in the form of change in a quantitative trait before and after treatment. There is some debate in recent literature regarding baseline adjustment, or inclusion of pre-treatment or baseline value as a covariate, in PGx genome-wide association studies (GWAS) analysis. Here, we provide a clear statistical perspective on this baseline adjustment issue by running extensive simulations based on nine statistical models to evaluate the influence of baseline adjustment on type I error and power. We then apply these nine models to analyzing the change in low-density lipoprotein cholesterol (LDL-C) levels with ezetimibe + simvastatin combination therapy compared with simvastatin monotherapy therapy in the 5661 participants of the IMPROVE-IT (IMProved Reduction of Outcomes: Vytroin Efficacy International Trial) PGx GWAS, supporting the conclusions drawn from our simulations. Both simulations and GWAS analyses consistently show that baseline-unadjusted models inflate type I error for the variants associated with baseline value if the baseline value is also associated with change from baseline (e.g., when baseline value is a mediator between a variant and change from baseline), while baseline-adjusted models can control type I error in various scenarios. We thus recommend performing baseline-adjusted analyses in PGx GWASs of quantitative change.

CauchyCP: A powerful test under non-proportional hazards using Cauchy combination of change-point Cox regressions.

Non-proportional hazards data are routinely encountered in randomized clinical trials. In such cases, classic Cox proportional hazards model can suffer from severe power loss, with difficulty in interpretation of the estimated hazard ratio since the treatment effect varies over time. We propose CauchyCP, an omnibus test of change-point Cox regression models, to overcome both challenges while detecting signals of non-proportional hazards patterns. Extensive simulation studies demonstrate that, compared to existing treatment comparison tests under non-proportional hazards, the proposed CauchyCP test (a) controls the type I error better at small α levels (<0.01); (b) increases the power of detecting time-varying effects; and (c) is more computationally efficient than popular methods like MaxCombo for large-scale data analysis. The superior performance of CauchyCP is further illustrated using retrospective analyses of two randomized clinical trial datasets and a pharmacogenetic biomarker study dataset. The R package CauchyCP is publicly available on CRAN.

VCSEL: PRIORITIZING SNP-SET BY PENALIZED VARIANCE COMPONENT SELECTION.

Single nucleotide polymorphism (SNP) set analysis aggregates both common and rare variants and tests for association between phenotype(s) of interest and a set. However, multiple SNP-sets, such as genes, pathways, or sliding windows are usually investigated across the whole genome in which all groups are tested separately, followed by multiple testing adjustments. We propose a novel method to prioritize SNP-sets in a joint multivariate variance component model. Each SNP-set corresponds to a variance component (or kernel), and model selection is achieved by incorporating either convex or nonconvex penalties. The uniqueness of this variance component selection framework, which we call VCSEL, is that it naturally encompasses multivariate traits (VCSEL-M) and SNP-set-treatment or -environment interactions (VCSEL-I). We devise an optimization algorithm scalable to many variance components, based on the majorization-minimization (MM) principle. Simulation studies demonstrate the superiority of our methods in model selection performance, as measured by the area under the precision-recall (PR) curve, compared to the commonly used marginal testing and group penalization methods. Finally, we apply our methods to a real pharmacogenomics study and a real whole exome sequencing study. Some top ranked genes by VCSEL are detected as insignificant by the marginal test methods which emphasizes formal inference of individual genes with a strict significance threshold. This provides alternative insights for biologists to prioritize follow-up studies and develop polygenic risk score models.