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

The isoTWAS pipeline

The isoTWAS framework prioritizes transcript-isoforms of genes in which the cis-genetic component of expression is significantly associated with a complex trait. Expression of different isoform transcripts of a gene can be modelled as a multivariate data object, with, on average, four GENCODE v40-annotated isoforms per gene (mean 4.05, standard deviation 7.28, median 1), with strong pair-wise correlations between these isoforms^33,35 Here, we assume that local genetic variants directly modulate expression of an isoform. In addition, we assume that the abundance of a gene (total gene expression) is the sum of the abundance of its isoforms (isoform expression) (Supplemental Figure S1), on the raw count scale. 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 associations. The advantages of isoTWAS over TWAS are summarized.

The isoTWAS pipeline contains three general steps. First, we build multivariate predictive models of isoform-level expression using well-powered functional genomics training datasets, such as GTEx³³ and PsychENCODE^23,26. Here, we trained and systematically compared 6 multivariate predictive frameworks: (1) multivariate elastic net penalized regression³⁶, (2) multivariate LASSO penalized regression with simultaneous covariance estimation³⁷, (3) multivariate Sum of Single Effects (SuSiE)³⁸, (4) Breiman and Friedman’s curds-and-whey procedure³⁴, (5) ordinary least squares with clustered standard errors³⁹, and (6) univariate regularized regressions of each individual isoform, as implemented in Gusev et al’s FUSION software^4,36,38,40 (see Methods and Supplemental Methods). Models were trained to predict the expression of heritable isoforms 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 between imputed isoform expression and 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⁴. isoTWAS uses a two-step hypothesis testing procedure to account for multiple comparisons: isoform associations 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, to identify a set of isoforms that explains the genetic association with the trait at the locus, probabilistic isoform-level fine-mapping can be applied (Figure 1b, Methods).

By considering isoforms of a gene as a multivariate data structure and aggregating isoform-level tests to the gene-level, isoTWAS boosts both prediction of gene expression and power to detect trait associations, especially when genetic effects on isoforms are different and isoforms have different effects on a phenotype. Not only does isoTWAS provide substantial power gains, as we show through simulations and empiric data, this last step also provides a testable hypothesis for the mechanism underlying a trait association (Figure 1b). isoTWAS is available as an R package through Github (https://github.com/bhattacharya-a-bt/isotwas).

isoTWAS improves prediction of total gene expression in 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^43,44. Using GTEx data from 13 brain tissues and 35 other tissues with sufficient sample sizes (N > 100, Supplemental Table S1; Methods), we demonstrate that, by accounting for isoform-level genetic architecture, isoTWAS increases prediction of total gene expression over TWAS by a median of 15-35%, both in- and out-sample (Figure 2). As the trait mapping step in TWAS and isoTWAS only tests models with a baseline prediction cross-validation R² > 0.01, this increase in accuracy leads to nearly 2-fold more testable gene models under the isoTWAS framework. Multivariate SuSiE provided the models with best prediction of isoform expression in 40-55% of models (Supplemental Figure S3), with multivariate regression with simultaneous covariate estimation showing the best prediction in 20-30% of models, underscoring the utility of modeling the covariance structure between isoforms to improve prediction.

• Download figure
• Open in new tab

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

(a) Schematic diagram for predictive modeling and comparison to gene-level expression using gene-level TWAS and isoTWAS models. (b) Barplot showing the number of genes passing heritability (error bar) and cross-validation (CV) cutoffs (bar). TWAS heritability (or CV) cutoff is nominal P < 0.05 (CV R² > 0.01). isoTWAS heritability (or CV) cutoff is a gene with at least 1 isoform that is heritabie at P < 0.05 (predictable at CV R² > 0.01) (c) For 13 GTEx brain tissues, the percent difference in CV R² of predicting total gene expression using isoTWAS versus TWAS models. The median percent increase is labelled and the red line marks 0% difference. (d) Using completely independent GTEx frontal cortex and adult PsychENCODE frontal cortex models, the percent difference in external predictive R² using either GTEx or PsychENCODE as the training sample. (e) Using 48 GTEx tissues, ratio of the number of cis-heritable isoforms over the number of cis-heritable genes in brain tissues (purple) and other tissues (brown). F-test shows the difference in means between the two groups, adjusting for sample size for the tissue.

The advantage for isoTWAS starts at the feature selection step, in which imputation models are filtered based on nominally significant positive cis-heritability (P < 0.05, Methods). Across the 13 brain tissues, there are nearly 3 to 4-fold more the number of genes with at least 1 cis-heritable isoform, compared to the number of genes that are cis-heritable themselves (mean of 6331 genes with at least 1 cis-heritable isoform vs. 1831 genes that are cis-heritable across 13 tissues, Figure 2b). When we consider genes that pass the cis-heritability feature selection cutoff, we trained twice as many isoTWAS models as TWAS models with CV R² ≥ 0.01 (Figure 2b). These trends extend for nearly all other tissues as well, with greater than 2.5-fold more cis-heritable and predictable genes on the isoform-level than total gene expression level across GTEx tissues (Supplemental Figure S4-5, Supplemental Table S2). The median ratio of numbers of heritable isoforms to genes was significantly larger for brain tissues compared to others (3.40 vs 3.08, respectively; P = 0.004), when accounting for sample size (Figure 2e; Supplemental Figure S8). Notably, the least complex isoform architecture was exhibited by whole blood samples, with only a 2:1 ratio of cis-heritable isoforms to genes. These results may reflect the more complex splicing patterns in the brain compared to other tissues and suggest that blood-based gene imputation panels may generalize poorly to the complex gene regulatory mechanisms present in the brain.

Even when restricting to the subset of genes with corresponding models at both gene and isoform levels, 60-75% of isoTWAS models showed better prediction of total gene expression than the corresponding TWAS model (Figure 2c, Supplemental Figure S6-7, Supplemental Table S3-4). Similarly, across all GTEx tissues, isoTWAS shows considerable improvements in prediction, even when considering tissues not derived from the brain (Supplemental Figure S6-7. Supplemental Table S3-4). In a completely independent eQTL reference panel derived from adult frontal cortex brain samples from the PsychENCODE Consortium (N = 850), multivariate isoTWAS models for adult frontal cortex trained in GTEx also outperformed univariate TWAS models in predicting total gene expression (median percent increase in adjusted R² of 21.9%; Figure 2d, Supplemental Table S5-6).

In sum, 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 will yield substantially improved discovery in trait mapping. Altogether, for 48 tissues in GTEx, we built significant predictive models for 26,000 to 33,000 isoforms across 5,000 to 6,000 unique genes per tissue (Supplemental Figure S9). Our results underscore that multivariate methods that incorporate the correlation between isoforms provide utility in prediction of both total gene and individual isoform expression, with strong implications on increased discovery in trait mapping⁴⁵.

isoTWAS increases discovery of trait associations across 5 neuropsychiatric disorders

To directly test our central hypothesis that isoform-centric multivariate predictive modeling will yield substantial power gains in complex trait mapping, particularly in brain tissues for relevant neuropsychiatric traits, we next directly compared isoTWAS and TWAS discovery across 5 neuropsychiatric GWAS results using models trained in larger functional reference panels derived from adult²⁶ and fetal²³ cortex from the PsychENCODE Consortium (Methods, Supplemental Figure S10). Our results demonstrate that isoTWAS increases discovery of trait associations without clear indication of elevated test statistic inflation, especially at genetic loci that house trait-associated genetic variants detected by GWAS (Figure 3).

• Download figure
• Open in new tab

Figure 3: 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) Mean percent increase in approximate χ²-test statistic (squared Z-score), which is proportional to increase in effective sample size, for isoTWAS trait associations over TWAS trait associations. (e) Number of gene-trait associations overlapping GWAS risk SNPs within 0.5 Mb using splicing-event-based TWAS (purple), isoTWAS (blue), or either (orange) at conservative permutation-based significance thresholds.

Using PsychENCODE data, we trained 7,530/3,402 models using isoTWAS/TWAS in adult cortex and 8,331/2,112 isoTWAS/TWAS models for fetal cortex, passing the heritability and cross-validation R² cutoffs mentioned previously (Methods; Supplemental Tables S5-6). We next integrated these models with 5 neuropsychiatric GWAS (Methods, Figure 3a): SCZ⁴⁶, BD⁴⁷, ALZ⁴⁸, ASD⁴⁹, and ADHD⁵⁰. At a conservative permutation-based significance thresholds, we detected far more trait associations with isoform-level mapping than gene-level (Supplemental Figure S11, Supplemental Data 1-4). Specifically, across all 5 neuropsychiatric GWAS and both adult and fetal tissue, isoTWAS identified significant associations within 92 unique GWAS loci compared with 27 detected by TWAS - nearly a 3.5-fold increase in discovery (Figure 3b; Methods). For SCZ, out of the 284 GWAS loci that have been identified to date⁴⁶, isoTWAS prioritizes genes within 81 unique loci across fetal and adult reference panels, compared with 22 from gene-level TWAS of which 19 were jointly identified. Over 89% (24/27) of gene-level TWAS associations within these GWAS loci were concordantly identified by isoTWAS. The standardized effect sizes across gene- and isoform-level mapping (for associations at nominal P < 0.05) are also positively correlated (r = 0.90, P < 2.2 × 10⁻¹⁶), with only 7.8% (103 out of 1220) of these effect sizes in opposite directions (Figure 3c). These results emphasize that isoTWAS recovers a large majority TWAS associations and significantly increases the GWAS signal that can be potentially attributed to underlying transcriptomic mechanisms. For example, out of the 284 GWAS loci that have been identified for SCZ⁴⁶, we find overlapping isoTWAS-identified genes for 46 and 52 loci in adult and fetal brain, respectively. These isoTWAS associations also recover 16 of 17 and 7 of 9 GWAS-overlapping TWAS associations from adult and fetal brain, respectively.

We next sought to ensure that this increase in trait mapping discovery reflected true biological signal rather than test statistic inflation due to the increased number of tests (approximately 7 to 8-fold increase in number of tests). Quantile-quantile plots of trait association gene-level P-values across isoTWAS and TWAS do not show any apparent differences in test statistic inflation (Supplemental Figure S12). Through empirical Bayesian methods⁵¹, we estimated inflation in gene-level test statistics across TWAS and isoTWAS with no marked differences in the 95% credible intervals for test statistic inflation across the 5 traits. In fact, the posterior median for estimated test statistic inflation for isoTWAS were generally smaller than for TWAS (Supplemental Figure S12, Supplemental Table S7). Using a heuristic to estimate increases in effective sample size (Methods), we observed an approximate increase in effective sample size of 15-35% when using isoTWAS compared to TWAS (Figure 2d, Supplemental Table S7). These analyses indicate that isoTWAS discovery is well calibrated in real data compared with TWAS

Lastly, in our comparison of methods in real data, we compared discovery using isoTWAS to discovery using splicing-event level trait mapping. Briefly, for the fetal brain cohort, 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, we found that isoTWAS discovers prioritizes features at far more independent GWAS loci (61 loci) than this splicing-event level trait mapping (31 loci), with 16 loci jointly identified (Figure 3e). This minimal overlap may suggest that these analyses could be used complementarily with careful interpretation of results, bridging annotation-based and -free quantifications of splicing patterns in the transcriptome. Taken together, isoTWAS’s specific focus on modeling isoforms of a gene provides important gains over considering only total gene expression or intron usage in identifying trait associations for genes and their transcripts.

Simulations support improved power and calibrated null of isoTWAS across genetic architectures

To systematically explore the effects of varied genetic architecture on the performance of isoTWAS and TWAS, we conducted extensive simulations using European-ancestry SNP data from the 1000 Genomes Project⁵⁴ (Methods, Figure 4a, Supplemental Figure S14). In general, simulations reinforce results from real data: isoTWAS increases prediction of total gene expression and power to detect trait associations, especially when genetic effects on expression are different across isoforms of a gene. First, we assessed scenarios in which isoTWAS improves prediction of gene expression over TWAS. Similar to results from real data, the multivariate prediction models (particularly multivariate SuSiE³⁸) provided the largest average adjusted R² for isoform expression prediction across most simulation settings (Supplemental Figure S15, Supplemental Data 5). At sparser isoQTL architectures, the optimal isoTWAS model greatly outperforms the optimal TWAS model, with the 25% and 75% quantiles of adjusted R² increasing by 1.3- and 14-fold (Figure 4b-c, Supplemental Figure S16-17, Supplemental Data 6). Performance gains decreased somewhat with denser isoQTL architectures, although real data reflects this 0.1-1% sparsity (i.e., 1-10 causal e- and isoQTLs per gene or isoform)³³. Similar performance gains were observed from simulations in a second simulated test locus around XRN2 (Methods, Supplemental Figure S18, Supplemental Data 7).

• Download figure
• Open in new tab

Figure 4: IsoTWAS models predict gene expression with more accuracy than TWAS models.

(a) Simulation setup to generate isoform expression with specified isoQTL architecture, controlled expression heritability, number of isoforms, and inter-isoform correlation structure. See Supplemental Figure S3 for more details. 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 (b) 0.10 and (c) 0.25. Here, we sum predicted isoform expression from isoTWAS models to form predicted gene expression and compare to predicted gene expression from a TWAS model.

We next introduced GWAS traits into our simulation framework to benchmark the power and false positive rate (FPR) of isoTWAS (Methods). First, we found that the FPR is controlled at 0.05 for isoform-level mapping using ACAT, similar to that of gene-level mapping using TWAS (Supplemental Figure S19, Supplemental Data 8). For a simulated trait, we modeled causal effect architectures for a genomic locus with two isoforms under three main scenarios (Methods; Figure 5, Supplemental Figure S20): (1) where the true trait effect is from 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 both isoforms are ‘effect’ isoforms, with varying magnitudes and directions of association. This first scenario shows clear increases in power for TWAS over isoTWAS, but this advantage of TWAS over isoTWAS decreases as the causal proportion of QTLs increases (Figure 5a, Supplemental Data 9). As effects on the trait varied for isoforms for the same genes (Figure 5b-c, Supplemental Data 10-11), we observed stark 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 carries a trait effect or when two effect isoforms of the same gene have different directions of effects. Our rigorous simulations reinforce the prediction and power advantages of isoTWAS over TWAS, suggesting varied isoQTL architecture and isoform effects on traits for isoforms of the same gene as key features that may drive advantages we observed in real data.

• Download figure
• Open in new tab

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

(a) Simulation settings for gene- and isoform-level associations with complex trait (Methods). (b) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10⁻⁶, Y-axis) across causal proportion of isoQTLs (X-axis), facetted by sparse and dense correlations between isoforms (right margin) and shared or different isoQTLs across isoforms (top margin). (c) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10⁻⁶, Y-axis) across proportion of gene expression explained by effect isoform (X-axis), facetted by sparse and dense correlations between isoforms and shared or different isoQTLs across isoforms. (d) Power to detect gene-trait association (proportion of tests with P < 2.5 × 10⁻⁶, Y-axis) across ratio of effect sizes for 2 effect isoforms, facetted by sparse and dense correlations between isoforms and shared or different isoQTLs across isoforms. Across all plots, isoform and total gene expression heritability is set to 0.10 with 2 simulated isoforms and causal proportion of 0.01, unless otherwise specified.

isoTWAS identifies biologically-meaningful trait associations undetectable by TWAS

Lastly, we applied probabilistic isoform-level fine-mapping to genetic loci where isoTWAS prioritized multiple isoform-trait associations to estimate sets of isoforms that are most likely to explain the trait association at the genetic locus (Methods). The full number of gene-level and isoform-level trait associations at each step of the testing pipeline are summarized in Figure 6a. In total, forming 90% credible sets explain the association with the trait at the genetic locus, using adult frontal cortex isoform expression, we detected trait associations at 25 genes (26 isoforms) for ADHD, 19 genes (21 isoforms) for ALZ, 17 genes (18 isoforms) for ASD, 28 genes (28 isoforms) for BD, and 136 genes (136 isoforms) for SCZ. Using fetal cortex isoform expression, we detected trait associations at 98 genes (105 isoforms) for ADHD, 20 genes (21 isoforms) for ALZ, 4 genes (4 isoforms) for ASD, 22 genes (22 isoforms) for BD, and 65 genes (66 isoforms) for SCZ (Supplemental Figure S21-S25, Supplemental Tables S8-9). Overall, we find, for isoforms associated with one of the 5 traits at nominal P < 0.05 using both adult and fetal frontal cortex expression, there is only a moderate correlation in isoTWAS effect sizes (r = 0.63, P = 4.4 × 10⁻¹¹, Figure 6b). This level of correlation may suggest that there is a difference in isoform-specific effects in the developing brain, compared to the adult brain, on these neuropsychiatric traits.

• Download figure
• Open in new tab

Figure 6: isoTWAS detects biologically-meaningful trait associations with ADHD, ALZ, ASD, BD, and SCZ using adult and fetal brain tissue.

(a) Number of gene and isoform (vertical margin) associations across various significance thresholds (color) in both adult and fetal brain (horizontal margin). (b) Scatter plot of isoTWAS effect sizes across adult (Y-axis) and fetal (X-axis) brain for isoform-trait associations. Black diagonal line shows the 45-degree line and the blue line shows the regression line. (c) Forest plot of isoTWAS effect sizes and 95% confidence intervals (X-axis) for isoforms that are within 10 kilobases of a significant GWAS SNP at P < 5 × 10⁻⁸. (d) Lollipop plot of enrichment ratio of various gene ontologies overrepresented by isoTWAS-identified genes at Bonferroni-adjusted P < 0.10.

Across all 5 traits, isoTWAS detected 15 trait-associated genes using fetal frontal cortex models and 42 genes using adult cortex models with gnomAD-calculated probability of loss-of-function intolerance (pLI) score of ≥ 0.8 (Supplemental Tables S8-9), suggesting that these genes are particularly intolerant to protein-truncating variation⁵⁵ and that isoTWAS detects biological-relevant signal due to its isoform-specific focus. Transcripts of many of these genes also house GWAS-significant variants within 10 kilobases, with effect sizes and 95% confidence intervals shown in Figure 6c. For example, we identified an association with SCZ with imputed fetal expression of ENST00000586144, an isoform of CXXC1 (P = 2.56 × 10⁻¹¹ in isoTWAS, pLI = 0.90, chromosomal location of 18q21.1), with a GWAS-significant variant within 10 kilobases upstream of the gene’s transcription start site. CXXC1 belongs to the SET1-COMPASS complex, implicated in transcriptomic pathways underlying schizophrenia risk⁵⁶. Another transcript whose imputed fetal brain expression was associated with SCZ risk is ENST00000511382, an isoform of the gene PDE4D (P = 2.4 × 10⁻⁸, pLI = 1, chromosomal location of 5q11.2-q12.1), with a GWAS SNP within the gene body (rs7701440). A previous GWAS has implicated PDE4D in associations with SCZ risk⁵⁷, but our analysis directly ties the genetically-regulated component of its isoform ENST00000511382 to SCZ risk through probabilistic fine-mapping. Imputed adult brain expression of ENST00000492146 (P = 6.9 × 10⁻¹², 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 (rs2710339). 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, overrepresentation analysis of isoTWAS-priortized genes for these 5 traits showed multiple relevant biological processes, cellular components, and pathways, at Bonferroni-adjusted P < 0.10: cation transportation, immune response ontologies, and plasma membrane and organelle organizations (Figure 6d). These results further emphasize that isoTWAS identifies sets of trait-associated genes with biological significance.

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 a few examples of isoforms prioritized by isoTWAS, all using adult frontal cortex tissues, which show limited evidence of QTLs on the total gene expression level. 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^59,60. At this locus, we found a strong isoQTL signal (P < 10⁻¹⁵) that intersects GWAS associations at P < 5 × 10⁻⁸; however, the eQTL signal failed to even reach P < 10⁻⁴ (Figure 7a). The lead isoQTL rs2125232 shows a strong positive effect on ENST00000492957 but an insignificant positive association with AKT3, as the number of alternative alleles at the SNP increased. This isoQTL effect was in the opposite direction as the negative associations with SCZ from GWAS.

• Download figure
• Open in new tab

Figure 7: isoTWAS implicates isoforms of AKT3, GIGYF2, and KMT5A in genetic associations with schizophrenia risk.

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. (bottom) Boxplots of gene and isoform expression by genotype of the lead isoQTL SNP and forest plot of the lead isoQTL’s effect and adjusted 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, GIGYF2 gene expression, ENST00000373566 isoform expression. (c) SCZ risk, KMT5A gene expression, ENST00000537270 isoform expression.

Similarly, we find a strong isoQTL signal for ENST00000373566 but a weak eQTL signal of GIGYF2 in the 2q37.1 locus (pLI = 1), in another association with SCZ (Figure 7b). GIGYF2 is a translation inhibitor, which, with the protein 4EHP, decreases translation of defective messenger RNAs to assist ribosome-associated quality control; decreased ribosomal quality is a hypothesized mechanism explaining neuropsychiatric disorders, as observed in both mice and humans with decreased GIGYF2 function⁶¹. Interestingly, the lead SNP rs2675972 has a positive effect on both isoform and a negative effect on SCZ. Lastly, our analysis implicated ENST00000537270, an isoform of KMT5A (pLI = 0.99, found in 12q24.31), in an association with SCZ risk (Figure 7c). There is a strong isoQTL signal at the locus with no significant eQTL signal; the lead isoQTL rs28577594 has a strong negative effect of the transcript-isoform expression and a protective effect on 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^46,63.

Lastly, our isoTWAS analysis detected an association between ALZ and ENST00000252491, an isoform of APOC1 in the 19q13.32 locus, a gene implicated in ALZ risk^64,65, with APOE and ACE, associated genes in the locus (Supplemental Figure S26). Interestingly, when we condition on APOE risk variants strongly associated with ALZ risk (rs7412, rs429358)^48,66 using analyses highlighted in Gusev et al^9,67, we find that the association between ENST00000252491 and ALZ persists (without adjustment for APOE variants: P = 1.9 × 10⁻⁷; with adjustment for APOE variants: P = 2.1 × 10⁻⁶). Despite limitations in these conditional analyses using GWAS summary statistics⁶⁸, our results may suggest an independent signal in this locus with an underlying transcriptomic mechanism through ENST00000252491.

These SNP-by-SNP views into the GWAS, eQTL, and isoQTL signal at a locus demonstrate why isoTWAS gives a clearer picture of the mechanism at the locus over TWAS. Integration with total gene expression cannot detect an overlapping genetic, transcriptomic, and phenotypic signal, especially at these sample sizes. When we specifically model the isoform transcripts belonging to the same gene family, we detect the more nuanced isoform-specific signals that may underlie a genetic association.

Training predictive models of isoform expression

Step-wise association 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, GReX 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,68. 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 potentially correlation in test statistics for isoforms of the same gene.

Accordingly, we adopt methods from the stageR framework (Supplemental Figure S2)⁷⁸. In the first step, for every isoform with a significant model, we generate a TWAS test statistic using either linear regression with direct GReX imputation 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 Bonferroni procedure, but the Benjamini-Hochberg procedure can also be applied for less conservative false discovery control. In the second step, for isoforms for genes with FDR-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 can use the null distribution to generate a permutation P-value for the original test statistic for each isoform. Additionally, we can 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. In addition, we employ gene-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.

GTEx processing and model training

We quantified GTEx v8³³ RNA-Seq samples for all 49 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³⁵. We then extracted reads from the available genome-aligned GTEx BAM files into FASTQ files and used the FASTQ files as input to Salmon with mapping validation and corrections for sequencing and GC bias. We estimated 50 inferential bootstraps for isoform expression using Salmon’s Expectation-Maximization algorithm. We then imported Salmon isoform-level quantifications and aggregate 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^80,81. 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.

SNP genotype calls were derived from Whole Genome Sequencing data for samples from individuals of European ancestry, filtering out SNPs that have less than 5% minor allele frequency and deviate from Hardy-Weinberg equilibrium at P < 10⁻⁵. We further filtered out SNPs with less than 1% minor allele 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. Before training genetics models of gene-level and isoform-level expression, we estimated the heritability of gene and isoform expression using all SNPs within 1 Mb with GCTA-GREML^82,83. We only considered genes and isoforms with positive heritability estimates with nominal P < 0.05. We borrowed techniques from FUSION⁴ to train gene-level models using the most predictive, based on McNemar’s adjusted 5-fold cross-validation (CV) R², of models trained using elastic net regression³⁶, linear mixed modeling⁴⁰, and SuSiE³⁸. 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; we only used multivariate elastic net, multivariate SuSiE, univariate modeling across isoforms, and multivariate LASSO regression with covariance estimation, as these four methods showed far better predictive performance across all 6 methods (Supplemental Figure S2). For gene expression prediction in TWAS, we employed univariate elastic net, linear mixed modeling, and SuSiE, as in FUSION^4,36,38,76, using all SNPs within 1 Mb of the gene body. All isoTWAS models generated are available at https://zenodo.org/record/6795947⁶⁹.

Adult and fetal PsychENCODE processing and model training

We quantified both adult (N = 850)^26,27 and fetal frontal cortex²³ (N = 205) RNA-Seq samples using Salmon v1.5.2³⁰ in mapping-based mode. We used the same Salmon indexed transcriptome used in the GTEx analysis. We then extracted reads from the available genome-aligned BAM files into FASTQ files and used the FASTQ files as input to Salmon with mapping validation and corrections for sequencing and GC bias. We estimated 50 inferential bootstraps for isoform expression using Salmon’s Expectation-Maximization algorithm. We then imported Salmon isoform-level quantifications and aggregate 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^80,81. 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 for fetal and adult brain, respectively) used in the original QTL analyses in Gandal et al and Walker et al.

Typed SNPs with non-zero alternative alleles, MAF > 1%, genotyping rate > 95%, HWE P < 10⁻⁶ were first imputed to TOPMed Freeze 5 using minimac4 and eagle v2.4^84,85. genotype calls were derived from Whole Genome We then filtered to SNPs that had imputation accuracy R² > 0.8, are biallelic, and have rsIDs. Furthermore, we filtered out SNPS with MAF < 0.05 and deviate from Hardy-Weinberg equilibrium at P < 10⁻⁶.

Details of the model training pipeline are summarized in Supplemental Figure S10 and are equivalent to those used to train models in GTEx data. All isoTWAS models generated are available at https://zenodo.org/record/6795947⁶⁹.

Gene- and isoform-level trait mapping

We conducted gene- and isoform-level trait mapping for five neuropsychiatric traits: (1) schizophrenia (SCZ; N_cases = 69,369; N_controls = 236,642)⁴⁶, (2) bipolar disorder (BD; N_cases = 41,917; N_controls = 371,549)⁴⁷, (3) attention deficit hyperactivity disorder (ADHD; N_cases = 20,183; N_controls = 35,191)⁵⁰, (4) autism spectrum disorder (ASD; N_cases = 18,381; N_controls = 27,969)⁴⁹, and (5) Alzheimer’s disease (ALZ; N_cases = 90,338; N_controls = 1,036,225)⁴⁸. 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. 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⁷⁹.

We estimate 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.

Simulation framework

We adopt techniques from Mancuso et al’s twas_sim protocol^79,87 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 CACNA1E (1,107 total SNPs) using an LD reference panel from European samples from 1000 Genomes Project⁵⁴. We use this ‘gene for these simulations as it displays interesting patterns of alternative isoform expression⁸⁸; however, we note that the selection of this particular gene is ultimately arbitrary. We extended our simulations to the locus around XRN2, as well, with similar results.

Next, we generate a matrix of SNP-isoform effects across different causal SNP proportions p_c, numbers of isoforms t, and “shared” or “different” SNP-isoform architectures. We define “shared SNP-isoform architecture as a setting where a SNP either has non-zero effects on all isoform or none at all, such as a case where a variant lies within a gene’s promoter and acts as a broad gene-level cis-eQTL. We define “different” SNP-isoform architecture as one in which variant(s) exhibit distinct (zero or non-zero) effects across different isoforms, such as the case where a variant is proximal to an alternative transcription start site shared only by a subset of isoforms. Lastly, we compare “sparse” (all off-diagonal elements are near 0) or “dense” (most or all off-diagonal elements are far from 0 and generally greater than 0.3-0.4) correlation structures for the isoform expression levels within a gene locus, using Joe’s C-vine method^89,90. Using matched moments followed by scaling to ensure a set isoform expression heritability of , we can generate the random error matrix that ensures the desired correlation between isoforms. We generate 10,000 simulations for each configuration of the simulation parameters, varying n ∈ {200, 500}, p_c ∈ {0.001,0.01,0.10}, , “shared” and different SNP-isoform effects, and “sparse” and “dense” correlation between isoforms. Full mathematical details are provided in Supplemental Methods and summarized in Figure 4a and Supplemental Figure S14.

We also generate traits under three distinct scenarios, each with 2 isoforms for a given gene:

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 GReX to 0 to generate 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.