An application of topological data analysis in predicting sumoylation sites

Dataset description

In this study, we considered two sumoylation site datasets. One was retrieved from GPS-SUMO (Zhao et al., 2014) and consisted of 510 proteins with 912 annotated sumoylation sites; this dataset was named dataset1. This dataset was designed to explore the utility and interpretability of the features constructed from TDA for identifying sumoylation sites. The other dataset, named dataset2, was obtained from iSumok-PseAAC (Khan et al., 2021) with 4,987 annotated sumoylation sites from 1,311 proteins. The applicability of our proposed model was validated using dataset2. As not all complete atomic coordinates of every protein were accessible, 471 proteins (Data S1) and 1,288 proteins (Data S2) with their complete coordinates and sequence information, were retained in dataset1 and dataset2, respectively. Relevant information is available on the AlphaFold database (Varadi et al., 2022) and the UniProt database (Wang et al., 2021).

The preprocessing procedures for each dataset were as follows. First, each lysine K from a protein was characterized by a peptide, P, whose length was selected based on the work of Khan et al. (2021). P was composed of 20 upstream and downstream residues, respectively, with K as the center. Missing residues were added dummy code X. The peptide, P, was considered to be a positive sample if its center, K, was experimentally annotated as a sumoylation site; otherwise, it was considered to be a negative sample. Further, we mixed the positive and negative samples and computed the pairwise sequence identity to avoid the bias of homology. If the sequence identity between two given samples was more than 40%, only one of them was retained while the other was ignored. Ultimately, the samples that were not redundant were divided into a positive subset or a negative subset according to the category to which they belonged. After going through these steps, we obtained 775 positive samples (Data S3) and 17,807 negative samples (Data S4) from dataset1, and 4,493 positive samples (Data S5) and 24,456 negative samples (Data S6) from dataset2.

Preparations

The features constructed from TDA for identifying sumoylation sites were obtained as follows. First, the original protein data was represented in an all-atom model. Then, simplicial complexes were constructed according to the represented data. Afterwards, the PH analysis was conducted to reveal the topological information of the proteins. Finally, features were extracted from persistence barcodes (PBs), which could visualize the results of the PH analysis. Additional information on TDA and PH can be found in Edelsbrunner, Letscher & Zomorodian (2000), Zomorodian & Carlsson (2004), Munkres (2018), and Wang, Cang & Wei (2020).

Representation of protein data

We used an all-atom model when dealing with the protein fragments. An all-atom model includes various types of atoms, such as O, C, N, S, and P, which are all of equal importance. Note that the hydrogen atoms were ignored during the PH analysis, as they created redundant barcodes that did not contribute much to feature construction. Figure 3 shows two protein fragments of O00429 (UniProt Entry) and their corresponding all-atom models. The central lysine of each fragment is K532 (sumoylation site) and K92 (non-sumoylation site), respectively.

Simplicial complexes

Figure 3 illustrates that the represented protein fragments are essentially the point cloud data of dimension 3. The Vietoris-Rips (VR) complex and Alpha complex were considered to characterize the point cloud data in this work. The following definitions also can be seen in Pun, Lee & Xia (2022).

Let X be a finite point set in Euclidean space ℝⁿ. The VR complex of X with parameter ϵ is the set of all σ⊆X, such that any pairwise distance of its points is at most 2ϵ.

To introduce the Alpha complex, we need some related concepts. Given a good cover U of X, i.e., X⊆∪_i∈IU_i. The nerve of U is defined as: (1)${\displaystyle N\left(U\right)=\left\{J\subseteq I\mid {\cap }_{j\in J}{U}_j\ne 0̸\right\}\cup 0̸.}$ A closed ball in ℝⁿ with center, x, and radius, δ, is denoted as B(x, δ). The union of closed balls with center points in X forms a cover of X, and the corresponding nerve creates a simplicial complex named the Čech complex, (2)${\displaystyle C\left(X,\delta \right)=\left\{\sigma \subseteq X\mid {\cap }_{x\in \sigma }B\left(x,\delta \right)\ne 0̸\right\}.}$

Given a point, x ∈ X, the Voronoi cell of x is defined as: (3)${\displaystyle V_x=\left\{y\in {\mathbb{R}}^n\mid |y-x|\le |y-x^{{}^{\prime }}|,\forall x^{{}^{\prime }}\in X\right\}.}$ The Voronoi diagram is the collection of all Voronoi cells. The dual graph of the Voronoi diagram forms a simplicial complex called the Delaunay complex. Let R(x, δ) be the intersection of the Voronoi cell, V_x, with the ball, B(x, δ), that is, R(x, δ) = V_x∩B(x, δ). The Alpha complex of X is defined as the nerve of cover ∪_x∈XR(x, δ), i.e.: (4)${\displaystyle A\left(X,\delta \right)=\left\{\sigma \subseteq X\mid {\cap }_{x\in \sigma }R\left(x,\delta \right)\ne 0̸\right\}.}$ Intuitively, the Alpha complex is a subcomplex of the Delaunay complex.

To illustrate the theory of simplicial complexes, we consider the example showing a set of points, S, of dimension 2 (Fig. 4A). With the increasing radius, the simplices contained in the two complexes have distinct differences. The VR complex is entirely determined by its 1-simplices, that is, if all the 1-faces of a simplex are in the VR complex, then so is the simplex. However, the Alpha complex is only suitable for subjects in Euclidean spaces, and its construction is more complicated. Based on the same dataset, diverse simplicial complexes can be constructed according to different rules. Therefore, it may be of value to combine these complexes for the PH analysis, as they may reveal different information about the same protein data.

(Element specific) persistent homology analysis

Typically in PH, a nested sequence of subcomplexes is constructed based on a filtration parameter. After analyzing the homology of each subcomplex, the topological information is characterized by different persistence time of homological generators with respect to the filtration parameter.

Given a simplicial complex, K, the filtration of K is a nested sequence of subcomplexes of K, i.e.: (5)${\displaystyle 0̸=K_0\subseteq K_1\subseteq \cdots \subseteq K_n=K.}$ The p-persistent k-th homology group at filtration time i can be represented as: (6)${\displaystyle H_k^{i,p}=Z_k^i/\left(B_k^{i+p}\cap Z_k^i\right),}$ where $Z_k^i$ is the k-th cycle group of K_i and $B_k^{i+p}$ is the k-th boundary group of K_i+p. The rank of $H_k^{i,p}$ is called the p-persistent k-th Betti number of K_i, denoted as ${\beta }_k^{i,p}$.

In the PH results, each topological generator is characterized by a pair of values that record when it appears and dies, named birth time (BT) and death time (DT), respectively. The PH results can be visualized as PBs using the endpoints of a bar to represent the BT and DT of each generator, respectively. Figure 5 shows the PBs of dimension 0, 1, and 2 of K532 and K92 from the filtration of the VR complex, respectively. With the increasing filtration values, atoms are pairwise connected according to different interactions of the protein. These connections induce the death of generators of dimension 0. More connections indicate a greater occurence of simplices of higher dimensions, which are related to the 1- and 2-bars. There is a clear difference between the 0-bars in [1.25, 1.5] Å of Figs. 5A and 5B that reflects the different interaction patterns of these fragments. Moreover, more types of 1- and 2-bars can be found in the lower right panels, indicating that K92 has a more complicated spatial structure, which is consistent with Fig. 3.

We adopted the element specific persistent homology (ESPH) analysis to reveal additional information about the biochemical and physical properties of the proteins. ESPH analyzes biomolecular data by specifying one or more types of elements, followed by the PH analysis (Meng et al., 2020). Distinguishing element types enables us to reduce biomolecular complexity and retain critical properties of the studied data (Cang & Wei, 2018). In this work, element C and N were considered for ESPH, and the software package GUDHI (Project, 2021) was used for the (ES) PH analysis.

Features based on topology

For a peptide sample, P, the features of P were divided into four categories according to different vectorization methods and filtrations as follows:

• TF1: features based on the 0-, 1-, and 2-bars of P from the filtration of the VR complex.

• TF2: features based on the 1- and 2-bars of P from the filtration of the Alpha complex.

• TF3: features based on each residue of P (excluding the central lysine K) from the filtration of the VR complex.

• TF4: features based on the local region of P from the filtration of the VR complex.

In TF1, 1-bins (unit: Å) were used to vectorize 0-, 1-, and 2-bars. The endpoints of the bars were taken from two consecutive values in the specific list for different dimensions. These lists were set to [1.2, 1.3, 1.4, 1.5, 1.6, 2.0], [1.5, 2.7, 3.5, 4.5, 5, 6.7], and [2.4, 2.9, 5.5, 6.7], respectively. For example, (1.2, 1.3) Å was used as a bin to vectorize 0-bars of P. Moreover, TF1 included the second and third longest bar lengths of dimension 0, the sum and mean of bar lengths of dimension 0, the onset value of the longest 1-bar, and the barcode statistics of 1- and 2-bars. Therefore, TF1 contained 48 features. The barcode statistics of 1- and 2-bars from the filtration of the Alpha complex resulted in 30 features in TF2.

In TF3, for each residue of P, the (1.25, 1.5] Å and (1.5, 1.75] Å bins were used. TF3 also contained the number of 0-bars, and the summation of bar lengths of dimension 0, 1, and 2. Note that each feature of the dummy code X was set to 0. Hence, the number of features in TF3 was equal to 240.

In TF4, the local region of P refered to the peptide fragment from the 2-th upside residue to the 2-th downside residue, according to the central lysine K. For the PH analysis, 1-bins were taken from the (1.2, 1.6) Å split by the fixed bin size 0.1 Å, and the barcode statistics were used to vectorize the 1-bars. For the ESPH analysis of element type C, 1-bins were taken from the (1.5, 3) Å split by the fixed bin size 0.5 Å, and the barcode statistics were used to vectorize the 1-bars; for the ESPH analysis of element type N, the number of 0-bars whose death times were less than 10 Å was considered. TF4 contained 38 features. All of these four categories gave rise to a total of 356 features for P.

Model training and validation

There was an imbalance of the sample size of positive and negative classes for each dataset, which may affect the learning process. To address this issue, an undersampling method named NearMiss (Lopez et al., 2020) was used. The “imbalanced-learn” package (Lemaître, Nogueira & Aridas, 2017) of Python was employed to balance each dataset. After undersampling, 775 and 4,493 negative samples were selected from dataset1 and dataset2, (Data S7 and S8), respectively.

The features constructed from TDA were fed into various binary classifiers, including the gradient boosting classifier (GBC), random forest classifier (RFC), and support vector classifier (SVC). GBC uses a boosting algorithm to make up for the shortcomings of the original model by building a weak learner at each step of the iteration. RF is an ensemble learning method based on the bagging algorithm, which independently integrates different decision trees during training. SVC is a classifier based on the support vector machine algorithm whose decision boundary is the maximum-margin hyperplane solved for the learning samples. All of these classifiers were implemented here with the “scikit-learn” package (Pedregosa et al., 2011) of Python. The “n_estimators” parameter was set to 370 and 500 for RFC and GBC, respectively. Moreover, for RFC, the “oob_score” parameter was “True”, and the “max_features” parameter was chosen as “sqrt”; for SVC, the “gamma” parameter was set to “auto”, and the “probability” parameter was “True”. All other parameters used their default values.

The K-fold cross validation and independent set test were adopted to evaluate the model performance here. The K-fold cross validation splits the original dataset into K disjoint subsets. Each subset is selected in turns for testing, while the rest of parts are used for training. Note that each data participates in training in the K-fold cross validation. However, the independent set test divides the original dataset into training and testing subsets at a given dividing ratio, where samples in the testing subset are not participants in model training.

Evaluation metrics

To evaluate the performance of predictors from different perspectives, we adopted multiple evaluation metrics as follows: (7)${\displaystyle \text{Sp}=\frac{\text{TN}}{\text{TN}+\text{FP}},}$ (8)${\displaystyle \text{Sn}=\frac{\text{TP}}{\text{TP}+\text{FN}},}$ (9)${\displaystyle \text{Acc}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}},}$ (10)${\displaystyle \text{MCC}=\frac{\left(\text{TP}\times \text{TN}\right)-\left(\text{FP}\times \text{FN}\right)}{\sqrt{\left(\text{TP}+\text{FP}\right)\left(\text{TP}+\text{FN}\right)\left(\text{TN}+\text{FP}\right)\left(\text{TN}+\text{FN}\right)}},}$ where “T” and “F” denote the true and false cases of prediction, “P” and “N” denote the positive and negative classes, respectively. Specifically, TP (true positive) counts the correctly predicted sumoylation sites, and TN (true negative), FP (false positive), and FN (false negative) are defined similarly.

The performance of predictor was also measured using the area under the receiver operating characteristic (ROC) curve. Different pairs of true positive rate (TPR, also known as Sn) and false positive rate (FPR, defined as follows) can be obtained by adjusting the classification threshold of a given classifier, which are the data points of ROC. The AUC of ROC refers to the probability that a classifier outputs a higher probability for the given positive sample being positive than for the given negative sample being positive, when randomly given a positive and a negative sample. It reflects the sorting ability of a classifier. A higher AUC implies a better classifier. (11)${\displaystyle \text{FPR}=\frac{\text{FP}}{\text{FP}+\text{TN}}.}$

In addition to evaluating model performance, we sought to determine the topological features that contributed to the prediction of sumoylation sites. The F-score value (Xu et al., 2016) measures the ability of features to distinguish among two classes. The F-score value of the ith feature is defined as: (12)${\displaystyle F_i=\frac{{\left({\overline{{\phi }_i}}^{+}-\overline{{\phi }_i}\right)}^2+{\left({\overline{{\phi }_i}}^{-}-\overline{{\phi }_i}\right)}^2}{\frac{1}{n_{+}-1}\sum _{k=1}^{n_{+}}{\left({\phi }_{k,i}^{+}-{\overline{{\phi }_i}}^{+}\right)}^2+\frac{1}{n_{-}-1}\sum _{k=1}^{n_{-}}{\left({\phi }_{k,i}^{-}-{\overline{{\phi }_i}}^{-}\right)}^2}\left(i=1,2,\dots ,\Omega \right),}$ where ${\phi }_{k,i}^{+}$ and ${\phi }_{k,i}^{-}$ are the ith feature value of the k-th positive sample and k-th negative sample, respectively. ${\overline{{\phi }_i}}^{+}$ is the mean of all ${\phi }_{k,i}^{+}$, i.e., ${\overline{{\phi }_i}}^{+}=\frac{1}{n_{+}}{\sum }_{k=1}^{n_{+}}{\phi }_{k,i}^{+}.{\overline{{\phi }_i}}^{-}$ can be obtained similarly. $\overline{{\phi }_i}$ denotes the mean of the ith feature values of all samples, that is, $\overline{{\phi }_i}=\frac{1}{n_{+}+n_{-}}\left({\sum }_{k=1}^{n_{+}}{\phi }_{k,i}+{\sum }_{l=1}^{n_{-}}{\phi }_{l,i}\right)$. The higher F-score value means the greater contribution to the classification. The codes of our work are available on the GitHub repository using the link: https://github.com/Xiaoxi-Lin/SUMO_TOP.git.