The protein sequence of the APOE gene was retrieved from Ensembl genome browser v107 in FASTA format (https://useast.ensembl.org/index.html). Python package Biopython was used to load and process the retrieved sequence. Only the reference isoform (APOE3) and the precursor APOE (pre-APOE) sequences were used in this study. The difference between pre-APOE and mature APOE was the addition of an 18-residue signal peptide at the beginning of the sequence. As a result, previously reported variants with respect to mature APOE, such as C112R and R158C, were reported as C130R and R176C, respectively, in this study.

The C130R variant was manually introduced to create a separate sequence representing APOE4, and R176C was manually introduced to create a separate sequence representing APOE2.

ESM-1b model was retrieved from GitHub (https://github.com/facebookresearch/esm) using PyTorch Hub. The same tokenizer as the original ESM model was used to encode input protein sequences. The variant effect for each amino acid variant (ESM score) was calculated as the log-likelihood ratio between the variant and the corresponding reference amino acid. To show a positive score, we multiplied each prediction score by −1. E S M   s c o r e = − log P V a r i a n t P R e f e r e n c e

The variant was predicted to be more damaging if it had a higher ESM score.

AlphaFold v2 model was run locally using a third-party implementation, namely, LocalColabFold (https://github.com/YoshitakaMo/localcolabfold) (Jumper et al., 2021; Mirdita et al., 2022). The algorithm first implements MMseqs2 (Steinegger and Söding, 2017) to retrieve MSA for the target protein. Then, it predicts the 3D protein conformation for the given sequence.

Due to the high computational cost of running AlphaFold, it was extremely time-consuming to run predictions (in silico mutagenesis) for all possible SAVs in APOE, which would require running AlphaFold 6,023 times (Supplementary Table S1). As a workaround, we retrieved all SAVs in APOE reported in ClinVar (Landrum et al., 2015). First, ClinVar database version 20220507 was downloaded from https://ftp.ncbi.nlm.nih.gov/pub/clinvar/. Second, only variants annotated as inside the APOE gene were kept (n = 69). Third, only non-synonymous single nucleotide variants were kept, and all insertions and deletions were excluded (n = 38). As a result, a total of 38 SAVs were retrieved, and a separate protein sequence was created for each SAV. The predicted 3D protein structure for the wild-type and each mutant sequence was compared using the root-mean-square deviation (RMSD) of atomic positions, which was commonly used as a distance measurement between two protein structures. A variant with a higher RMSD score was expected to have a greater impact on the protein structure. Therefore, the RMSD score was used as a surrogate for AlphaFold’s prediction of the variant’s impact.

Besides the two computational tools described previously, we used two additional web services to measure/predict the stability of the protein with and without the variants. First, the Missense3D database for APOE was retrieved from http://missense3d.bc.ic.ac.uk:8080 (Khanna et al., 2021), which contains 307 pre-calculated predictions in APOE. Second, DynaMut2 was used to predict user supply variants (Rodrigues et al., 2020). The same SAVs retrieved from ClinVar were used and submitted to the DynaMut2 web service at: https://biosig.lab.uq.edu.au/dynamut2/.

To evaluate the performance of the main predictor (ESM-1b model), we retrieved population allele frequencies from gnomAD (https://gnomad.broadinstitute.org/news/2020-10-gnomad-v3-1/). Maximum population frequencies were retrieved for the same SAVs retrieved from ClinVar as described previously.

An Evolutionary conservation score, GERP++ (Davydov et al., 2010), was retrieved from the dbNSFP v4.3a database (Liu et al., 2011; Liu et al., 2020), available at https://sites.google.com/site/jpopgen/dbNSFP.

Additionally, we have retrieved 3 popular tools for predicting protein stability change upon mutation, namely, FoldX (Schymkowitz et al., 2005), DDGun (Montanucci et al., 2019) and Maestro (Laimer et al., 2015). First, FoldX was downloaded from https://foldxsuite.crg.eu/ using the academic license. The “Stability” command was used to calculate the Gibbs energy of protein folding for all 38 potential SAVs. The difference in folding energy between wild-type and mutant sequences was calculated and their absolute values were used to represent each SAV’s impact predicted by FoldX, since both stabilizing and destabilizing mutations may all have substantial impacts on the function of the protein. Second, the DDGun web service, available at: https://folding.biofold.org/ddgun/index.html, was used to make predictions on protein stability change given a list of mutations. Specifically, the wild-type sequence of APOE with a list of IDs for all 38 SAVs was uploaded. A global Delta Delta G (DDG) value was predicted for each of the SAVs, and its absolute value was used to represent each SAV’s impact predicted by DDGun. Third, Maestro v1.2.35 Linux executable file was downloaded from https://pbwww.services.came.sbg.ac.at/?page_id=477. All 38 SAVs were submitted as input for the Maestro program with the wild-type 3D structure of APOE obtained using AlphaFold2. Similarly, the DDG values were obtained from the prediction and their absolute values were used to represent each SAV’s impact predicted by Maestro.

To evaluate the correlation between allele frequencies of the variants and predictions made by computational tools, Pearson’s correlation coefficient was calculated using the Python library SciPy (https://scipy.org/). We calculated the area under the receiver operating characteristic curve (auROC) and average precision scores to evaluate each predictor’s ability in prioritizing potential clinically relevant variants. Specifically, an auROC was calculated by measuring the predictor’s true positive rate (TPR) and false positive rate (FPR) using different score cutoffs. Similarly, the average precision (AP) score was calculated by measuring the predictor’s precision and recall (same as TPR) using different score cutoffs. The formulas for calculating TPR, FPR, and precision are: T P R   R e c a l l = T P T P + F N F P R = F P T N + F P P r e c i s i o n = T P T P + F P

Where TP refers to the number of true positives (correctly predicted ClinVar pathogenic variants), FN refers to the number of false negatives (incorrectly predicted ClinVar pathogenic variants as benign), FP refers to the number of false positives (incorrectly predicted ClinVar benign variants as pathogenic), and TN refers to the number of true negatives (correctly predicted ClinVar benign variants). Both auROC and AP scores were calculated using the Python library sklearn with functions roc_auc_score and average_precision_score, respectively.

Additionally, PyTorch was used to calculate ESM model predictions, and Tensorflow was used to calculate AlphaFold model predictions.

As illustrated in Figure 1, the entire length of the APOE protein was predicted by the ESM-1b model, and all potential amino acid variants were evaluated as the log odds ratio between the mutant and wild-type predictions. Variants with lighter colors indicate a low predicted likelihood of the existence of a variant at this position, which implies their functional importance. Key functional domains, including a signal peptide, receptor binding domain and lipid binding domain showed higher importance, as illustrated by light color bands. Interestingly, these regions of high importance showed higher conservation scores (GERP++ score), as illustrated by the top panel. In contrast, amino acids from positions 18–45 showed both low conservation and low predicted functional importance. This observed concordance of the ESM prediction with annotated functional domains and evolutionary conservation demonstrated the model’s ability to capture important regions in the APOE gene, given that the gene is only moderately conserved and is quite tolerant to missense variants (Lek et al., 2016).

Additionally, multiple clustering patterns were observed in the prediction heatmap, as illustrated by regions with high predicted values. One of these regions was amino acids 1–18, representing the signal peptide region. While few studies have tried to evaluate the functional importance of variants residing in this region, it is clear that multiple variants can be extremely harmful to the protein’s function.

To illustrate if the scores predicted by ESM-1b can truly reflect function importance at the variant level, we next evaluated allele frequencies observed in a large-scale population cohort, namely, gnomAD, and see if the model’s predictions show correlations with allele frequencies (AFs) of the variants in general populations. Due to purifying selection, variants with lower AFs are more likely to be deleterious, whereas variants with higher AFs are more likely to be tolerated (benign). As illustrated in Figure 2, predictions by ESM-1 showed statistically significant positive correlations with −log₁₀(AFs) in the general population, which indicates its capability of identifying truly functional variants that have undergone purifying selection.

Next, we investigated the predictions made by AlphaFold. AlphaFold gives a per-residue confidence metric called pLDDT (predicted Local Distance Difference Test) score for all alpha-carbon atoms (Mariani et al., 2013). Regions with high pLDDT scores usually have fewer clashes and structural violations. As shown in Figure 3, pLDDT scores correlate with conservation scores (Spearman correlation coefficient = 0.43, p-value = 1.94 × 10⁻¹⁵), which was expected, as AlphaFold prediction relies on MSA data as input, which primarily utilizes conservation data. Additionally, in regions with high pLDDT scores (pLDDT >70), for example, the amino acid’s approximate position from 45–170, only pathogenic variants and no benign variants were reported. Their potential structural importance could explain this observed pattern.

We compared the correlation of the predictions made by four popular computational frameworks, namely, ESM-1b, AlphaFold, Missense3D, and DynaMut2, which measure protein properties from different perspectives (Figure 4). Specifically, ESM model studies comprehensive protein sequence features from millions of protein sequences using a language model. AlphaFold model studies protein sequences and tries to predict 3D protein structures using sequences and available templates. Missense3D adopts a bioinformatics pipeline and evaluates a wide range of structural impacts of an SAV. DynaMut2 model predicts protein stability by learning a series of biochemical and biophysical features from the target proteins. We have examined additional popular in silico tools that can predict protein stability (Caldararu et al., 2020; Pan et al., 2022), including FoldX (Schymkowitz et al., 2005), DDGun (Montanucci et al., 2019) and Maestro (Laimer et al., 2015), but all of them showed inferior performance compared to DynaMut2 in APOE (Supplementary Figures S1, S2; Supplementary Table S2). Therefore, we chose only DynaMut2 as the representative tool for protein stability prediction. Interestingly, benchmarked using ClinVar labels, the DynaMut2, as the best individual predictor among protein stability methods, outperformed predictors from other categories, including ESM and AlphaFold. Therefore, we provided predictions from DynaMut2 for all possible SAVs in APOE in Supplementary Table S3.

All these four selected tools showed no or very weak correlations with each other, which can be both concerning and useful. On one hand, if the inconsistency arises from methodological flaws, their ability to capture useful information is extremely limited. Users should be cautious when adopting these tools in their workflow. On the other hand, if this inconsistency arises from the differences in methodological preferences and their ability to capture different aspects of protein functions, then these tools can provide valuable orthogonal information.

To evaluate the usefulness of the previously described tools and illustrate if their low correlation can be beneficial to explaining variant effects, we obtained top candidates from multiple predictions and examined their biological relevance as a means of validation. The top candidates were obtained based on the predictions made by each of the four tools. We consider variants to be potentially pathogenic if predictions from two or more tools showed indicative of a disruptive effect. Using this ensemble approach (majority vote), five candidate variants were obtained. As shown in Table 1, variants C130R, R163C, and R132C are most likely to be functional. Importantly, DynaMut2, which predicts the stability of the mutant protein sequence, showed destabilizing effects for all these 3 variants. AlphaFold models also predict the top 2 variants to disrupt the key functional domains. Interestingly, the most promising variant, C130R, is the variant that separates the transcript that carries the variant gene (APOE4) from the wild-type transcript (APOE3). The variant replaces Cysteine with Arginine, which was predicted to change a residue state from buried to exposed (Figure 5). The functional importance of the C130R variant was validated by previous studies, which reported the variant to be associated with an elevated risk of AD (Husain et al., 2021; Martens et al., 2022). This observation highlighted the ability to combine multiple functional prediction tools in finding key functional variants.

Next, using a similar approach, we identify one variant to be potentially benign. As shown in Table 2, all four tools predicted the variant to be non-functional. ClinVar reported a conflicting interpretation of pathogenicity for this SAV, meaning that multiple clinical laboratories reported contradictory interpretations for the same variant. Specifically, some studies reported it to be benign while others report it to be uncertain significance, according to the 2015 ACMG-AMP guidelines (Richards et al., 2015). Given its previous uncertain annotations and the fact that all orthogonal in silico methods showed concordant prediction, its function is worth investigating in future studies to confirm whether the SAV is truly benign. All calculated scores for all 38 SAVs analyzed in the study are provided in Supplementary Table S2.