Isoform-level transcriptome-wide association uncovers extensive novel genetic risk mechanisms for neuropsychiatric disorders in the human brain

The isoTWAS framework

isoTWAS seeks to prioritize genes with transcript-isoforms of genes whose cis-genetic component of expression is significantly associated with a complex trait. We extend the traditional gene-level TWAS approach by jointly modeling the expression of distinct transcript-isoforms of a given gene as a matrix while accounting for the pair-wise correlations between these isoforms^{21, 37, 39, 40}. Here, we assume that local genetic variants directly modulate expression of an isoform. In addition, we assume that the abundance of a gene is measured as the is the sum of the abundance of its isoforms, computed as transcripts per million, or TPM (Supplemental Figure S1)^{32, 33, 41, 42}. Flexibly integrating isoform-level expression into trait mapping may lead to novel discoveries in disease mapping and prioritize specific isoforms that explain genetic associations. Accordingly, gene-level trait mapping using traditional TWAS methods may not necessarily detect a trait association on the gene-level if a gene has multiple isoforms but only one is associated with the trait (Figure 1a). This scenario may be particularly relevant in the human brain, where some genes may express up to hundreds to thousands of unique isoforms²⁵. By modeling the genetic architectures of isoforms of a gene simultaneously, isoTWAS provides a deeper understanding of potential transcriptomic mechanisms that underlie genetic associations.

• Download figure
• Open in new tab

Figure 1: Isoform-centric approach for complex trait mapping and prioritization of disease mechanisms at a genetic locus.

(a) Motivation for isoTWAS. Gene G has three isoforms associated with a trait but only one has an effect on the trait. Gene G itself does not show an association with the trait. Studying genetic associations with an isoform-centric perspective will prioritize Gene G, but not with a gene-centric perspective. (b) Schematic comparison of isoTWAS and TWAS. First, isoTWAS differs from TWAS by training a multivariate model of isoform expression, while TWAS models total gene expression, the sum of isoform expression. Second, isoTWAS maps isoform-trait associations through a step-wise hypothesis testing framework that provides gene-level false discovery control and isoform-level family-wide error rate control. TWAS only maps gene-trait association.

The isoTWAS framework contains three general steps (Figure 1). First, we build multivariate predictive models of isoform-level expression using well-powered functional genomics training datasets, including GTEx³⁷ and PsychENCODE^{24, 27}. Here, we trained and systematically compared 4 multivariate predictive frameworks: (1) multivariate elastic net penalized regression⁴³, (2) multivariate LASSO penalized regression with simultaneous covariance estimation (MRCE)⁴⁴, (3) multivariate elastic net regression with stacked generalization (joinet)⁴⁵, and (4) sparse partial least squares (SPLS)⁴⁶. As a baseline for comparison, we also modeled each individual isoform independently with univariate regularized regressions, as implemented in Gusev et al’s FUSION software^{4, 43, 47, 48} (see Methods and Supplemental Methods). Models were trained to predict isoform expression using the set of cis-SNPs within 1 Megabase (Mb) of the gene body (Methods, Figure 1b). Model performance was assessed via 5-fold cross-validation, using McNemar’s adjusted R² between observed and predicted expression.

Next, we use these models to impute isoform expression into an external GWAS cohort and quantify the association with the target GWAS phenotype (Figure 1b). If individual-level genotypes are available, isoform expression can be directly imputed as a linear combination of the SNPs in the models, and these associations can be estimated through appropriate regression analyses. If only GWAS summary statistics are available, imputation and association testing is conducted simultaneously through a weighted burden test⁴. Finally, isoTWAS performs step-wise hypothesis testing procedure to account for multiple comparisons and control for local LD structure. Isoform-level P-value are first aggregated to the gene-level to prioritize a gene using the Aggregated Cauchy Association Test (ACAT)⁴⁹, where false discovery rates are controlled, and then individual isoforms of prioritized genes are subjected to post-hoc family-wise error control⁵⁰ (Supplemental Figure S2, Methods). After this step, a set of isoforms are identified whose cis-genetic component of expression are associated with the trait of interest⁴. For these isoforms, we apply a rigorous permutation test whereby the SNP-to-isoform effects are permuted 10,000 times to generate a null distribution; this permutation test assesses how much signal is added by isoform expression, given the GWAS architecture of the locus, and controls for large LD blocks⁴. Lastly, we can isoform-level Bayesian fine-mapping at each locus with a significant trait association to identify the minimal credible set of isoforms that contains the ‘casual’ isoform at a 90% confidence level and to assign individual posterior inclusion probabilities (Figure 1b, Methods). isoTWAS is available as an R package through Github (https://github.com/bhattacharya-a-bt/isotwas).

isoTWAS improves prediction of isoform and gene expression in simulations and real data

Previous work has demonstrated that isoform-level quantifications, when propagated to the gene-level, can lead to more accurate gene expression estimates and differential expression results with short read RNA-seq data^{41, 42}. We therefore hypothesized that our multivariate SNP-based imputation of isoform expression, when summed to the gene-level, would outperform traditional gene-level TWAS models. To systematically evaluate the performance of TWAS versus isoTWAS models on prediction of total gene expression across a variety of genetic architectures, we conducted an extensive set of simulations across 22 different gene loci using European-ancestry data from the 1000 Genomes Project⁵¹ (Methods, Figure 2a). At each gene locus, we controlled gene expression heritability and simulated 2-10 distinct isoforms, varying the proportion of causal isoQTLs (p_causal) and their sharing between isoforms (p_shared). We then trained cross-validated, multivariate predictive models of isoform expression (isoTWAS) or univariate models of gene expression (TWAS). For isoTWAS, of the specific multivariate prediction models tested, multivariate elastic net⁴³ demonstrated the greatest CV prediction of isoform expression across most simulation settings (Figure 2b, Supplemental Figure S3, Supplemental Data 1). For gene expression prediction, the optimal isoTWAS models (in sum) outperformed the optimal TWAS model, particularly at sparser isoQTL architectures, with median absolute increase in adjusted R² of 0.6-3.5% (Figure 2c, Supplemental Figure S4, Supplemental Data 2). Performance gains decreased somewhat with denser isoQTL architectures, although real data is consistent with 0.1-1% sparsity (i.e., 1-10 causal e- and isoQTLs per gene or isoform)³⁷. In simulations, we found that isoTWAS prediction of gene expression also increases as the proportion of shared non-zero effect SNPs across isoforms decreases (Figure 2b-c, Supplemental Figure S4, Supplemental Data 2).

• Download figure
• Open in new tab

Figure 2: IsoTWAS models predict gene expression with more accuracy than TWAS models in simulated data.

(a) Simulation setup to generate isoform expression with specified isoQTL architecture, controlled expression heritability, number of isoforms, and inter-isoform correlation structure. (b) Proportion of simulations where the isoTWAS model has the maximum adjusted R² for marginal isoform prediction (Y-axis) across shared isoQTL proportion (X-axis), colored by isoTWAS method, facetted by causal isoQTL proportion (top margin) and proportion of isoform expression variance attributed to shared non cis-genetic effects (right margin). (c) Boxplots of difference in adjusted R² in predicting gene expression between isoTWAS and TWAS models from simulations with sample size 500 where isoform and gene expression heritability are set to 0.05, across causal isoQTL proportion (X-axis) and colored by number of transcripts per gene, facetted by proportion of shared isoQTLs (top margin) and proportion of variance explained by shared non cis-genetic effects (right margin).

Next, we assessed predictive performance in real data from 48 tissues (13 brain) with sufficient sample sizes (N > 100) in GTEx for all genes with multiple expressed isoforms (Supplemental Table S1; Methods). Altogether, for 48 tissues in GTEx, we built significant predictive models for 50,000 to 80,000 isoforms across 8,000 to 12,000 unique genes per tissue (Supplemental Figure S11, Supplemental Table S2). We considered 3 main criteria to evaluate the performance of both the multivariate and isoform-centric approaches of isoTWAS: (1) the number of isoforms whose expression can be imputed using multivariate/univariate models with cross-validation (CV) R² > 0.01; (2) the number of unique genes with at least one isoform that can be imputed at CV R² > 0.01; and (3) the number of unique genes in which total gene expression can be imputed at CV R² > 0.01 using isoTWAS (summed) or TWAS models.

At the isoform level (criterion 1), through joint multivariate modeling of isoform expression, we trained 2.3-2.5-fold more models at cross-validation (CV) R² > 0.01 across the 48 tissues, compared to traditional univariate approaches (Figure 3a, Supplemental Figure S5). Using these multivariate models, we improve prediction for 79-82% of isoforms with a median increase of ∼1.8-2.4 fold increase in adjusted R² (Supplemental Figure S6, Supplemental Table S2). Concordant with simulations, we found that the multivariate elastic net overwhelmingly outperformed other multivariate (and univariate) methods, indicating that leveraging the shared genetic architecture between isoforms of the same gene greatly aids in prediction of each individual isoform (Supplemental Figure S7, Supplemental Table S2). Notably, we observed that multivariate models were particularly important for brain tissues compared with non-brain tissues in GTEx, which showed significantly improved performance compared with univariate models (Figure 3b; P = 0.011 from OLS regression of median percent increase in CV R² for multivariate/univariate models against tissue type, adjusted for sample size), suggesting more shared isoQTL architecture in brain tissues than others which can be leveraged by isoTWAS for improved prediction. These gains in prediction accuracy directly translate into increased power in the trait association step⁵².

• Download figure
• Open in new tab

Figure 3: Multivariate isoform-level predictive models improve upon gene-level predictive models in predictive gene-level expression.

(a) Barplot showing the number of isoforms with CV R² > 0.01 (Y-axis) using multivariate (cream) and univariate (blue) modelling methods across brain tissues (X-axis). (b) Boxplot of median percent difference in predicting isoform expression (Y-axis) using multivariate compared to univariate method by brain and other tissue (X-axis). P-value corresponds to difference in median across groups, adjusted for sample size. (c) Barplot showing the number of genes passing CV thresholds (Y-axis) using TWAS (red) and isoTWAS (blue) across brain tissues. (d) Barplot showing the number of genes with CV R² > 0.01 (Y-axis) using TWAS (red) and isoTWAS (blue) across brain tissues (X-axis). (e) Boxplot of percent difference in R² (Y-axis) of predicting gene expression (isoTWAS – TWAS) in external datasets. X-axis shows the training and imputation datasets. The median percent difference in labelled.

At the gene level (criteria 2 and 3), isoTWAS also increased the number of genes with testable models in the trait mapping step and improved prediction of total gene expression. The number of unique genes with at least 1 isoTWAS model at CV R² > 0.01 (inclusion criterion for isoTWAS trait mapping) was 1.9-2.5 times larger than the number of unique genes with TWAS models achieving CV R² > 0.01 for gene expression prediction (Figure 3c, Supplemental Figure S8, Supplemental Table S2). For a given gene, isoTWAS models (summed) outperformed TWAS models in prediction of total gene expression by a median of 25-70% in cross validation (Supplemental Figure S9) with a 50-80% increase in the number of genes that are predicted at CV R² > 0.01 (Figure 3d, Supplemental Figure S10). We replicated these gains in total gene expression prediction using an independent, out of sample QTL reference panel of adult cortex from PsychENCODE/AMP-AD (Methods). Multivariate isoTWAS models outperformed univariate TWAS models in predicting total gene expression with a 15.2% median percent increase in adjusted R² when training in GTEx and testing in PsychENCODE/AMP-AD, and 23.9% vice versa; Figure 3e, Supplemental Table S3). As predictive performance is positively related to power to detect trait associations⁵¹, both the increased number and accuracy of trainable imputation models using isoTWAS’s multivariate predictive framework have strong implications on increased discovery in trait mapping⁵¹.

Performance of isoTWAS across distinct genomic contexts

As genes can differ widely with respect to the number and expression patterns of their constituent isoforms, as well as by other potentially relevant features such as gene length and SNP density, we next sought to characterize the impact of these factors of isoTWAS performance using GTEx models as evaluated by the 3 criteria outlined above (Methods, Supplemental Note, Supplemental Figures S12-19, Supplemental Data 3-4). Overall, we observed an increase in performance of isoTWAS multivariate modeling of both isoforms and genes with increasing number of isoforms per gene, although there was less conclusive of a pattern with increasing dominant isoform fraction⁵² (Supplemental Figures S12-13, Supplemental Note, Methods). We also noticed trends in the performance gain using isoTWAS multivariate modeling with respect to both isoform and gene expression prediction across gene length, SNP density at the gene locus, and sample size (Supplemental Figures S14-16, Supplemental Note). Finally, as the proportion of non-zero effect SNPs in the isoTWAS model that are shared across isoforms increased (Supplemental Note), we found an increasing trend in the gain in prediction of gene expression using isoTWAS compared to TWAS models (Supplemental Figure S17), reflecting a similar observation from simulation. Interestingly, as this proportion increased, we found an increase in the gain in prediction of isoform expression, reinforcing the utility of multivariate modelling in marginal prediction of isoform expression.

Lastly, we investigated how the robustness of isoform abundance estimation from short-read RNA-seq impacted performance gains of isoTWAS, compared with TWAS. Isoform abundance was initially quantified from probabilistic point estimates using Salmon, guided by Gencode annotations. We then assessed the performance of isoTWAS across loci binned by quantification variance measured across 50 inferential replicates from Salmon³². In general, we found that for isoform-level prediction, multivariate modeling in isoTWAS substantially outperformed univariate approaches as quantification variance increased. However, comparing isoTWAS isoform-centric and TWAS gene-centric models, there were no discernable trends in prediction of gene expression as the mean count and/or quantification variance of genes increased (Supplemental Figure S18-19, Supplemental Note). Finally, we evaluated the impact of reference transcriptome annotation fidelity by generating a synthetic dataset quantified using a reference annotation masking the dominant isoforms for a set of genes. As expected, performance of both isoTWAS and TWAS models declined when isoforms failed to be detected in expression quantification (Supplemental Note, Supplemental Figure S20).

Simulations support improved power and calibrated null across genetic architectures

We next introduced GWAS data for complex traits into our simulation framework to benchmark the false positive rate (FPR) and power of isoTWAS (Methods). First, we found that the FPR is controlled at 0.05 for isoform-level mapping using ACAT, consistent with gene-level mapping using TWAS (Supplemental Figure S21, Supplemental Data 4). For a simulated trait, we modeled causal effect architectures for a genomic locus with 2-10 isoforms under three main scenarios (Methods; Figure 4, Supplemental Figure S22): (1) where the true trait effect is from only total gene expression, (2) where there is only one isoform with a non-zero effect on the trait, called the “effect isoform”; and (3) where 2 isoforms are effect isoforms, with varying magnitudes and directions of association. This first scenario showed clear increases in power for TWAS over isoTWAS, but this advantage declined as the causal proportion of isoQTLs increased and as the proportion of shared isoQTLs increased (Figure 4a, Supplemental Data 5). For scenarios (2) and (3), as effects on the trait varied across isoforms of the same gene (Figure 4b-c, Supplemental Data 6-7), we observed large increases in power for isoTWAS over TWAS. Here, across most scenarios and causal effect architectures, isoTWAS demonstrated improved power compared with gene-level TWAS, particularly when only one isoform of a gene carried a trait effect or when two effect isoforms of the same gene had different directions of effects. However, when the effect sizes of two effect isoforms of the same gene converged, TWAS and isoTWAS demonstrated similar power to detect a gene-trait association (Figure 4c).

• Download figure
• Open in new tab

Figure 4: IsoTWAS improves power to detect gene-trait associations, especially when genetic effects differ across isoforms, in simulations.

(a) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10^-6, Y-axis) across causal proportion of isoQTLs (X-axis), facetted by proportion of shared isoQTLs (top margin) and proportion of variance explained by shared non cis-genetic effects (right margin). (b) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10^-6, Y-axis) across proportion of gene expression explained by effect isoform (X-axis). (c) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10^-6, Y-axis) across ratio of effect sizes for 2 effect isoforms. Across all plots, isoform and total gene expression heritability is set to 0.05 and causal proportion of 0.01 in (b-c). Points are colored by method and shaped by the number of isoforms per gene.

Finally, we assessed the performance of probabilistic fine-mapping in identifying the true effect isoform in our simulation framework of genes with 5 or 10 isoforms (Methods, Supplemental Figure S23, Supplemental Data 9). In general, the sensitivity of 90% credible sets (proportion of credible sets that contained the true effect isoform) was under-calibrated, likely due to difficulties in fine-mapping when QTL horizontal pleiotropy is extremely high⁵³. We found that the sensitivity of 90% credible sets decreased and the mean set size increased with increasing proportion of shared isoQTLs. Our simulation results suggest that varied isoQTL architectures and isoform-trait effects for isoforms of the same gene are key features that influence power gains in isoform-centric modeling.

isoTWAS increases discovery of trait associations across 15 neuropsychiatric disorders

To explore our central hypothesis that isoform-centric multivariate predictive modeling will improve discovery for complex trait mapping, particularly for brain relevant traits, we next sought to apply and directly compare isoTWAS and TWAS results across 15 neuropsychiatric disorders and traits. To maximize discovery, we trained both isoTWAS and TWAS models using a greatly expanded adult brain functional genomics reference panel (N = 2,115), comprised of frontal cortex samples from PsychENCODE and AMP-AD Consortia^{27, 54}, as well as using a developmental⁵⁵ pre-frontal cortex (N = 205) dataset (Methods; Figure 5). In the adult cortex, we trained models for 15,127 genes using isoTWAS passing the CV R² > 0.01 cutoff, compared with 14,283 genes using gene-level TWAS. In the developing cortex, despite a smaller sample size, isoTWAS models were successfully trained for 16,504 genes, compared with 10,535 genes using TWAS (Methods; Supplemental Table S1).

• Download figure
• Open in new tab

Figure 5: Isoform-level trait mapping increases discovery of genetic associations over gene-level trait mapping.

(a) Schematic diagram for trait mapping using gene-level TWAS and isoTWAS using PsychENCODE data. (b) Number of gene-trait associations overlapping GWAS risk SNPs within 0.5 Mb using gene-level TWAS (red), isoTWAS (blue), or either (green) at conservative permutation-based significance thresholds. (c) Scatterplot of standardized effect sizes (Z-scores) using isoTWAS and gene-level TWAS, considered associations with nominal P < 0.05 using both TWAS and isoTWAS. Gray line shows the 45-degree line and the green line shows an ordinary least squares regression. (d) Empirical Bayes estimate of test statistic inflation (X-axis) for TWAS (red) and isoTWAS (red) gene-level associations across 15 traits (Y-axis). (e) Mean percent increase in approximate X²-test statistic (squared Z-score), which is proportional to increase in effective sample size, for isoTWAS trait associations over TWAS trait associations. (f) Number of gene-trait associations overlapping GWAS risk SNPs within 0.5 Mb using splicing-event-based TWAS (purple), isoTWAS (blue), or either (green) at conservative permutation-based significance thresholds.

We next applied these models to perform trait associations using summary statistics from 15 neuropsychiatric and brain-related GWAS^56–70 (Methods, Figure 5a, Supplemental Figure S24), using the conservative stepwise hypothesis testing procedure (Figure 1a; FDR-adjusted P < 0.05 and within-locus permutation P_ACAT\< 0.05). We detected far more trait-associated genes with isoTWAS compared with TWAS, across adult (2,595 vs 1,589 genes) and developmental (4,062 vs 890 genes) reference panels, respectively (Supplemental Figure S25, Supplemental Data 10-13). In total across both reference panels and all 15 traits, isoTWAS detected 3,436 unique gene and 5,377 unique isoform-trait associations (Supplemental Figure S26). In addition, of the 1,335 genes with multiple isoform-trait associations, 661 gene exhibited distinct isoform-level associations in different directions.

We next sought to compare the performance of isoTWAS/TWAS in prioritizing candidate mechanisms within independent, high-confidence GWAS-significant loci⁷³. Across a combined 1,149 GWAS loci, isoTWAS identified significant associations within 323, compared with 201 detected by TWAS, a ∼60% increase in discovery (Figure 5b; Methods, Supplemental Table S4). For example, of the 287 GWAS loci identified for SCZ⁷⁴, isoTWAS prioritized genes within 70 and 86 unique loci across adult and developmental reference panels, respectively, compared with 56 and 29 loci for TWAS (Figure 5b). Further, across the 15 traits, 96% of gene-level TWAS associations (193/201) were concordantly identified by isoTWAS. Likewise, the standardized effect sizes for significant gene- and isoform-level associations were highly correlated (r = 0.84, P < 2.2 × 10^-16; Figure 5c). Finally, to explore whether these isoTWAS-specific associations were capturing true disease signal, we compared the rate at which each method prioritized constrained genes (probability of loss-of-function intolerance, pLI ≥ 0.9; Supplemental Tables S5-S8), which are known to be substantially enriched for disease associations⁷⁵. Across adult and developmental panels for 15 traits, respectively, isoTWAS prioritized 724 and 385 constrained genes compared with 106 and 200 with TWAS, a significant increase (adult: P=0.048, developmental: P = 1.23 × 10^-5, Fisher’s Exact test). Altogether, these results emphasize that isoTWAS not only recovers the vast majority of TWAS associations but also greatly increases discovery of candidate GWAS mechanisms, particularly for genes intolerant to protein-truncating variation⁷⁶.

In our evaluation of methods in real data, we sought to compare discovery using isoTWAS to discovery using local splicing-event based trait mapping. Briefly, for the developmental brain dataset, we first calculated intron usage using Leafcutter⁷⁷ and transformed these usage percentages to M-values⁷⁸. Then, for all introns mapped to a given gene, we used all SNPs within 1 Mb of a splicing-event to predict its usage and mapped trait associations for these splicing events using isoTWAS’s multivariate framework (Methods). Overall, when aggregated to the gene-level, across 15 traits, we found that isoTWAS prioritized features at ∼40% more independent GWAS loci (167 loci) than splicing-event based trait mapping (119 loci), with 108 loci (90.7%) jointly identified (Figure 5f), using the same developmental brain reference panel. Taken together, isoTWAS’s specific focus on modeling isoforms of a gene provided considerable gains over considering only total gene expression or intron usage in identifying trait associations for genes and their transcript-isoforms.

To investigate whether this increase in trait mapping discovery reflected true biological signal, rather than test statistic inflation due to the increased number of tests (∼4-fold increase in number of tests), we next compared the null distributions across methods for results across the 15 traits (Supplemental Figure S27). As the genomic inflation factor is not a reliable measure in TWAS settings⁷⁹, we estimated inflation in gene-level test statistics using an empirical Bayes approach (Methods). Collapsing across all 15 traits, there were no significant differences between TWAS and isoTWAS in the 95% credible intervals for test statistic inflation (Figure 5d). Using a heuristic to estimate increases in effective sample size (Methods), we observed an approximate increase in effective sample size of 10-20% when using isoTWAS compared to TWAS (Figure 5e, Supplemental Table S9). These analyses indicate that isoTWAS discovery is both well-calibrated to the null and facilitates increased discovery in real data compared to gene-level TWAS.

We also empirically compared probabilistic fine-mapping⁵⁵ of results from isoTWAS and gene-level TWAS (Methods). Here, we conducted fine-mapping on significant trait-associated genes/isoforms (adjusted P < 0.05 and permutation P < 0.05) and are within 1 Mb of one another; we term a locus with overlapping genes/isoforms a risk region. Overall, the mean number of genes in a risk region using TWAS was 3.15 compared to 3.90 using isoTWAS (Supplemental Figure 28a); the mean number of genes in a 90% credible set using TWAS was 1.33 compared to 1.25 using isoTWAS (Supplemental Figure 28a). On average, there were 1.54 isoforms per gene in a risk region and 1.27 isoforms per gene in a 90% credible set (Supplemental Figure 28b). Isoform-centric modeling presents unique challenges for fine-mapping due to potentially high levels of horizontal pleiotropy, and remains an important and open question for the field. Nevertheless, isoTWAS identified a relatively comparable number of genes in risk regions compared with TWAS, and the combination of conservative two-step trait mapping, permutation testing, and probabilistic fine-mapping were critical for maintaining a narrow credible set size.

isoTWAS identifies biologically-meaningful trait associations undetectable by TWAS

Overall, isoTWAS prioritized dozens of novel candidate risk genes and mechanisms in the developing and adult brain for 15 neuropsychiatric traits. These isoTWAS-prioritized genes were enriched for multiple relevant pathways consistent with the biology of the underlying trait, including cell proliferation pathways for brain volume, calcium channel activity for schizophrenia and neuroticism, and proteasome regulation in Alzheimer’s Disease, among others (Supplemental Figure S29). We highlight below several examples of novel trait associations in which isoTWAS prioritized a highly constrained gene within a GWAS locus (Supplemental Table S5-8). For example, we identified an association for SCZ with developmental expression of ENST00000519133, an isoform of SNAP91 (adjusted P = 6.06 × 10^-7 in isoTWAS, pLI = 0.99, chromosomal location of 6q14.2), with a GWAS-significant variant rs217291 within the gene body. SNAP91 has been predicted to affect clathrin and phosphatidylinositol biding activity and synaptic vesicle recycling and has previously shown to impact synaptic development⁸⁰. In addition, we found that imputed developmental cortex expression of ENST00000476671, an isoform of KMT2E, was associated with CDG risk (adjusted P = 7.63 × 10^-3, pLI = 1, chromosomal location of 7q22.3); the GWAS SNP rs2385537, associated in a meta-analysis of ADHD, ASD, BP, and SCZ⁸¹, is within the gene body. KMT2E regulates post-translational histone methylation of histone 3 on lysine 4A, and KMT2E heterozygous variants are associated with risk of neurodevelopmental disorders⁸². Imputed adult brain expression of ENST00000492146 (adjusted P = 1.14 × 10^-6, pLI = 0.94, chromosomal location of 3p21.1), an isoform of SFMBT1, showed a strong association with SCZ risk and contains a GWAS SNP within the gene body (rs2071044). In a recent gene-based analysis of GWAS data for SCZ and BD, decreased expression of SFMBT1 was associated with increased risk of both disorders⁸³, consistent with the effects for the isoform we identify in this isoTWAS analysis. Lastly, our analysis implicated ENST00000537270, an isoform of KMT5A (pLI = 0.99, found in 12q24.31), in an association with SCZ risk. KMT5A, a H4K20 methyltransferase, has been previously implicated in GWAS for SCZ but cis-eQTLs of the gene have not be colocalized with the GWAS signal previously^{74, 84, 85}.

As evidenced by our simulation analyses, a main advantage of isoTWAS over TWAS is the identification of trait associations for isoforms of genes, where the gene itself is not associated with the trait. We illustrate several examples of isoforms prioritized by isoTWAS, all in adult cortex, for genes with limited or distinct expression QTLs (Figure 6, Supplemental Figure S30, Supplemental Data 14). First, we detected a strong SCZ association with ENST00000492957, an isoform of AKT3 (pLI = 1) at 1q43-q44, which encodes a serine/threonine-protein kinase that regulated many processes like metabolism and cell growth, proliferation, and survival. AKT3 has shown effects on anxiety, spatial and contextual memory, and fear extinction in mice, such that loss-of-function of AKT3 causes learning and memory deficits^{86, 87}. Within the GWAS locus, there was a strong overlapping isoQTL signal (P < 10^-50); however, the eQTL signal includes only a single SNP that reaches P < 10^-6, which is in low LD with the GWAS-significant SNPs in the locus (Figure 6a). The lead isoQTL (rs4430311) showed a significant, negative association with ENST00000492957, but not AKT3 expression, with increasing SCZ risk alleles. However, this SNP has only a nominally significant positive association with AKT3, as the number of alternative alleles at the SNP increased. Interestingly, a different isoform of AKT3 (ENST00000681794) was prioritized in an association with BV, which also has a GWAS hit at this locus (Supplemental Figure S30). The two distinct isoforms of AKT3 have distinct 3’ transcript structures, particularly close to the lead isoQTL of ENST00000681794. These results suggest a complex role of AKT3 isoforms with brain-related traits to be explored further computationally and experimentally.

• Download figure
• Open in new tab

Figure 6: isoTWAS implicates isoforms of AKT3, CUL3, HSPD1, and PCLO in genetic associations with psychiatric traits.

In each plot, (top) Manhattan plots of GWAS risk, total gene eQTLs, and isoQTLs, colored by LD to the lead isoQTL and lead isoQTL shown in a triangle and labelled. LD is based on the 1000 Genomes European reference. (bottom) Boxplots of gene and isoform expression by genotype of the lead isoQTL SNP and forest plot of the lead isoQTL’s effect and 95% confidence interval on the trait, gene, and isoform, with P-values labelled. (a) SCZ risk, AKT3 gene expression, and ENST00000492957 isoform expression. (b) SCZ risk, CUL3 gene expression, ENST00000409096 isoform expression. (c) SCZ risk, HSPD1 gene expression, ENST00000678969 isoform expression. (d) CDG risk, PCLO gene expression, ENST00000423517 isoform expression.

Similarly, we found a strong isoQTL signal for ENST00000409096 but a weak eQTL signal of CUL3 in the 2q36.2 locus (pLI = 0.99), in another association with SCZ (Figure 6b). CUL3 is a component of Cullin-RING E3 ubiquitin ligase complexes, involved in protein ubiquitylation, cell cycle regulation, protein trafficking, signal transduction, and transcription. Previous work has implicated CUL3 dysregulation as a pathological mechanism for both SCZ and ASD risk⁸⁸. Next, an isoform ENST00000678969 of HSPD1, encoding a mitochondrial heat shock protein, was associated with SCZ risk (adjusted P = 8.79 × 10^-12, pLI = 0.99, 2q33.1) and showed a similar pattern across GWAS, eQTL, and isoQTL signals (Figure 6c). HSPD1 was identified in 3 independent secondary analyses of SCZ GWAS data, is among multiple non-MHC immune genes implicated in SCZ, and has roles in brain hypomyelination⁸⁹. Our work adds to these previous studies by implicating a specific isoform of HSPD1 at the locus. Lastly, ENST00000423517, an isoform of PCLO, was associated with multiple traits in the CDG analysis (meta-analysis of ADHD, BP, MDD, and SCZ; adjusted P = 1.83 × 10^-4, pLI = 1). Again, we found a strong isoQTL but not eQTL signal, with the CDG risk allele negatively associated with isoform expression. PCLO is involved in the presynaptic cytoskeletal matrix, establishing active synaptic zones, and synaptic vesicle trafficking; rare variants of PCLO in diverse populations have been recently implicated in risk of SCZ and ASD^{90, 91}. Altogether, these results highlight the importance of incorporating isoform-level regulation for prioritizing novel candidate GWAS risk mechanisms, as implemented in our isoTWAS framework.

Training predictive models of isoform expression

Trait association and step-wise hypothesis testing

We note that the tests of association in isoTWAS are similar to tests in differential transcript expression analyses, as TWAS tests of association are analogous to tests in differential gene expression analyses. isoTWAS and TWAS are distinct, as these methods consider imputed isoform and gene expression, respectively, as predicted by the trained expression models. If individual-level genotypes are available in the external GWAS panel, isoform expression can be directly imputed by multiplying the SNP weights from the predictive model with the genotype dosages in the GWAS panel. If only summary statistics are available, we adopt the weighted burden test from Gusev et al with an ancestry-matched linkage disequilibrium panel^{4, 98}. Compared to TWAS, a main distinction for isoTWAS association testing is the increased number of tests (approximately 4-fold the number of isoforms than genes)²¹ and the potential correlation in test statistics for isoforms of the same gene.

Accordingly, we perform a two-step hypothesis testing framework (Supplemental Figure S2)¹⁰². In the first step, for every isoform with a trained model, we generate a TWAS test statistic using either linear regression for GWAS with individual-level genotypes or the weighted burden test for GWAS with only summary statistics⁴. Given the t test statistics T₁, …, T_t for a single gene, we used an omnibus test to aggregate the t test statistics into a single P-value for a gene. We benchmark different omnibus tests in simulations, but the default omnibus test in isoTWAS is aggregated Cauchy association test (ACAT)⁴⁹. We control for false discovery across all genes via the Benjamini-Hochberg procedure, but the Bonferroni procedure can also be applied for more conservative false discovery control. In the second step, for isoforms for genes with an adjusted omnibus P < 0.05, we employ Shaffer’s modified sequentially rejective Bonferroni (MSRB) procedure to control the within-gene family-wide error rate. At the end of these two steps, we identify a set of genes and their isoforms that are associated with the trait.

Control for false positives within GWAS loci

In TWAS and related methods, association statistics have been shown to be well calibrated under the null of no GWAS association. However, within loci harboring significant GWAS signal, false positive associations can result when eQTLs and GWAS coincide within overlapping LD blocks. To address this, we have adopted two conservative approaches to control for type 1 error within GWAS loci, namely (1) permutation testing and (2) probabilistic fine mapping. The permutation testing approach, adopted from Gusev et al⁴, is a highly conservative test of the signal added by the SNP-transcript effects from the predictive models, conditional on the GWAS architecture of the locus. Briefly, we shuffle the SNP-transcript effects in the predictive models 10,000 times and generate a null distribution for the TWAS test statistic for each isoform. We then use this null distribution to generate a permutation-based P-value for the original test statistic for each isoform. Additionally, we aggregate these null distributions using an omnibus test and compare the omnibus P-value to this distribution to generate permutation P-value for each gene. Finally, we can employ isoform-level probabilistic fine-mapping using methods from FOCUS⁵⁵ to generate credible set of isoforms that explain the trait association at a locus. We only run isoform-level fine-mapping for significantly-associated isoforms in overlapping 1 Megabase windows.

Simulation framework

We adopt techniques from Mancuso et al’s twas_sim protocol¹⁰³ to simulate multivariate isoform expression based on randomly simulated genotypes and environmental random noise. First, for n samples, we generate a matrix of genotype dosages for the SNPs within 1 Megabase of 22 different genes (1 per chromosome) using an LD reference panel of European subjects from 1000 Genomes Project⁵¹.

Next, we generate a matrix of SNP-isoform effects across different causal SNP proportions p_c, numbers of isoforms t, and p_s proportion of the SNP-isoform effects being shared across isoforms of the same gene. We then add two matrices of random noise U and ϵ. The first matrix U noise represents non cis-genetic effects on isoforms that are correlated between samples and isoforms; we control the proportion of variance explained in isoform expression attributed to U using a parameter ᓂ_h. The second matrix ϵ is a matrix of random noise that is independent for each isoform, such that where . We generate 10,000 simulations for each configuration of the simulation parameters, varying , and . Full mathematical details are provided in Supplemental Methods and summarized in Figure 2.

We also generate traits under three distinct scenarios, with a continuous trait with heritability and a GWAS sample size of 50,000 (Supplemental Methods):

1. Only gene-level expression has a non-zero effect on trait. Here, we sum the isoform expression to generate a simulated gene expression. We randomly simulate the effect size and scale the error to ensure trait heritability.

2. Only 1 isoform has a non-zero effect on the trait. Here, we generate a multivariate isoform expression matrix with 2 isoforms and scale the total gene expression value such that one isoform (called the effect isoform) makes up p_g ∈ {0.10, 0.30, 0.50, 0.70, 0.90} proportion of total gene expression. We then generate effect size for one of the isoforms and scale the error to ensure trait heritability.

3. Two isoforms with different effects on traits. Here, we generate a multivariate isoform expression matrix with 2 isoforms that make up equal portions of the total gene expression. We then generate an effect size of α for one isoform and p_eα for the other isoform, such that p_e ∈ {−1, −0.5, −0.2,0.2,0.5,1}. We then scale the error to ensure trait heritability.

Lastly, to estimate the approximate false positive rate (FPR), we followed the same simulation framework to generate eQTL data and GWAS data. In the GWAS data, however, we set the effect of gene- and isoform-level imputed expression to 0 to generate a simulated trait under the null, where the gene and isoforms are not associated with the trait. We then estimate the FPR by calculating the proportion of gene-trait associations at P < 0.05 under this null across 20 sets of 1,000 simulated GWAS panels. We also assess isoform-level fine-mapping using FOCUS in a scenario with a gene with 5 or 10 isoforms and a single effect isoform. We compute the sensitivity of 90% credible sets of isoforms (proportion of credible sets that contain the effect isoform) and the number of isoforms in the 90% credible set.

GTEx processing and model training

We quantified GTEx v8³⁷ RNA-Seq samples for 48 tissues using Salmon v1.5.2³² in mapping-based mode. We first built a Salmon index for a decoy-aware transcriptome consisting of GENCODE v38 transcript sequences and the full GRCh38 reference genome as decoy sequences²¹. Salmon was then run on FASTQ files with mapping validation and corrections for sequencing and GC bias. We computed 50 inferential bootstraps for isoform expression using Salmon’s Expectation-Maximization algorithm. We then imported Salmon isoform-level quantifications and aggregated to the gene-level using tximeta⁴¹.

Using edgeR, gene and isoform-level quantifications underwent TMM-normalization, followed by transformation into a log-space using the variance-stabilizing transformation using DESeq2^{104, 105}. We then residualized isoform-level and gene-level expression (as log-transformed CPM) by all tissue-specific covariates (clinical, demographic, genotype PCs, and expression PEER factors) used in the original QTL analyses in GTEx. We calculated the quantification variance across inferential replicates using the computeInfRV() function from the fishpond package¹⁰⁶. We computed the isoform fraction using the isoformToIsoformFraction() function from the IsoformSwitchAnalyzeR package⁵⁴.

SNP genotype calls were derived from Whole Genome Sequencing data for samples from individuals of European ancestry, filtering out SNPs with minor allele frequency (MAF) less than 5% or that deviated from Hardy-Weinberg equilibrium at P < 10^-5. We further filtered out SNPs with MAF less than 1% frequency among the European ancestry samples in 1000 Genomes Project⁵¹.

Details of the model training pipeline for GTEx are similar to those summarized in Supplemental Figure S10. Gene-level univariate models were trained using elastic net regression⁴³, best linear unbiased predictor (BLUP) in a linear mixed model⁴⁸, and SuSiE⁴⁷, using all SNPs within 1 Mb of the gene body^{4, 43, 47, 101}. For each gene, the best performing model was chosen based on McNemar’s adjusted 5-fold cross-validation (CV) R². We selected only genes with CV R² ≥ 0.01. We applied multivariate modeling outlined in isoTWAS to train isoform-level predictive models, selecting only those isoform models with CV R² ≥ 0.01 (Supplemental Figure S2). All isoTWAS models generated are publicly available (see Data Availability).

Developmental brain reference panel processing and model training

We quantified developmental frontal cortex²⁴ (N = 205) RNA-Seq samples using Salmon v1.8.0³² in mapping-based mode. We used the same indexed transcriptome as in the GTEx analysis and ran Salmon with mapping validation and corrections for sequencing and GC bias. We computed 50 inferential bootstraps for isoform expression using Salmon’s Expectation-Maximization algorithm. We then imported Salmon isoform-level quantifications and aggregated to the gene-level using tximeta⁴¹. Using edgeR, gene and isoform-level quantifications underwent TMM-normalization, followed by transformation into a log-space using the variance-stabilizing transformation using DESeq2^{104, 105}. We then residualized isoform-level and gene-level expression (as log-transformed CPM) by covariates (age, sex, 10 genotype PCs, 90 and 70 hidden covariates with prior, respectively). Typed SNPs with non-zero alternative alleles, MAF > 1%, genotyping rate > 95%, HWE P < 10^-6 were first imputed to TOPMed Freeze 5 using minimac4 and eagle v2.4^{107, 108}. We then retained biallelic SNPs with imputation accuracy R² > 0.8, with rsIDs. Finally, we filtered out SNPs with MAF < 0.05 or that deviated from Hardy-Weinberg equilibrium at P < 10^-6.

Adult brain reference panel processing and model training

Matched genotype and RNAseq data from adult brain cortex tissue from n = 2,365 individuals were compiled and processed from the PsychENCODE Consortium²² and the Accelerating Medicines Partnership Program for Alzheimer’s Disease (AMP-AD)⁵⁶, consisting of the individual studies BipSeq, BrainGVEX, CommonMind Consortium (CMC), CommonMind Consortium’s National Institute of Mental Health Human Brain Collection Core (CMC HBCC), Lieber Institute for Brain Development-szControl (LIBD_szControl), UCLA-ASD, Religious Orders Study and the Memory and Aging Project (ROSMAP), Mount Sinai Brain Bank (MSBB) and MayoRNAseq.

Typed genotypes were lifted over to the GRCh38 build using CrossMap v.0.6.3¹⁰⁹ and then filtered to remove variants where the reference allele matched any of the alternate alleles. Genotype data from whole genome sequencing (BrainGVEX, UCLA-ASD, ROSMAP, MSBB, and MayoRNAseq) was further filtered to variants present on the Infinium Omni5-4 v1.2 array in order to satisfy the imputation server’s maximum limit of 20,000 typed variants per 20Mb. All genotype data was further processed with PLINK v1.90b6.21¹¹⁰, removing variants with HWE P < 10^-6, MAF< 0.01 or missingness rate > 0.05, and removing samples with missingness rate > 0.1 across typed variants or missingness rate > 0.5 on any individual chromosome. Genotype data was prepared for imputation using the McCarthy Group’s HRC-1000G-check-bim-v4.3.0 tool against freeze 8 of the Trans-Omics for Precision Medicine (TOPMed) reference panel¹¹¹. The tool removes A/T and G/C SNPs with MAF > 0.4, variants with alleles that differ from the reference panel, variants with an allele frequency difference > 0.2 from the reference panel and variants not in the reference panel. Additionally, the tool updates strand, position and reference/alternate allele assignment to match the reference panel.

Genotypes were then passed into the TOPMed Imputation Server by individual array batch¹¹². The genotypes were phased with Eagle v2.4 and imputed with Minimac4 using the TOPMed reference panel^{107, 108}. Further QC was performed on the imputed genotypes using bcftools v1.11 and PLINK. The imputed genotypes were filtered to well-imputed variants with R2 > 0.8. The arrays were merged after filtering to variants that were well-imputed in all arrays to be merged. Only arrays with at least 400,000 variants after pre-imputation QC were merged in order to prevent too many variants from dropping out. The merged genotype data was then converted to PLINK 1 binary format and further processed with PLINK, removing variants with duplicates, HWE P < 10^-6, MAF < 0.01 or missingness rate > 0.05 and removing samples with missingness rate > 0.1. Samples from the same individual were identified by calculating the genetic relatedness matrix using SnpArrays.jl and finding sets of samples with relatedness > 0.75. From each set of replicates, only the genotyped sample from the array with the most variants after pre-imputation QC was kept. For model training, only SNPs annotated in HapMap3 were retained¹¹³.

RNAseq paired reads from each study were sorted by name and then converted to FASTQ format using samtools v1.14¹¹⁴. The reads were then quantified using salmon v1.8.0 in mapping-based mode using a full decoy indexed from GENCODE v38 transcriptome and GRCh38 patch 13 assembly³². Quantification was run using standard EM algorithm with library type automatically inferred and estimates adjusted for sequence-specific and fragment-level GC biases. Bootstrapped abundance estimates were calculated using 50 bootstrap samples. Isoform-level expression was summarized to the gene-level using tximeta¹¹⁵. Only isoforms with 0.1 TPM for more than 75% of samples were retained. The resulting expression was normalized using the variance-stabilizing transformation from DESeq2¹⁰⁵. Samples with WGCNA network connectivity scores of less than -3 were removed as outliers, resulting in a total of 2,115 samples¹¹⁶. Isoform- and gene-level expression was then batch-corrected using ComBat, using study site as the batch¹¹⁷. Lastly, age, sex, squared age, 10 genotype PCs, and 200 (for gene expression) and 175 (for isoform expression) hidden covariates with prior (using sequencing metrics from PicardTools to estimate the prior) were removed from the expression^{118, 119}. The number of hidden covariates with prior (HCP) were selected by optimizing the number of nominal cis-gene-level expression QTLs and cis-isoform-level expression QTLs at Bonferroni-corrected P < 0.05, respectively, on a grid from 100 to 300 HCPs, as detected by QTLtools¹²⁰.

Details of the model training pipeline are summarized are equivalent to those used to train models in GTEx data. All isoTWAS models generated are available at https://zenodo.org/record/679594792.

Gene- and isoform-level trait mapping

We conducted gene- and isoform-level trait mapping for 15 neuropsychiatric traits: attention-deficit hyperactivity disorder (ADHD, N_cases = 20,183/N_controls = 35,191)⁵⁹, Alzheimer’s disease (ALZ, 90,338/1,036,225)⁶⁰, anorexia nervosa (AN, 16,992/55,525)⁷², autism spectrum disorder (ASD, 18,381/27,969)⁵⁸, bipolar disorder (BP, 41,917/371,549)⁶¹, brain volume (BV, N = 47,316)⁶², cross-disorder (CDG, 232,964/494,162)⁶³, cortical thickness (CortTH, N = 51,665)⁶⁴, intracranial volume (ICV, N = 32,438)⁶⁵, major depressive disorder (MDD, 246,363/561,190)⁶⁶, neuroticism (NTSM, N = 449,484)⁶⁷, obsessive compulsive disorder (OCD, 2,688/7,037)⁶⁸, panic and anxiety disorders (PANIC, 2,248/7,992)⁶⁹, post-traumatic stress disorder (PTSD, 32,428/174,227)⁷⁰, and schizophrenia (SCZ, 69,369/236,642)⁷¹. For gene-level trait mapping, we used the weighted burden test, followed by the permutation test, as outlined by Gusev et al⁴. For isoform-level trait mapping, we used the stage-wise testing procedure outlined in the isoTWAS method. In-sample LD from the QTL reference panels was used to calculate the standard error in the weighted burden test. For isoforms, irrespective of their corresponding genes, passing both stage-wise tests and the permutation test, we employed isoform-level probabilistic fine-mapping using FOCUS with default parameters⁵⁵. These methods are summarized in Supplemental Figure S24.

We estimated the percent increase in effective sample size by employing the following heuristic. We convert gene-level association P-values into χ² test statistics with 1 degree-of-freedom. For χ² > 1, we then calculate the percent increase for isoTWAS-based associations versus TWAS-based associations. As the mean of the χ² distribution is linearly related to power and sample size¹²¹, we can use this percent increase in test statistic as a measure of power or effective sample size. We defined independent genome-wide significant SNPs in GWAS by LD clumping with lead GWAS SNP < 5 × 10^-8 with P-value used for ranking and a R² threshold of 0.2.