A Premalignant Cell-Based Model for Functionalization and Classification of PTEN Variants

A large number of sequence variations exist in the population. For example, in the first Canadian personal genome project, 207 million sequence variants were identified in the whole genome sequences of just 56 participants (1). Although databases such as ClinVar (2) have been established to curate sequence variations, only 19 out of a total of 1,606 disease-associated variants identified in this study have been identified as clinically relevant. This is the case because, apart from nonsense mutations, it is often difficult to predict whether the identified mutations alter protein stability or function, or have no effect. Collectively, such variants, which are in the overwhelming majority, are known as “variants of uncertain significance” (VUS). The large number of VUS located within disease-specific genes creates a huge information burden for clinicians attempting to craft personalized medical treatments. This is clearly the case in most, if not all, cancer syndromes (3). Therefore, there is a desperate need for methods to systematically determine the pathogenicity of VUS in humans.

A number of computational models have been developed to predict the functional impacts of variants. Nonmachine learning models such as SIFT (4) and Site Directed Mutator (SDM; ref. 5) focus on changes in amino acid properties such as charge, size, or polarity between the native and mutant residues to make predictions. Machine learning models generally utilize evolutionary information (e.g., multiple sequence alignments) and features pertaining to protein structure and function for training and prediction. These include SNAP2 (6), CADD (7), and PolyPhen-2 (8). Computational approaches can quickly process a large number of variants; however, assessment of their prediction accuracy is difficult, as it is limited by the benchmarking dataset. Evidently, their accuracy in clinical test cases is often lower than initially published (9, 10). The fact that different computational models often generate conflicting results further confounds their utility.

Here, we focused on developing a method to experimentally functionalize variants of PTEN at scale. PTEN is a frequently mutated tumor-suppressor gene whose mutations can arise sporadically in a significant proportion of breast, prostate, and endometrial carcinomas and glioblastomas (11, 12). In addition, germline PTEN mutations drive the formation of hyperplastic cellular overgrowths known as hamartomas. Affected patients are at high risk of developing malignant tumors in numerous tissues including the thyroid, breast, endometrium, kidney, and colon (13). Therefore, an ability to identify all pathogenic PTEN variants would be invaluable for clinicians working to implement timely treatment and appropriate clinical management of the affected patients.

A major challenge in functionalizing PTEN variants is that the molecule is a dual-specificity phosphatase. As a lipid phosphatase, PTEN hydrolyzes PI(3,4,5)P3 into PI(4,5)P2, thereby antagonizing PI3K activity in the Akt signaling pathway. Misregulated Akt pathway facilitates the emergence of proliferative, invasive, and survival phenotypes that can contribute to tumor formation and progression (14). In addition, PTEN also shares significant homology with the protein tyrosine phosphatase (PTP) family. The functional consequence of PTEN's protein phosphatase function is less clear, although new findings highlight a role in cell adhesion that can affect cell survival (15). One of its protein substrates is itself; PTEN can dephosphorylate its C-terminal tail at T366, which is important for PTEN's ability to suppress cancer cell invasion (16). Point mutants that selectively disrupt the protein and lipid phosphatase functions show they each inhibit proliferation in glioblastoma cells, and both are required to inhibit invasion (17). Thus, both phosphatase activities of PTEN have been implicated in cancer and hence, a model for functionalizing PTEN variants should ideally encompass both activities.

Several methods have been developed to functionalize PTEN variants. One example is a humanized yeast assay that has been used to rapidly probe the functional effects of approximately 86% of all possible single residue PTEN mutations (18). Another assay utilizes HEK cells and is designed to reveal functionally compromised variants due to reduced protein stability, and is generalizable to other genes (19). Both assays are high throughput, but each has its limitations. The yeast assay is limited to testing only the lipid phosphatase activity of PTEN, whereas the HEK assay does not measure the functional consequences of PTEN mutation. Indeed, Matreyek and colleagues estimated that a significant proportion of known pathogenic PTEN variants do not have reduced stability (19).

We sought to develop a cell-based assay that closely captured PTEN's multiple roles in tumor suppression. Experimentally, the mammary epithelium is exquisitely sensitive to PTEN dosage and function (20, 21). We developed a 3D suspension culture–based assay that utilized a premalignant, anchorage-independent phenotype of normal mammary epithelial cells lacking PTEN to characterize the functional significance of variants in PTEN. This assay was quantitative, reproducible, and had a short run time of 3 to 4 days. We also determined the relative stability of variants as part of the workflow. As a proof of concept, we functionalized 47 variants spanning its phosphatase and C2 domains. Our method successfully identified a mutational hotspot in the well-characterized phosphatase domain of PTEN. We built a machine learning model trained using ClinVar data to classify variants based on the functional scores obtained in our assay. Our model achieved 87.5% accuracy at predicting whether variants were benign or pathogenic and outperformed current computational models.

MCF10A parental and PTEN^−/− cells were purchased from Horizon Discovery and verified by Western blotting. Cells were cultured according to published protocols (22) and were maintained in a 37°C incubator with 5% CO₂. Cells were passaged 35 times, at which point, a previously frozen vial was thawed. Mycoplasma was tested monthly by direct DNA staining with DAPI.

GFP-P2A-PTEN was PCR amplified using the following primers and cloned into AgeI/XhoI sites of pmaxGFP plasmid (Lonza): forward 5′-GCGCaccggtATGGTGTCCAAAGGAGAAGA-3′ and reverse 5′-GCGCctcgagCCCACTTTGTACAAGAAAGC-3′. This plasmid became known as the wild-type (WT) plasmid. Subsequent variants were cloned into AscI/XhoI sites of the WT plasmid using the following primers: forward 5′-CCTggcgcgccGCCACCAGCTGCAATGACAGCCATCATCAAAGAGATCGTTAGCAGAAAC-3′ and reverse 5′-CGCGctcgagTCAGACTTTTGTAATTTGTGTATGCTGATCTTCATCAAAAGGTTCATTCT-3′. Transfection was carried out using Lipofectamine 2000 (Thermo Fisher Scientific) according to the manufacturer's protocols.

Images were acquired on the Leica DMI-4000B microscope equipped with the DFC-345FX camera and N-PLAN 10x objective. All images were further processed using the ImageJ software (NIH).

Cells were sorted using a BD Influx Cell Sorter configured with a 100 μm nozzle at 18 psi at a flow rate <5,000 evt/s. GFP was excited with a 488 nm laser, and the resulting emission was detected using a 530/40 bandpass filter.

Cells were seeded in 6-well ultralow attachment plates (Corning Costar, #3471) at 3 × 10⁴ cells per well. Spheroids were allowed to form in 24 to 36 hours. Spheroid formation was quantified using a 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (M2128, Sigma). MTT was added to each well at 1 mg/mL for 2 hours after cells were pelleted and lysed with DMSO. A spectrophotometer plate reader was used to measure OD 560 nm.

Cells were lysed using RIPA buffer (50 mmol/L Tris, pH 7.5, 150 mmol/L NaCl, 1% Nonidet P-40, and 0.5% sodium deoxycholate) containing protease inhibitors (Roche cOmplete ULTRA tablets) before separation by SDS-PAGE. Western blotting was performed according to the manufacturer's guidelines (BioRad). PTEN primary antibody (138G6 Rabbit mAb #9559, Cell Signaling Technology) was used at 1:2,000. GFP primary antibody (R&D Systems, AF4240) was used at 1:1,000. Goat anti-mouse and goat anti-rabbit secondary antibodies (BioRad, 1721011 and 1706515, respectively) were used at 1:2,000. The chemiluminescent substrate was purchased from Thermo Scientific (34580). Densitometry was performed using ImageJ.

We quantified the PTEN and GFP Western blots of each variant using imageJ. The PTEN:GFP ratios of each variant is normalized to those of wild-type, and each variant is done in at least three biological replicates.

The method for measuring PTEN variant stability in HEK cells is essentially as previously published (23). In brief, HEK 293T TetBxb1BFP Clone 4 cells were transfected with the pCAG–NLS–HA–Bxb1 expression vector and attB-EGFP-PTEN_IRES-mCherry-2A-BlastR plasmids encoding WT or variant PTEN proteins. Cells were switched to Dox-containing media 2 days after transfection and analyzed by flow cytometry at least a week after transfection. Cells were fixed in 4% formaldehyde, and flow cytometry was performed with a BD LSRII flow cytometer configured to analyze BFP, eGFP, and mCherry. Analysis was performed with the FlowJo software. Before analysis of fluorescence, live, single cells were gated using FSC-A and SSC-A (for live cells) and FSC-A and FSC-H (for single cells). We next gated for recombined BFP⁻/mCherry⁺ cells. We created a new ratiometric parameter dividing GFP by mCherry fluorescence and took the geometric mean of this population to ascribe each variant a single score. This process was repeated for four independent replicates.

All bioinformatic analyses were carried out using custom Python (3.6) scripts with Pandas (0.24.2) and numpy (1.16.3) packages. Machine learning, cross-validation, and model evaluation were done using Sci-Kit Learn (0.20.3). All plots were generated using Matplotlib (3.0.3).

When short-listing variants, we focused on sampling different locations of the PTEN protein for coverage, with an emphasis on the phosphatase domain. Thus, we tested variants regardless of their origin, because PTEN loss of function (LOF) is important in both hereditary and somatic cancers. Many of our tested variants are both germline and somatic, whereas some are unique to each category (Supplementary Table S1). When building the support vector machine (SVM) model, variants were labeled as “pathogenic” or “likely pathogenic” using ClinVar (germline) classifications whenever available. When a variant had the designation “conflicting interpretations of pathogenicity,” we reviewed all of its submissions on ClinVar and leveraged somatic data from The Catalogue of Somatic Mutations In Cancer (COSMIC) entries to make a final decision on its label. Variants designated “uncertain significance” were not labeled. Population controls were curated using the following criteria: (i) they were represented in the gnomAD noncancer or control datasets; (2) they have not been implicated in any COSMIC tumor. The last date of access to the above databases was September 20, 2019.

We used the SVM classifier in a multiclass strategy. The decision function shape is set to one-vs-one. The reference dataset is divided into training and testing data in a 70–30 split. We also performed 5-fold cross-validation and reported the average accuracy score. The ROC curves were also generated using the SVM classifier with the linear kernel and plotted using Matplotlib. All codes are available on GitHub at https://github.com/jessecanada/Chao-PTEN-functionalization-2019.

Tumorigenesis is a multifaceted and dynamic process that is difficult to recapitulate using a traditional 2D/monolayer cell culture, particularly with respect to the acquisition of true anchorage independence. Thus, we searched for a 3D culture system to better assess a universally important aspect of the transformed phenotype. The MCF10A cell line is unique as it is a spontaneously immortalized human breast epithelial cell line that is nontumorigenic. Thus, these cells are entirely anchorage dependent, and they undergo anoikic apoptosis in suspensions due to the loss of adhesive cues. Conversely, PTEN^−/– MCF10A cells survive in suspension, and form 3D spheroids that arise due to the heterotypic adhesive aggregation with neighboring cells (24). Therefore, we reasoned that the ability to reverse the survival of PTEN null cells in suspension could form the basis of a reconstitution-based assay to test the function of PTEN variants.

As proof of principle, we focused on testing the function of missense variants that contained amino acid substitutions within the coding region of PTEN. Missense variants prove to be difficult to predict computationally and account for over half of known clinical variants in PTEN (25). We curated a testing library of variants suitable for comprehensive evaluation of our assay (Fig. 1A). We strove to generate a library that included variants already presumed to be pathogenic or likely pathogenic as well as variants that have not been reviewed or were found to be of uncertain or conflicting significance. We surveyed the following public variant databases: (i) The Genome Aggregation Database (gnomAD) has a harmonized collection of large-scale exome and genome sequencing data for germline mutations; (ii) COSMIC database comprehensibly contains clinically submitted and expert-curated somatic variants and includes information from The Cancer Genome Atlas and the International Cancer Genome Consortium; and (3) ClinVar curates clinically obtained germline variants and their impacts in human diseases. Using these databases, we selected 47 variants to be tested (Supplementary Table S1).

The variants in our testing library were sampled from the entire length of PTEN, covering both its phosphatase and C2 domains, as well as the linker regions (Fig. 3D). Among these variants, 18 had not been reviewed, 8 had uncertain significance, and 2 had conflicting significance (Fig. 3D; Supplementary Table S1). It is important to note that even though a variant was found in cancerous tissues according to COSMIC, this does not necessarily represent a causative relationship. We used the following criteria to select non–disease-associated variants: variants that were not known to be associated with any tumor in COSMIC and not known to be pathogenic in ClinVar were considered controls, meaning they were not likely to be cancer causing. Because PTEN is a dual-specificity phosphatase, we deliberately chose variants that affect its lipid and protein phosphatase activities. The G129E variant is lipid phosphatase deficient, but protein phosphatase functional (26). The Y138L variant is the opposite: lipid phosphatase functional, but protein phosphatase deficient (17). The G129R and C124S variants are devoid of all phosphatase activities (27).

To reconstitute PTEN^−/− MCF10A cells with WT PTEN and each of the 47 variants, we engineered a construct in which a CMV promoter drives the expression of a cassette containing GFP-P2A-PTEN (Fig. 1B). The P2A self-cleaving peptide allows for the stoichiometric expression of independently folded GFP and PTEN full-length proteins (Fig. 2B). Thus, we could monitor efficient expression of the cassette and purify PTEN-expressing cells from a heterogeneous pool of transiently transfected cells by gating for GFP using FACS (Fig. 1C). Reconstituted cells were then grown on ultralow attachment plates. The rationale here was that fully functional variants should suppress the anchorage-independent cell adhesion/survival phenotype of PTEN^−/− cells such that they do not form spheroids in a manner equivalent to that observed when WT PTEN is used in the reconstitution. Conversely, nonfunctional pathogenic variants should not rescue the PTEN^−/− spheroid phenotype. In addition to assessing spheroid formation optically, we used the MTT assay as a quantitative indicator of spheroid formation, which is a common approach (Figs. 1C and 2; ref. 28). The entire protocol, from transfection to phenotypic analysis, takes 4 days.

We devised the following plans for quality control. If the transfection efficiency was below 10%, we did not proceed with additional experiments. On average, the transfection efficiency for each variant was approximately 25% to 35%. Also, if the PTEN^−/− control cells did not robustly produce spheroids in a given experiment, we excluded that experiment from further analysis. We observed that FACS could induce cell death, leading to low readings in the MTT assay and false negative results. To circumvent this issue, we calculated cell viability immediately after FACS prior to suspension culture and excluded the experiment if viability was below 20%. We carried out a minimum of three independent biological replicates for each individual variant tested.

Next, we conducted pilot experiments (Fig. 2) to verify the variant functionalization pipeline using PTEN^−/− cells transfected with a control vector or vectors that contained either a functional WT PTEN or a previously characterized LOF, phosphatase-dead C124S variant (27) as well as parental MCF10A cells. After transfection and FACS, cells were maintained in suspension culture for 24 hours and assessed optically and by MTT-based staining (Fig. 2A). As expected, parental MCF10A cells, which are anchorage-dependent (24), remained almost exclusively as singlets that were largely MTT negative and hence not viable. In contrast, PTEN^−/− cells aggregated into robust spheroids that were approximately 100 μm in diameter and viable based on strong MTT staining. Reconstitution of PTEN^−/− cells with WT PTEN appeared to fully restore these cells to the parental MCF10A phenotype, and the cells remained mostly as singlets and were nonviable. Reconstitution of PTEN^−/− cells with the C124S phosphatase-dead variant showed no rescue of the PTEN^−/− phenotype.

We verified the expressions of GFP and PTEN by Western blotting (Fig. 2B). Although the steady-state protein levels were generally similar for both GFP and PTEN, there was a slight increase in the transgenic forms in PTEN^−/− cells compared with the endogenous protein in MCF10A cells.

Finally, we quantified viable cell spheroid formation in suspension by extracting and solubilizing MTT and assessing it photometrically (Fig. 2C). Similar to what we observed optically, we found PTEN^−/− cells that formed robust spheroids showed 3-fold greater MTT activities when compared with parental cells that did not form spheroids. Expressing WT PTEN also lowered the MTT activity while the C124S variant did not (Fig. 2C). Therefore, using this workflow, we could quantitatively evaluate the effects of these variants on PTEN function in a relevant tumor-suppressive context based on loss of anchorage dependence.

We screened 47 missense variants (PBD ID: 1D5R; Fig. 3A) and calculated LOF scores by scaling MTT assay measurements such that fully rescued PTEN^−/− cells with WT PTEN were given a score of 0, and the empty vector had a score of 1 (Fig. 3B). Plotting the LOF scores from low to high revealed our method detected a continuous range of PTEN function with strong linearity (R² = 0.97; Fig. 3B). The 90th and 10th percentile scores were 0.90 and 0.18, respectively, with a median score of 0.54, indicating our method was sensitive enough to detect moderate LOF variants. We plotted LOF scores for variants in order of their amino acid position (Fig. 3C). The C124S variant (red bar) as before showed no activity (LOF = 0.88). Interestingly, the Y138L variant, which is only protein phosphatase–dead (red circles), had an intermediate LOF of 0.55, indicating our assay detected PTEN protein phosphatase activity in addition to lipid phosphatase activity. Control variants (yellow circles) as expected all had low LOF scores, indicating they were likely functional. Variants listed in ClinVar and other cancer-related databases were color-coded in green and blue respectively and showed a range of LOF scores.

We compared the trend in LOF scores with PTEN domains and ClinVar annotations in more detail (Fig. 3D). Even though ClinVar annotations were sparse, we noticed variants that have been labeled “likely pathogenic” and “pathogenic” were also LOF in our assay. Thus, our data validated some functional assumptions that have been made based on clinical correlation (29). Next, we closely examined variants in the catalytic pocket, which contains a conserved P-loop with a motif (HCXXGXXR, residues 123–130) that is commonly found in PTPs and dual-specificity phosphatases (Fig. 3E; ref. 30). Variants devoid of all phosphatase activities were strongly LOF, including the two variants of the catalytic cysteine at position 124 (C124S, C124R) and G129R (27). In addition, G129E, which lost its lipid phosphatase but retained the protein phosphatase activity, also had a high LOF score (Fig. 3C; ref. 26). We found variants in the P-loop tended to have higher LOF scores (>0.5), whereas those in the TI-loop and WPD-loop, which are also conserved among PTPs, had more moderate scores (∼0.2–0.5), suggesting variants in the TI- and WPD-loops had a tendency to be more tolerated (Fig. 3E; refs. 30, 31).

As decreased protein abundance caused by protein destabilization due to amino acid substitution could be a mechanism of pathogenicity, we determined the relative abundance of each variant expressed in MCF10A cells by comparative Western blotting. We calculated ratios of PTEN:GFP from Western blots as a measure of PTEN stability for each variant (Fig. 4A; Supplementary Fig. S1 and Supplementary Table S2). We performed logistic regression on the measurements and identified two classes: “destabilized,” mean = 0.22; and “stable,” mean = 0.89 (Fig. 4B). The six destabilized variants were I101T, G129R, Y155C, R173H, P246L, and T382S (Fig. 4A and C), which also showed a substantial degree of LOF (LOF scores > 0.5; Fig. 4C). This was in contrast to stable variants that showed a full range of variant function, having LOF scores ranging from 0.1 to 1. The R173H variant was an outlier in that it was classified as destabilized, but its LOF score was in the low range (LOF score < 0.3). We postulated R173H would need to retain considerable function to produce such a score. Indeed, several functional studies of R173H indicate it may be WT-like. A yeast assay of PTEN PI3K activity reported this variant to be only mildly LOF (18). R173H did not induce genomic instability in contrast to a nonsense variant (R130*) in an endometrial cancer cell line (32). In addition, R173H correlated strongly with WT PTEN in a gene expression profiling experiment comparing the transcript levels of 978 landmark genes, and did not contribute to in vivo tumor formation in mice (32). However, R173H has been found in glioblastoma and PTEN hamartoma tumor syndrome (PHTS; ClinVar, COSMIC), suggesting it may be functional in in vitro assays, but requires the proper genetic context to give rise to a measurable phenotype in humans.

We now compared our variant stability measurements made in MCF10A cells with those obtained using an established method that we will refer to as the HEK-cell stability assay. In the HEK-cell assay, an integration cassette containing GFP-tagged PTEN variants and an mCherry reporter were recombined directly into the genome of HEK cells (23). Thus, by quantifying the GFP:RFP ratio of individual cells by FACS, a stability score for each variant can be determined. We utilized this method to obtain stability scores for 11 variants (Fig. 4D). We selected variants to cover the major domains of the protein, including the phosphatase and C2 domains, the linker region separating these domains and the C-terminal tail. Variants also represented different classes including population controls (Y180H, P354Q), biochemical controls (C124S), ClinVar (A79T, G129R, G132D, T167N, and Y176C), and other databases (P38H, I101T, and D268E). We found the two sets of stability scores were generally comparable (Pearson's r = 0.589, P = 0.06). Importantly, the HEK-cell assay confirmed that I101T and G129R were destabilized, independently verifying the results from MCF10A cells. One variant, G132D, appeared to be destabilized in the HEK-cell assay but was stable according to our MCF10A data. The LOF score for G132D was moderately high (∼ 0.65; Fig. 3C), suggesting the decreased function of this variant in MCF10A cells was not solely due to destabilization of the protein. The difference in stability in the HEK-cell assay may have arisen as a result of differences in folding of this variant in HEK cells and/or due to the variants having a potentially destabilizing N-terminal GFP tag in that assay. Subsequently, the throughput of the HEK-cell assay was significantly improved by the use of massively parallel sequencing (called VAMP-seq), which measured the stability of 4,112 PTEN variants (19). VAMP-seq confirmed that a third variant P246L was also destabilized, supporting our MCF10A classification (Y155C, R173H, and T382S were not analyzed).

We next compared our empirical stability measurements with those predicted by SDM, a popular in silico modeling suite that predicts the effects of amino acid substitution by calculating the free energy difference of the mutant protein (5). We obtained stability scores (as ΔΔG) for 43 of the 47 variants we tested using SDM; K6E, K6I, P354Q, and T382S were outside of the amino acid range of the structural model (PDB ID: 1D5R). We scaled the ΔΔG scores for the 43 variants (between 0 and 1) and compared them with our normalized stability scores (Fig. 4E). Although the two sets of scores were comparable for a few variants including T167N, Q171E, R173H, N262S, and T232A, there was no overall correlation (Pearson's r = 0.02). For the variants we concluded to be LOF due to destabilization, SDM predicted I101T (ΔΔG = −3.12), Y155C (ΔΔG = −1.31), and R173H (ΔΔG = −1.48) to be destabilized, but failed to identify G129R (ΔΔG = 1.24) and P246L (ΔΔG = 2.33; T382S could not be analyzed). Although overall the correlations were low, SDM predicted three out of possible five destabilized variants.

Unstable proteins are generally cleared from the cytosol by the ubiquitin proteasome system (33). To confirm the destabilized variants observed in our assay were cleared by this pathway, we investigated whether inhibiting the proteasome would increase their abundance. We chose to test the variant I101T because it was identified as destabilized by both our assays and SDM (Fig. 4E). MG132 is a highly potent and cell membrane–permeable inhibitor of the 26S proteasome complex (34) and we treated cells transfected with the I101T variant (and WT PTEN) with increasing concentrations of MG132. MG132 increased the amount of I101T detected in the cell lysate by approximately 1.7-fold relative to GFP but had little/no effect on WT PTEN (Fig. 4F). Thus, these data supported PTEN mutations that lead to decreased protein levels were likely a result of protein destabilization and subsequent degradation.

A yeast-based assay was developed to indirectly measure the lipid phosphatase activity of human PTEN (18). In this assay, overexpression of the human PI3K p110α causes slow growth through depletion of its substrate PI (4, 5) P₂ (via its conversion to PI (3, 4, 5) P₃), which can be reversed by overexpression of human PTEN. This approach was used to measure 7,244 PTEN missense, nonsense, or in-frame deletion mutations, and found 1,789 missense mutations with reduced lipid phosphatase activity (18). Because functional annotations for PTEN variants are not widely available, we compared the functional scores from our assay with those from the Mighell and colleagues dataset (herein referred to as the yeast assay; ref. 18).

We compared the two datasets and found all of the variants we tested except for M35V were empirically measured in the yeast assay (Fig. 5A). We noticed variants in the active site of the phosphatase domain were commonly LOF in both assays, consistent with both assays measuring PTEN's lipid phosphatase activity. However, there was more variability outside of this region with the MCF10A assay appearing more sensitive than the yeast assay toward the damaging effects of these variants. For example, toward the N-terminus, A79T, H93Q, H93R, I101F, and N117S were more damaging according to the MCF10A assay. In the C-terminus, variants I135V to T382S also generally were more deleterious in the MCF10A than the yeast assay. This is likely due to the fact that the yeast assay exclusively measures the lipid phosphatase activity of PTEN, whereas the spheroid phenotype is a broader cell function-based assay that assesses the effects on tumor-suppressive anchorage dependence. In addition, although none of these variants were destabilized in MCF10A cells (Fig. 4A), it is possible they exhibit differential stabilities in the yeast assay.

Next, we performed clustering (Fig. 5B) to observe correlations between the two assays (Fig. 5C). The hierarchical clustering algorithm ordered the variants into six clusters, which we used to color-code the scatter plot comparing the MCF10A scores with the yeast scores (Fig. 5C). Overall, we observed a Pearson's r = 0.48 (P = 0.00078), which indicated a moderate linear correlation between the two datasets. A moderate correlation was perhaps anticipated because the yeast assay is limited to measuring the lipid phosphatase activity of PTEN. The most correlative variants were represented by the teal, blue, and purple clusters (Fig. 5C), which corresponded to variants Y180H-R173H, P38H-Q171E, and C124S-G32D, respectively, in the dendrogram (Fig. 5B). These variants likely all affected the lipid phosphatase activity of PTEN. This is further supported by the presence of the lipid phosphatase–specific mutant G129E (purple cluster). The remaining orange, red, and green clusters consisted of variants that had proportionally higher MCF10 LOF scores than yeast LOF scores. This suggested these variants had deficiencies in addition to decreased lipid phosphatase activity, which were potentially important for PTEN's role as a tumor suppressor. The orange cluster in particular contained variants with LOF scores as high as approximately 0.5 in our assay but as low as 0.05 in the yeast assay (e.g., A79T). The protein phosphatase–specific mutant Y138L was one of these variants (∼0.45 MCF10A vs. ∼0.15 yeast), supporting this cluster contained variants with deficiencies other than lipid phosphatase activity. The green cluster (variants K342N to G127R in the dendrogram) also stood out because the LOF scores for these variants were high in the MCF10A assay but low in the yeast assay, suggesting these variants also had defects beyond lipid phosphatase activity. These analyses strongly suggested our assay monitored multiple facets of PTEN function.

To better understand the differences between the two assays, we investigated whether PTEN function beyond its lipid phosphatase activity was important for its ability to suppress spheroid formation. PTEN via its lipid phosphatase activity antagonizes PI3K and the absence of PTEN leads to an overactive PI3K-Akt pathway, which is a major contributing mechanism to tumorigenesis (35). In our assay, if spheroid formation was entirely due to overactive PI3K-Akt pathway, then inhibiting PI3K in PTEN^−/− cells should completely prevent spheroid formation on par with the phenotype of MCF10A parental cells.

We tested a small-molecule inhibitor of PI3K, LY294002, in preventing spheroid formation in PTEN^−/− cells (36). Although LY294002 suppressed spheroid formation in a dose-dependent manner (Fig. 5D), the maximal suppression at saturating concentration did not match the level of MCF10A control cells (Fig. 5D) or of PTEN^−/− cells with WT PTEN (Fig. 1D), indicating increased PI3K-Akt pathway was not solely responsible for the spheroid phenotype. Consistently, the protein phosphatase–only dead mutant, Y138L, did not fully rescue PTEN^−/− cells (LOF = ∼0.5). Thus, although the lipid phosphatase activity of PTEN plays a major role in its anchorage-dependent tumor-suppressive function, other activities of PTEN are important and were captured by our assay.

We sought to derive clinically relevant predictions for PTEN variants using our empirically determined LOF scores, so we created a machine learning model using the SVM algorithm. To train the SVM model, we used specific criteria to construct the reference dataset (see Materials and Methods; Supplementary Table S3). For classification, we chose to use three classes: “benign,” “likely pathogenic,” and “pathogenic.” We reasoned that having three classes would take advantage of the continuous range of LOF scores generated in our assay and would better align with the current annotation standards for PTEN VUS (37).

Classification of variants using this SVM model achieved an impressive accuracy score of 0.875 (Fig. 6A). This model classified variants with LOF scores less than approximately 0.3 as “benign,” and those greater than approximately 0.8 as “pathogenic.” Variants with LOF scores in between these values were categorized as “likely pathogenic.” Using this model, we classified all tested variants (Fig. 6B), including 37 variants that previously did not have a clear ClinVar interpretation (Fig. 6B). We were able to predict only eight variants were “benign,” whereas the majority of the ClinVar VUS were classified as “likely pathogenic” or “pathogenic.” In only one case did our model predict a variant as “benign” when it was annotated as likely pathogenic in ClinVar (R173H). However, as discussed earlier, our result for R173H is consistent with three other functional studies, which all indicate R173H has near WT activity.

Many computational tools for predicting the functional consequences of variants have been developed. We compared our results with the four commonly used, highly rated computational models: CAAD, SNAP2, SIFT, and PolyPhen-2 (for an overview of these tools, see ref. 25). We compared our LOF scores with computational predictions using common statistical metrics by using the labels in the training dataset to benchmark the performance (Supplementary Table S3). We then plotted receiver curves that allowed us to visualize the diagnostic abilities of different systems by analyzing the relationships between true positive rate (sensitivity) and false positive rate (1- specificity) at various thresholds. Thus, a model with no skill (e.g., making predictions by chance) would fall on the diagonal and have an AUC of 0.5. In a binary classification task (benign vs. pathogenic), our assay outperformed all computational models and the yeast assay with an impressive AUC of 0.97 (Fig. 6C). Because MCF10A, CADD, and SNAP2 scores have a continuous range, we tested if they would perform well in a multiclass classification task in differentiating benign, likely pathogenic, and pathogenic variants (Polyphen-2 and SIFT are binary classifiers and were not included in this comparison). The AUC-ROC indicated MCF10A scores performed the best in a multiclass task with an AUC of 0.833 compared with 0.715 for SNAP2 and 0.645 for CADD (Fig. 6D). In this test, the yeast assay performed comparably to our assay.

We next visualized the correlations between MCF10A scores and each set of computational scores (Fig. 6E). The scatter plots were color-coded according to the classification given to the variant by our SVM model. We first observed that the CADD and SNAP2 scores showed a continuous range and correlated only modestly with MCF10A scores (Pearson's r = 0.4483 and 0.4533, respectively). CADD proposes (but does not recommend) a genome-wide cutoff of 15, above which, a variant is likely to be deleterious (7), although alternative methods such as MSC (38) and GAVIN (39) give a PTEN-specific cutoff of 18.97 and 17.33, respectively. Using a lower cutoff of 18.5, only two variants were considered benign according to CADD. These were also categorized as benign by SVM. The remainder of variants classified as benign by SVM had CADD scores ranging from 22 to 34. GAVIN indicates CADD scores > 29.39 are pathogenic. Only half of variants determined to be pathogenic by SVM had scores > 29, and the majority of variants classified as likely pathogenic by SVM had scores below 29. This suggested CADD poorly distinguished between functional and dysfunctional variants consistent with its low AUC (0.645) in the multiclass classification.

Next, we examined SNAP2 scores that appeared to have outperformed CADD according to their AUC. SNAP2 recommends scores above 50 as strongly deleterious, whereas scores below -50 as no effect/neutrality (6). From the scatter plot, besides 2 variants classified as benign by SVM, all variants with SNAP2 scores above 50 were classified as pathogenic or likely pathogenic by SVM. Three variants with scores below -50 were classified by SVM as likely pathogenic, whereas the remainder of variants classified as benign were captured, but three pathogenic variants were not classified as deleterious by SNAP2.

PolyPhen-2 calculates the Naïve Bayes posterior probability that an amino acid change is damaging, but the program also gives a qualitative output (8). As such, PolyPhen-2 scores > 0.908 are predicted to be probably damaging, between 0.908 and 0.446 as possibly damaging, and <0.446 as benign. Overall, PolyPhen-2 performed well (AUC 0.818 in binary classification) and classified all variants designated pathogenic by SVM as probably damaging. No variants were classified as possibly damaging. Many variants classified as likely pathogenic by SVM were classified as benign by PolyPhen-2.

SIFT scores are the normalized probability of an amino acid change, and a score below 0.05 is predicted to be deleterious, whereas all others are tolerated (40). For ease of comparison, we inverted the computed SIFT scores for plotting such that values above 0.95 are deleterious and values below are tolerated. Accordingly, two variants classified as benign by SVM were deleterious according to SIFT, whereas all other benign variants were tolerated. One pathogenic variant was classified as tolerated, whereas all others were deleterious. Five likely pathogenic variants were classified as tolerated according to SIFT. Overall, these four computational methods did reasonably well at predicting the most damaging variants. However, they were less effective at categorizing variants with intermediate LOF scores designated as likely pathogenic by SVM.

Here, we developed the MCF10A spheroid assay, which tests the ability of PTEN variants to suppress anchorage independence, a premalignant phenotype, in order to assess the contribution of these variants to cancer. This assay captured multiple mechanisms of PTEN inactivation in cancer, including loss of PTEN's lipid and protein phosphatase activities and reduced protein stability (41). Accordingly, the assay offers several improvements over other functional assays for PTEN, which either measure its protein abundance only (19) or test exclusively its lipid phosphatase activity (18). Also, the finding that a PI3K inhibitor did not fully suppress spheroid formation of PTEN^−/− cells (Fig. 5D) suggests our use of a premalignant human breast epithelial cell–based model harnessed phosphatase-independent tumor-suppressing functions of PTEN, such as regulation of p53 and cell motility (42, 43). We propose the assay has clinical utility for the functional classification of both germline and somatic PTEN variants. For germline mutations, our assay can be used to assist with the molecular diagnosis of PHTS. For somatic mutations, given the universal importance of anchorage independence in solid tumor progression, and that somatic PTEN missense mutations are regularly found in gliomas (11, 44), prostate cancers (45), and malignant melanomas (46) among others, our assay should also be well suited for assessing the contribution of PTEN variants to tumor progression.

We extended the utility of our assay by building a machine learning model trained on genotypically annotated sequencing data to classify the pathogenicity of variants according to their LOF scores. This allowed us to make predictions about the disease relevance of variants. The linearity of scores in our assay indicated that PTEN variants comprised not only fully functional and complete LOF proteins, but also a full range of partial LOF proteins. Because of this, we chose to classify variants into three classes instead of just two as is often done in computational models. Having the third class “likely pathogenic” allows for a degree of uncertainty that these variants cause disease. In addition, because these variants show partial LOF, their ultimate phenotypic effects may also be more susceptible to genetic, epigenetic, and tissue-specific contexts, making pathogenicity more difficult to predict. Nevertheless, we believe SVM is a powerful approach for multiple reasons. First, the model allows us to interpret LOF scores in relation to disease rather than simply setting arbitrary thresholds, e.g., functional versus nonfunctional. Second, the decision boundaries for the predictions will adapt to increasing data points as LOF scores for more variants are determined. Third, we can fine-tune the SVM model with improved clinical training data, because as more clinical evidence becomes available for model retraining the SVM decision boundaries will be honed to reflect these improvements, ideally reducing the number of likely pathogenic variants and more confidently distinguishing benign versus pathogenic variants.

According to the ClinGen PTEN–specific guidelines, our assay is closely related to the PS3 criteria in supporting pathogenic classifications. At present, the definition of the PS3 criteria is “to include studies showing phosphatase activity < 50%.” We argue that this is problematic because it has not been defined that a 50% loss in lipid phosphatase activity is sufficient to cause cancer, and it does not take into account nonlipid phosphatase activities of PTEN. We propose the MCF10A assay and SVM model could be used to refine the ClinGen PS3 criteria or be used as an additional criterion for classification of germline variants. If implemented as an additional criterion, it could be used to elevate the role for functional data in classification. This is particularly important for the majority of PTEN variants as these are missense variants and are excluded from the “very strong criteria” in the guideline (PVS1), which applies only to predicted null mutations. Thus, because the MCF10A assay is designed to quantitatively measure multiple aspects of PTEN function in malignant transformation, it should provide enhanced functional information for missense variants for classification by ClinGen (see below).

We noticed three variants, K342N, N117S, and A79T, which are of both germline and somatic origins to be drastically less functional in our assay compared with the yeast assay (Fig. 5A), suggesting they were compromised in functions other than the lipid phosphatase activity. SVM predicted K342N to be pathogenic and N117S and A79T to be likely pathogenic (Fig. 6B). K342N has been implicated in Cowden Syndrome (VUS, ClinVar) and reported in several somatic cancers (COSMIC), and is predicted to be damaging by computational models. Consistent with the yeast assay, when tested in vitro, K342N has WT activity against PtdIns(3,4,5)P3, but shows reduced activity toward soluble Ins(1,3,4,5)P4 (47). Ins(1,3,4,5)P4 can trigger calcium influx through the plasma membrane, which plays a pivotal role in cell migration and tumor metastasis (48, 49), properties that may be detected by our assay. N117S is a germline mutation found in patients with PHTS (VUS, ClinVar), but has also been reported as a somatic mutation in lung carcinoma (COSMIC). N117S is predicted to be likely tolerated by computational models (Supplementary Table S3), which would qualify it as likely benign according to ACMG/ClinGen guidelines. A79T is a germline mutation found in patients with PHTS/Cowden Syndrome (ClinVar) and is also found as a somatic mutation in glioma (50), malignant melanoma (46), and gastric cancer (51), and is predicted to be tolerated by computational models. A79T is currently classified as likely benign by the ClinGen PTEN Expert Panel, primarily because its allele frequency is >0.1% (BS1 criterion). However, there is one pathogenic submission on ClinVar due to the patient's clinical presentation showing an 80% match with Cowden Syndrome. Thus, for K342N, N117S, and A79T, the clinical picture may be more complex given that these variants largely retain their lipid phosphatase activity, but according to our assay are compromised in other cancer-relevant functions.

In addition, we noticed our assay detected stronger deleterious effects than the yeast assay for somatic cancer–only variants, H93Q, A126D, V133I, T232A, and T382S, which have little to moderate loss in lipid phosphatase activity (Fig. 5). This suggests these variants affect the tumor-suppressing functions of PTEN outside of its lipid phosphatase activity, which may be important in somatic cancers. Because these variants are somatic-only, they are not annotated in ClinVar, and computational predictions are mixed, predicting H93Q, A126D, and V133I as deleterious and T232A and T382S as tolerated (Supplementary Table S3). In contrast, by SVM all except V133I (classified as pathogenic) were classified as likely pathogenic, suggesting our assay not only detects increased risk of malignant transformation for germline variants, but can also reveal pathogenicity of somatic variants. A likely reason is that the spheroid model mimics key aspects of somatic cancer progression, including abnormal cell adhesion, the ability to evade apoptosis, and the formation of a hypoxic core (28, 52, 53).

Although the MCF10A assay captures multiple functions of PTEN, incorporating additional features into the assay would help improve the classification of variants, especially in cases where there are conflicting results between it and other assays. Using the powerful 3D spheroid model, it should be possible to study the subcellular targeting of PTEN, activation of the PI3K pathway, and measure PtdIns(3,4,5)P3 levels (and hence PTEN lipid phosphatase activity) by immunofluorescence microscopy for every variant studied. This expanded dataset can be integrated with the SVM model and through capturing these additional features should enhance its predictive power. For example, PTEN nuclear accumulation has been observed in a variety of cancers, presumably due to its decreased availability at the plasma membrane to control PtdIns(3,4,5)P3 levels and PI3K signaling (54). Hence, studying the subcellular localization of PTEN variants could reveal novel mechanisms of PTEN dysfunction in cancer in addition to direct effects on its lipid phosphatase activity.

At the time of writing, over two thirds of PTEN missense variants in ClinVar currently do not have a clear interpretation of significance. It has been shown uncertain or conflicting variant interpretations can cause psychologic distress in patients (55). Thus, the ability to reclassify these VUS will have significant clinical benefits for patient well-being in addition to diagnosis and treatment of PHTS and somatic cancers. We demonstrate our functional assay correlates highly with current ClinVar classifications, and we provide functional evidence to reclassify 28 variants with uncertain interpretation in ClinVar. Further efforts to scale-up and expand the repertoire of data captured in the MCF10A spheroid assay will generate an informative dataset that will help accelerate the reclassification of PTEN VUS and aid in the development of treatments for PTEN-related cancers.

No potential conflicts of interest were disclosed.

Conception and design: J.T. Chao, D.M. Fowler, C.D. Roskelley, C.J.R. Loewen

Development of methodology: J.T. Chao, R. Hollman, W.M. Meyers, F. Meili, P. Dean, C.D. Roskelley, C.J.R. Loewen

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.T. Chao, R. Hollman, K.A. Matreyek, K. Haas, C.J.R. Loewen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.T. Chao, D.M. Fowler, C.D. Roskelley, C.J.R. Loewen

Writing, review, and/or revision of the manuscript: J.T. Chao, W.M. Meyers, K. Haas, C.D. Roskelley, C.J.R. Loewen

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.T. Chao, W.M. Meyers, K. Haas, C.J.R. Loewen

Study supervision: D.M. Fowler, C.D. Roskelley, C.J.R. Loewen

We thank Werner Chao, MASc, for advice on computational approaches and the UBC Life Sciences Institute FACS Facility for providing cell sorting services. This work was supported by grants from the Canadian Institutes of Health Research (PJT-152967) and the Simons Foundation (SFARI award #573845).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.