Learning with phenotypic similarity outperforms both conventional baseline and state-of-the-art models

We first trained naïve baseline models, one channel-specific support-vector machine (SVM) for variants from every single voltage-gated sodium channel (Dirac), and one global SVM for variants from all voltage-gated sodium channels (Union). We then implemented a method from our previous work: Multi-task learning (MTL) SVM embeds explicit prior knowledge on the structural and evolutionary pairwise similarity between ion channels in the feature space induced by our kernel function, which enables the model to assign more weight when learning from variants of similar channels [28]. MTL-SVM outperformed the baseline models across all metrics (Table 1). Finally, our novel multi-task multi-kernel learning model (MT)MKL approach for GOF/LOF prediction was trained with the Jaccard similarity measure and uniform kernel weights. Including phenotypic similarity information in the machine learning model further improves predictive performance (Table 1 and Fig 1). For comparison to current state-of-the-art methods, we trained a gradient boosting machine as used in previous work by Heyne et al. and the decision rule introduced by Brunklaus et al. [27,29]. For internal validation, these methods were trained and tested on the same data splits of our pre-processed training data using the k-fold nested cross-validation procedure described above. For external validation, all three models were used to obtain predictions for a list of novel missense variants in voltage-gated sodium channels (S1 Table). Results are shown in Table 2. In both validation procedures, our multi-task multi-kernel learning model performed best.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Phenotypic multi-task multi-kernel learning outperforms both conventional and state-of-the-art models.

a: Cross-validation model performance estimates of mean accuracy for each task (channel). b: ROC curves. The AU-ROC for each model is shown in Table 1. Abbreviations: MTL–multi-task learning; MKL–multi-kernel learning; ROC–receiver operator characteristics.

https://doi.org/10.1371/journal.pcbi.1010959.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Nested cross-validation (outer folds) model performance estimates for naïve baseline models (Dirac, Union), our previous state-of-the-art multi-task learning model (MTL), and phenotypic multi-task multi-kernel learning (MKL).

Metrics are reported as mean ± standard deviation (SD).

https://doi.org/10.1371/journal.pcbi.1010959.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Performance estimates for other state-of-the-art models and phenotypic multi-task multi-kernel learning (MKL).

For internal validation, each model was assessed by nested k-fold cross-validation. For external validation on a set of novel missense variants, the web-based tools for the Gradient boosting and decision rule method were used (see ‘Validation’). Note that AU-PRC and AU-ROC are not available for the decision rule as it is not a probabilistic classifier. Phenotypic learning outperforms other models across all metrics and is shown in bold. Abbreviations: ACC–accuracy; AU-PRC–area under the precision recall curve; AU-ROC–area under the receiver operator characteristics curve.

https://doi.org/10.1371/journal.pcbi.1010959.t002

Learning with phenotypic similarity is robust across different similarity measures and largely insensitive to phenotypic noise or sparsity

We compared three commonly used similarity measures: Jaccard, Lin, and Resnik. The Jaccard similarity measure always results in a positive semi-definite kernel matrix, which extends naturally to kernel-based learning methods [30]. Measures based on information content (IC), such as Lin and Resnik, lead to indefinite kernel matrices for the human phenotype ontology, as it allows for multiple parent terms. This invalidates the underlying theory and would lead to a non-convex SVM optimization problem, which is computationally expensive to solve and is not guaranteed to find the global optimal solution. For similarity-based classification, three methods are available to approximate the nearest positive semi-definite similarity matrix to be used as kernel matrix: spectrum clip, flip and shift [31]. The choice of matrix correction method did not have an impact on MKL model performance (Fig 2A and Table 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Phenotypic multi-task multi-kernel learning performs well across different similarity measures or matrix correction methods and is largely insensitive to sparse or noisy phenotypic information.

a: Dot chart of MKL model mean accuracy and standard deviation for each similarity measure and methods for approximating the nearest PSD matrix. The Jaccard index does not require matrix correction but is shown for comparison. b: Dot chart of MKL model mean accuracy and standard deviation for each similarity measure across levels of simulated phenotypic sparsity and noise. Noise is shown as human phenotype ontology term ratio (e.g. noise ratio 0.25 removes 3 out of 4 terms from the set, while noise ratio 4 adds 3 random terms for each original term in the set).

https://doi.org/10.1371/journal.pcbi.1010959.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Nested cross-validation (outer folds) model performance estimates.

Model performance of semantic similarity measures over methods of approximating the nearest PSD matrix. The Jaccard index does not require matrix correction but is shown for comparison.

https://doi.org/10.1371/journal.pcbi.1010959.t003

Information content measures are insensitive to varying taxonomic link distances and adapt the static knowledge structure of an ontology to different contexts by including probabilistic estimates [32]. We hypothesized that these measures may perform better in imperfectly phenotyped cases. Traditional phenotyping is expensive and prone to misclassification error. Recent studies have introduced ‘silver standard’ phenotyping from data mining of electronic health records and demonstrated advances in automated phenotyping pipelines [33,34]. We simulated sparse phenotypes by removing terms from the set of all terms for each variant at random before computing semantic similarity. As more terms were removed, the information from the phenotype kernel matrix was reduced and the mean accuracy of the MKL model approached that of the MTL model (Fig 2B and Table 4). Likewise, we simulated noisy phenotypes by sampling (without replacement) additional terms from the set of all possible terms and adding these to the set of terms for each variant before computing semantic similarity. We observed that all semantic similarity measures resulted in models with a stable accuracy even at high levels of simulated noise (Fig 2B and Table 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Model performance of semantic similarity measures across levels of simulated phenotypic sparsity and noise (see Fig 2).

Abbreviations: ACC–accuracy; AU-ROC–area under the receiver operator characteristics curve.

https://doi.org/10.1371/journal.pcbi.1010959.t004

Localized multi-task multi-kernel learning with hierarchical decomposition offers biological insight

Previously, we trained the model under the assumption that each kernel matrix (modality) contributes equally by assigning uniform global weights. The global weights for the best convex combination of kernel matrices can also be conveniently learnt from the data with wrapper methods, which have previously been demonstrated as promising tools to work with multi-source genomic data [35]. Here, we observed no significant difference in model performance between uniform weights and learnt global weights (Table 5). This is not surprising, as learning weights is most useful when we are working with a large set of kernel matrices and aim to achieve a sparse solution. If we have the prior knowledge to select a small set of adequate kernel matrices, uniform weights are usually sufficient. For theoretical guarantees, see Bach et al. [36].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Nested cross-validation (outer folds) model performance estimates for each multi-kernel learning method.

Uniform weights correspond to the sum of kernel matrices. For global weights, a wrapper method is used to learn the optimal weighted sum of kernel matrices. Hierarchical decomposition learns both localized weights and an optimal weighted sum of subset kernel matrices. The highest mean metric for each method is shown in bold.

https://doi.org/10.1371/journal.pcbi.1010959.t005

Global weights assign the same weight to a kernel for all pairs of samples. Therefore, localized weights for each task in the input space may be more appropriate for precision medicine approaches [37]. We implemented hierarchical decomposition multi-task multi-kernel learning, which makes efficient use of our prior knowledge on channel taxonomy to learn localized weights for each task and each subset of tasks that are leaves of subtrees rooted at each node [38]. Thus, we could identify tasks and subsets of related tasks where the weight of our phenotypic similarity kernel was high. These tasks were likely to contain phenotypic information that contributed towards classification between GOF and LOF, which may be useful for downstream analysis. A dendrogram of voltage-gated sodium channels and their local weights is shown in S3 Fig. The model correctly identified channels with known or shared genotype-phenotype correlations, e.g. SCN1A, SCN2A, SCN3A, and SCN5A. A weighted linear combination of these subset kernel matrices was then used as a composite kernel matrix for the model, where the weights were learnt with wrapper methods. Here, the weight vector favoured the kernel matrix of the subset of all tasks (root, wt = 0.93). Thus, the model recognized that voltage-gated sodium channels were similar enough to be learnt together. Overall, hierarchical decomposition offered a modest improvement in model performance (Table 5).

Finding difficult-to-predict variants in the reproducing kernel Hilbert space

Kernel-SVMs are considered ‘black box’ machine learning models. However, we may still be able to indirectly gather insight on model decisions. Fig 3A shows the histogram of projections for our model, as introduced by Cherkassky et al. [39]. Training data were projected onto the normal vector and plotted by their feature space distance to the optimally separating hyperplane. Samples within the range were support vectors. If a sample is close to the hyperplane, model prediction confidence decreases. The set of samples within an arbitrary distance of (“low confidence”) was compared to the set of samples outside this range (“high confidence”) by feature correlation with an unpaired two-samples Wilcoxon test. We adjusted p-values for multiple testing with the Benjamini-Hochberg method. Features that were significantly (α = 0.05) associated with variants that were hard to predict (“low confidence”) are shown in Fig 3B and 3C. This included variants in hotspots along the channel family sequence alignment corresponding to repeat domain IV, as well as variants with a high accessible surface area. For the final model, “low confidence” variants (n = 171) were predicted with a mean accuracy of 0.760, while “high confidence” variants (n = 204) were predicted with a mean accuracy of 0.931.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

a: Cherkassky histogram of projections for the phenotypic multi-task multi-kernel learning model. b-c: Features that were significantly different between variants predicted with low confidence (distance to hyperplane ≤ 0.5) and high confidence (distance to hyperplane ≥ 0.5). The accessible surface area feature is based on surrogate predictions by NetSurfP-2.0 and is given in Å [40].

https://doi.org/10.1371/journal.pcbi.1010959.g003

Data

Recently, Brunklaus et al. conducted a systematic literature search of variants in voltage-gated sodium channels that were functionally assessed in whole-cell patch-clamp recordings of mammalian cells [29]. We refer to their Supplementary Methods for the criteria used to assign functional labels to each variant. Of 437 variants, 62 variants were similar to wildtype function or did not have a clear overall gain- or loss-of-function. We removed these variants from further analysis. Our training data set included 375 unique non-synonymous missense variants across the 9 voltage-gated sodium channels (Na_V1.1-Na_V1.9). 164 variants were labelled as GOF, 211 variants were labelled as LOF.

Features

Each variant was identified by its channel, sequence position, and amino acid substitution. Sequence- and structure-based features were extracted and pre-processed as previously described [28]. These 62 features included amino acid physicochemical properties, radicality of substitution, paralog and ortholog conservation, secondary and tertiary structure, residue accessibility, disordered regions, and protein-protein interaction or binding sites (S2 Table). Variant-level information on primary disease and reference for phenotypes were extracted from the training data set and mapped onto Phenotype MIM numbers of the Online Mendelian Inheritance in Man (OMIM) [55]. Each OMIM ID is associated with a set of human phenotype ontology (HPO) terms which represent a standardized vocabulary of phenotypic features in human disease [44]. These mappings were requested via the HPO REST API (last accessed 22/June/2022 13:35 pm, release v2022-06-11) and resulted in an average of 9 HPO terms per variant (median 9, range 3–42 terms). All parents of the HPO terms in the directed acyclic graph from the parent ← child relationships were added to the respective HPO term sets, which resulted in 446 unique terms across the set of all variants, with each variant enriched to an average of 40 terms (median 38, range 12–226). Ontology analysis was performed with the ontologyX suite [56].

We considered three different measures for semantic similarity analysis. For the Jaccard index, we take the ratio of the number of elements in the intersection and the number of elements in the union of the two variant sets. A more rigorous definition and proof of positive definiteness of the function has been provided by Gärtner et al. [30]. For the Resnik measure, the similarity between two variants is given as the information content (IC) of their most informative common ancestor (MICA) term [32]. The IC is defined as −log₂(f) where f is the term frequency over all sets and the MICA is the common ancestor term with the highest IC. Lastly, the Lin measure extends the Resnik measure by incorporating the IC of the individual terms [57]. Further detail and a visual overview on phenotypic feature extraction and processing are shown in S1 Fig.

Both Resnik and Lin measures do not lead to positive semi-definite (PSD) similarity matrices in ontologies with multiple ancestors such as the HPO, where terms may have more than one parent term each. Before proceeding with kernel-based learning, we approximated the nearest PSD matrix with three methods suggested by Chen et al. [31]. In short, these are: i) clip, where each negative eigenvalue in our similarity matrix S is set to zero; ii) flip, where we flip the sign of each negative eigenvalue, which is equivalent to replacing the original eigenvalues of S with its singular values; iii) and shift, where we shift the spectrum of S by the absolute value of the minimum eigenvalue. For more rigorous definitions and an in-depth comparison of each method, see Chen et al. [31].

Models

When we consider phenotypic semantic similarity analysis as a kernel function and assert its positive semidefinite properties, we satisfy Mercer’s condition and can make implicit use of the reproducing kernel Hilbert space (RKHS) induced by our so-called phenotype kernel. We can avoid explicit representation of each variant in a Euclidean space and instead consider only their inner product, also thought of as pairwise similarity in the high-dimensional phenotypic feature space. This similarity-based learning represents a natural extension of semantic similarity analysis towards kernel-based machine learning methods.

For each variant, the vector of sequence- and structure-based features is transformed into pairwise similarities with a radial basis function (RBF) kernel K_n with x_j∈X [58]. Task similarity, i.e. the hierarchical taxonomy of voltage-gated sodium channels across the channel superfamily, is represented as an alignment-based tree distance measure and enables multi-task (parallel transfer) learning. Both methods have been previously described [28]. We can consider each of these three data input types of our variants (channel, sequence/structure, and phenotype) as different measurement modalities with individual kernel matrices for each modality. First, each kernel matrix is normalized in feature space:

Multiple kernel learning then allows us to represent these modalities in a single kernel matrix K, either by assigning uniform weights or by learning a linear combination of matrices K_m with weights β_m such that:

Global weights assign the same weight over the whole input space, which may not be desirable. Several methods of localized multiple kernel learning have instead been proposed [37]. We implemented hierarchical decomposition multi-task multi-kernel learning (MTMKL) by Widmer et al. [38]. This method makes use of prior knowledge in the form of a tree structure G, which in our setting is the taxonomy that relates voltage-gated sodium channels, where each channel is a task t_i corresponding to a leaf in G. Given the set of tasks T = {t₁, t₂,…,t_T}, let the subsets of tasks S_i⊆T be the leaves of the sub-tree rooted at each node n of G, that is leaves(n) = {l|l is descendant of n} [38]. The following definition for the kernel function is given:

A visualization of these kernel matrices is shown in S2 Fig. This method also offers additional insight [38]. First, the localized kernel weights of each task allow us to see which modality is important for which task. Secondly, we can refine our prior task similarity by learning latent task similarity from the data. If two tasks t_k and t_l are present in subsets S_i that are assigned weights β_i, we can define their latent similarity by summing up the weights of their shared task sets where the set of tasks sets containing task t_k is defined as , where T is the set of all tasks, such that the pairwise similarity of tasks γ_k,l is given as:

We trained a support-vector machine (C-SVM) on the resulting kernel matrix. Hyperparameter tuning was performed via grid search across cost C = {1e−4, 1e−2,…,1e4}, the RBF kernel hyperparameter σ = {1e−5, 1e−4,…,1e0}, and hierarchical multi-task learning baseline similarity parameter a = {1, 3, 5, 10, 100}. The choice of semantic similarity measure and PSD matrix approximation were included in the set of hyperparameters. Hyperparameter tuning and model assessment were done with nested k-fold cross validation, k = 5, where the model is trained and tested on different random sub-partitions of the data throughout multiple exhaustive iterations to obtain an estimate of the generalization performance. For illustrative purposes, a pseudocode and graphical representation of k-fold cross validation is provided in our previous work [28]. We report the following metrics: i) Accuracy, the number of correct classifications divided by the number of all classifications; ii) Sensitivity, or true positive rate (TPR); iii) Specificity, or true negative rate (TNR); iv) Matthews correlation coefficient (MCC), a symmetric correlation coefficient between true and predicted classes; iii) Cohen’s Kappa (κ), a measure of agreement between true and predicted classes; iv) F₁ score, the harmonic mean of precision and recall; v) The area under the receiver operator characteristics (AU-ROC) curve; vi) and the area under the precision-recall curve (AU-PRC).

Further model information is available in our AIMe Registry report (aime.report/xmSKHj).

Validation

For baseline comparison, we trained one SVM on variants from each task (channel-specific model, Dirac) and one SVM on variants from all tasks (global model, Union). For internal validation, we implemented gradient boosting as described by Heyne et al. [27], training and testing their method on the same splits of our data set. We also implemented the decision rule model proposed by Brunklaus et al., assigning labels based on missense pairs of identical paralog variants across the sodium channel family [29]. For external validation, we conducted a systematic literature search of PubMed with the search string (((SCN1A) OR (SCN2A) OR (SCN3A) OR (SCN4A) OR (SCN5A) OR (SCN8A) OR (SCN9A) OR (SCN10A) OR (SCN11A)) AND ((gain-of-function) OR (loss-of-function)) OR ((GOF) OR (LOF)) AND ((electrophysiology) OR (patch clamp) OR (voltage clamp))) AND (("2021/04/28"[Date—Publication]: "2022/07/13"[Date—Publication])). This query captures all publications after the last update of the training data set. 55 studies were identified and screened for missense variants in voltage-gated sodium channels with experimental evidence of whole-cell patch-clamp recordings that yielded an overall gain- or loss-of-function effect (n = 91). Variants that were present in the training data set, i.e. those known from previous functional experiments, were removed from further analysis (n = 61). Predictions were obtained from our own model and the web interfaces of Heyne et al. (funNCion, http://funncion.broadinstitute.org/), and Brunklaus et al. (SCN viewer, https://scn-viewer.broadinstitute.org/).

Software

This study was carried out in the R programming language, version 4.1.0, with RStudio, version 1.4.1106.