Analysis of nonlinear ordinary differential equations with the generalized Mittag-Leffler kernel

1.   Introduction

Various chemical modifications, including cytosine modification, uridine isomerization, and adenosine methylation, have been found in cellular RNA ^[1] and have been linked to important biological and physiological functions in cells ^[2]. У modification is a common posttranscriptional RNA modification known as the fifth base in RNA ^[3]. It is commonly present in a variety of species, and research has revealed that tRNA and rRNA contain large amounts of it ^[4]. Numerous biological processes have shown У to be crucial, and distinct У modifications serve different purposes at various places ^[5,6,7]. Therefore, the discovery of У sites in RNA sequences is crucial for both fundamental and applied biological research.

Initially, researchers identified У modification sites based on biochemical experiments. At first, researchers used paper chromatography to find У modification sites in the RNA of yeast, which was achieved by using RNA decomposition enzymes to decompose RNA and electrophoresis to separate out column chromatography on the upper layer of paper ^{[3,4,5,6,7,8]}. Later researchers successively used high-performance liquid chromatography and mass spectrometry to detect modification sites ^[9]. With the growing interest in this field, researchers have proposed a variety of high-throughput sequencing technologies, including Ψ-seq ^[10,11], PseudoU-seq ^[12] and CeU-Seq ^[13], and successfully used them to detect У sites. However, the methods described above are reliant on time-consuming, expensive, and difficult biochemical experiments, which are susceptible to environmental factors, and the sequencing process becomes increasingly difficult as the sequence length increases. Therefore, robust, fast, and inexpensive calculation methods are needed to predict У sites in RNA sequences.

First, Panwar and colleagues proposed a tRNAmod model to predict У sites in tRNA ^[14]. Then, a web server (PPUS) based on support vector machine (SVM) was proposed by Li et al. to identify У sites in S.cerevisiae and H.sapiens ^[15]. The frequency composition of nucleotides and the pseudo K-tuple nucleotide composition (PseKNC) were merged for feature representation in the iRNA-PseU model that Chen et al. created ^[16]. Subsequently, He et al. developed the SVM model (PseUI) to identify У sites in H.sapiens, S.cerevisiae and M.musculus, which combined a variety of feature extraction techniques including position-specific dinucleotide propensity (PSDP) ^[17]. Later, utilizing convolutional neural networks, Tahir et al. created a predictor (iPseU-CNN) ^[18]. Extreme gradient boosting (XGboost) was used by Liu et al. to create a new model known as XG-PseU ^[19]. Lv et al. also proposed a method called RF-PseU, which utilizes the LGBM algorithm for feature selection while combining the random forest algorithm for classification ^[20]. Saad et al. proposed a convolutional neural network model MU-PseUDeep ^[21], which combines sequence and secondary structural features to predict У sites. Li et al. built the model Porpoise by utilizing multiple type features and inputting them into the stacked ensemble learning framework ^[22]. Although the aforementioned techniques have proven successful in correctly identifying У sites in RNA sequences, they might still use more work in comparison to high-performance predictors ^{[23,24,25,26,27,28]}.

In this study, we build a У site identification model (iPseU-TWSVM) based on TWSVM, and Figure 1 depicts the model construction process. The model combines multiple feature representation methods, including Kmer, ENAC and EIIP. To obtain the best subset of features, the mRMR approach is utilized. The model is then evaluated using 10-fold cross-validation (10-CV) and independent testing (Ⅰ-testing). The average Ⅰ-testing accuracy of the iPseU-TWSVM is 3.4% higher than that of current advanced predictors, demonstrating the better generalization performance of our model. Therefore, iPseU-TWSVM may become an effective tool for У site identification.

2.   Materials and methods

2.1.   Datasets

In this work, we train and evaluate our models using datasets created by Chen et al. ^[29]. The steps of constructing the benchmark dataset were as follows: 1307 positive samples and 33,280 negative samples were obtained at first, and then the subset-balancing treatment was adopted to reduce the number of negative samples according to Euclidean distance. The obtained distance values were sorted in ascending order, and the first 1307 negative samples were selected to form the negative subset. The training datasets contained data from three species, namely, H.sapiens, S.cerevisiae and M.musculus. The H.sapiens training dataset included 495 positive samples and 495 negative samples; the S.cerevisiae dataset included 314 positive samples and 314 negative samples; and the M.musculus dataset included 944 samples, half of which were positive samples. There were just two species in the Ⅰ-testing datasets: H.sapiens and S.cerevisiae. Each of them included 200 samples, of which only half were positive and half were negative.

2.2.   Feature extraction

Different types of features reflect biological significance from different perspectives, including sequence composition and physicochemical properties. In this work, a variety of types of features are used to comprehensively consider the composition, distribution and physicochemical properties of nucleotides in the sequence from various aspects to further improve the prediction performance of subsequent work.

2.2.1.   Kmer

One effective technique for extracting RNA sequence characteristics is Kmer, which reflects the frequency of k adjacent nucleotides in the sequence. The frequencies of the k-neighboring nucleotides are used to generate the feature vector ^[30]. The method is provided by the web server Pse-in-One2.0 (http://bioinformatics.hitsz.edu.cn/Pse-in-One2.0/) ^[31].

2.2.2.   PC-PseDNC-General

The approach offers 22 different physicochemical properties to create the pseudo-dinucleotide composition ^[32,33,34]. It overwrites the local sequence order and the global sequence order information into the feature vector. The relevant features are expressed in this form:

with

where ${q}_{i}\left(i = \mathrm{1, 2}, \cdots, 16\right)$ represents the 16 dinucleotides' normalized frequency of occurrence; $\alpha \left(0\le \alpha \le 1\right)$ is the weight factor; and $\lambda$ is the highest counted rank. ${\rho }_{\mathrm{k}}$ is the k-tier correlation factor.

${\rho }_{\mathrm{k}}$ displays the relationship between the sequence orders of all neighboring dinucleotides along a specific RNA sequence, which can be written as

where $C\left({R}_{j}, {R}_{j+k}\right)$ indicates the correlation function expressed as

where parameter $\sigma$ is the number of physicochemical properties studied; ${P}_{g}\left({R}_{j}\right)$ and ${P}_{g}\left({R}_{j+k}\right)$ are the related values of the g^th property for the dinucleotides ${R}_{j}$ at position j and ${R}_{j+k}$ at position j+k.

2.2.3.   Nucleotide chemical properties (NCP)

The coding method reflects that each nucleotide in the sequence has different chemical structures and binding properties. The ring structures of the four RNA nucleotides (ACGU) differ from one another, hydrogen bond, and functional group. Based on these differences, they may be represented with a 3D coordinate ^[35].

2.2.4.   Accumulated nucleotide frequency (ANF)

The method incorporates data on each nucleotide's distribution in the RNA sequence as well as its frequency ^[35]. We can calculate the density ${d}_{i}$ of an RNA sequence's i^th prefix subsequence. It is defined as

where $i$ is the length of the sliding string and ${x}_{j}$ represents the nucleotide at the j^th position.

2.2.5.   Electron-Ion Interaction Pseudopotentials (EIIP)

The EIIP values represent the energy of the delocalized electron in the nucleotide. The nucleotides in the DNA sequence have previously been denoted by the EIIP values of A, G, C and T ^[36]. In the RF-PseU method ^[20], each nucleotide in an RNA sequence was also coded by EIIP feature vectors.

2.2.6.   Enhanced Nucleic Acid Composition (ENAC)

Using a fixed length window, the approach was used to determine the nucleotide composition ^[20,35]. Afterward, RNA sequences were converted into equally long feature vectors. Sequence length and sliding window size are two factors that affect the dimension of ENAC coding.

where N is the sliding window size and n is the coding dimension.

2.2.7.   Nucleotide binary profiles (NBP)

Binary profiles provide the position specific composition of nucleotides in RNA fragments ^[35,36]. A four-digit binary vector is used to encode each nucleotide. Dibinary profiles are different from binary profiles in that they are encoded for 16 dinucleotides, i.e., AA is denoted by (0, 0, 0, 0).

2.3.   Feature selection

2.3.1.   Max-relevance and min-redundancy (mRMR)

mRMR ^[37] is a commonly used feature selection method for compressing feature vector space. The goal of this technique is to identify a subset of features from the initial feature set that have the lowest correlation between features and the highest correlation with the output result. It considers the connection between features as well as the association between features and labels. The mechanism of feature selection is as follows.

The mutual information is used to find the feature subset S containing m features first, so that the m features found have the maximum correlation with the category c. The correlation between the feature subset S and the category c is defined by the average value of all mutual information between each feature and category as shown in (7).

where $I({x}_{i}; c)$ is mutual information; $S$ is a subset of features of length m; ${x}_{i}$ is the i^th feature in $S$ and $c$ is category variable.

Then the features selected by the maximum correlation may be redundant, and (8) is used to eliminate the redundancy among m features.

The final feature subset S is obtained by combining the maximum correlation D with the minimum redundancy R.

Compared with other feature selection methods, the proposed algorithm considers the redundancy among features, further optimizes the feature subset, and solves the problem that the maximum dependency is difficult to achieve. However, only approximate optimal solutions can be obtained in practical applications.

2.4.   Classifiers

2.4.1.   Twin support vector machine (TWSVM)

Consider the binary classification issue using the training datasets

where ${u}_{i}\in {R}^{n}, i = \mathrm{1, 2}, \cdots, m+n.$

Let $T = {\left({u}_{1}, {u}_{2}, \cdots, {u}_{m}\right)}^{T}\in {R}^{m\times n}, F = {\left({u}_{m+1}, {u}_{m+2}, \cdots, {u}_{m+n}\right)}^{T}\in {R}^{n\times n}$ and $l = m+n$.

TWSVM ^[38] looks for a pair of nonparallel hyperplanes in the linear case.

where ${w}_{+}\in {R}^{n}, {w}_{-}\in {R}^{n}, {b}_{+}\in R, {b}_{-}\in R$ by solving the following pair of QPPs:

and

where ${c}_{1}, {c}_{2}$ are the penalty parameters, ${e}_{+}, {e}_{-}$ are all 1 vectors (${e}_{+}, {e}_{-} = {[1\cdots 1]}^{T}$) whose dimensions are the same as the number of positive and negative samples respectively, and ${\xi }_{+}, {\xi }_{-}$ are slack vectors of appropriate dimension.

Minimizing the objective function means making a hyperplane as close as possible to one type of data, and the constraint requires that the distance between the hyperplane and the other type of data is at least greater than 1. Their corresponding Lagrange dual problems are

and

where

The solution to the primary problem can be acquired by addressing the dual problem, which can be obtained by

Therefore, an unknown point $u\in {R}^{n}$ is predicted to the Class by

where $\left| \cdot \right|$ is the perpendicular distance of point $u$ from the planes ${w}_{s}u+{b}_{s} = 0, s = -, +.$

This method not only divides a large quadratic programming problem into two small quadratic programming problems, which improves the training speed, but also is not very sensitive to noise.

2.5.   Model evaluation metrics and methods

Five indicators were widely used to assess how well the built models performed ^[39,40,41], accuracy (ACC), sensitivity (SN), specificity (SP), Matthew correlation coefficient (MCC), and integral area under the receiver operating characteristic curve (auROC), which were calculated using the following equations.

where TP, TN, FP, and FN represent true positive, true negative, false positive and false negative, respectively.

We use 10-CV for comparison ^{[42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62]}. The training datasets are equally divided into ten subsets. The remaining one subset is tested after the proposed model has been trained using nine subsets. After each subset is tested once, the procedure is repeated ten times, and the average results represent the final performance. Finally, Ⅰ-testing was used in the testing datasets to evaluate the training model.

3.   Results and discussion

3.1.   Comparison with different feature combinations

Feature extraction affects the results of subsequent sequence classification. To obtain better performance, this paper studied seven different features, including Kmer, PC-PseDNC-General and ANF, EIIP, ENAC, NCP + NBP from the RF-PseU method ^[20]. These features were first used in the experiment separately, and then multiple features were selected for different combinations according to the test results to obtain better experimental results.

Table 1 lists the results of feature combination for the H_990 dataset using the TWSVM methods. The first six rows are the performance of single features, with Kmer, EIIP, ENAC and NCP + NBP returning the best results, which are roughly distributed in the range of 0.57–0.59. Since the test results of single features were lower than those of the RF-PseU predictor, we combined several features that perform well and used the mRMR method to select the best feature for model construction. For the H_990 dataset, the combined characteristics listed in Table 1 have four results, including Kmer + ENAC, Kmer + EIIP, Kmer + NCP + NBP and Kmer + PC-PseDNC-General + ANF + EIIP + ENAC + NCP + NBP. Using the TWSVM method, the result of combined features was usually approximately 1–3% higher than that of single features, with a maximum ACC value of 0.65 and the best feature combination being KMER + PC-PseDNC-Generel +ANF + EIIP + ENAC + NCP + NBP. Compared with the RF-PseU predictor, the accuracy was improved by 0.7%. We chose this feature combination for the 10-CV of the training set H_990 and applied it to the Ⅰ-testing of the testing set H_200.

Table 2 lists the test results of different feature combinations on the S_628 dataset using the TWSVM method. Similar to the results in Table 1, the test results of Kmer, EIIP, ENAC and NCP + NBP were better. We combined those features that reported better performance. Except for Kmer + NCP + NBP, the results of other combined features were improved compared with the single feature results. Among them, the performance of the feature combination Kmer + PC – PseDNC – Generel + ANF + EIIP + ENAC + NCP + NBP was the best, with test results improved by approximately 6% compared to other combinations. We also used the independent test set S_200 to test this feature combination.

Table 3 shows the test results for the dataset M_994. The feature combination Kmer + PC-PseDNC-Generel + ANF + EIIP + ENAC + NCP + NBP had the best performance, and all test indicators were higher than the rest of the feature combinations. The combination was increased by approximately 5%, the MCC was increased by approximately 10%, and the AUC was also significantly improved.

3.2.   Comparison of different feature selection methods

In this study, we contrast the mRMR approach and the LightGBM method ^[64] since feature selection is a crucial component of model construction. Figure 2 shows their accuracy on the training datasets of the three species. The findings demonstrate that the performance of the mRMR technique is superior, which further enhances the classification accuracy of the model. The accuracy of the mRMR method on the three species is greater than that of the LightGBM approach.

3.3.   Optimization with different feature dimensions

The accuracy of classification results may be successfully increased by feature selection. We initially utilized the mRMR technique to pick feature subsets with high correlation with class labels and low feature redundancy to obtain the optimum feature dimension. To further obtain the feature dimension with the best precision, the incremental feature selection approach was applied. After many experiments, we found that the accuracy of Ⅰ-testing and 10-CV fluctuates as the number of characteristics rises and the highest accuracy mostly appeared within 100 or 150 dimensions, as illustrated in Figure 3. The accuracy of each species initially increases rapidly as the feature dimension increases, and then fluctuates continuously. For H.sapiens species, the highest 10-CV accuracy of 0.65 was obtained when the feature dimension reached 33, while the highest independent test accuracy of 0.763 was obtained at relatively low dimensions. The highest 10-CV accuracy of 0.722 and independent test accuracy of 0.825 for S. cerevisiae species were between 60–80 dimensions, obtained in 72 and 76 dimensions, respectively. M. musculus species only showed 10-CV results with the highest value at a feature dimension of 62.

Tables 4–6 show the changes of feature dimensions after feature selection and the distribution of feature subsets after optimization for the three species. It can be found that ENAC and EIIP occupy a large number in the optimized feature subset of the three species, followed by NCP and NBP, ANF and Kmer occupy a small number, and there is no PC-PseDNC-General in the optimized feature subset. It indicates that each feature has different contributions in the model, and ENAC and EIIP play an important role in the model.

3.4.   Comparison with SVM classifier

Since many previous researchers built У sites recognition models based on support vector machines, we employed SVM ^[65] as a classifier in the same feature space to compare the performance of TWSVM with that of SVM. Figure 4 displays how it performed. The ACC, MCC, and AUC based on the TWSVM model were found to be larger than those based on the SVM model for the 10-CV results of the three species, while the independent test results may have more clearly indicated the difference between the two. All of the evaluation metrics outperformed the SVM model. As a result, we concluded that the TWSVM model performs much better than the SVM model, suggesting that it may be better suited for identifying У sites in RNA sequences.

3.5.   Comparison with previous methods

The effectiveness of iPseU-TWSVM was also evaluated in comparison to other advanced predictors, such as iRNA-PseU ^[16], PseUI ^[17], iPseU-CNN ^[18], XG-PseU ^[19] and RF-PseU ^[20]. The 10-CV and Ⅰ-testing results of the advanced У site predictors using iPseU-TWSVM are contrasted in Tables 7 and 8, respectively. The 10-CV results reveal that the accuracy of iPseU-TWSVM on H.sapiens is 1.7% less accurate than that of the best predictor iPseU-CNN on this species, and the accuracy on S.cerevisiae and M.musculus is 0.722 and 0.728, respectively, which is 2.6 and 2.0% less accurate than that of the best predictor RF-PseU. Although iPseU-TWSVM does not perform optimally on the training set, iPseU-TWSVM has higher accuracy than other predictors on all species in terms of Ⅰ-testing. H.sapiens and S.cerevisiae are 1.3 and 5.5% more accurate in Ⅰ-testing than the best predictor RF-PseU, with corresponding accuracy values of 0.763 and 0.825, respectively. We also calculated the average accuracy of several species so that we could compare the predictors' performance in depth. As shown in Table 9, the 10-CV accuracy of iPseU-TWSVM is 1.3% less than that of RF-PseU. In terms of Ⅰ-testing, iPseU-TWSVM is significantly improved by 3.4% compared with RF-PseU. iPseU-TWSVM performs much better overall than the other predictors. The findings demonstrate that iPseU-TWSVM, a very practical technique, has greater generalization performance and is more appropriate for recognizing У sites in RNA sequences.

4.   Conclusions

This work proposes the use of a novel model called iPseU-TWSVM to identify RNA У sites across various species. We have used an efficient feature selection method to obtain the best feature subset and selected TWSVM as the classifier to increase recognition accuracy. Finally, we compared advanced predictors and found that iPseU-TWSVM significantly improved the independent test accuracy by 3.4%, while the accuracy of cross validation was lower by 1.3%. Through comprehensive analysis, it was concluded that the relatively poor performance of the training datasets was due to the following two reasons. One is that the features used by the best predictor are different, and the other is that the classifier of the model is different. The above results indicate that iPseU-TWSVM had better generalization performance and could more accurately identify У sites from RNA sequences. It is anticipated that iPseU-TWSVM will be effective in identifying RNA У sites.

The contribution of this work has the following three aspects: (ⅰ) the model uses TWSVM as a classifier, which improves the accuracy of the model and improves the training speed; (ⅱ) the model has good generalization performance and can be applied to the prediction of other sites in the sequence; (ⅲ) further accurate identification of У sites in the sequence lays the foundation for disease control and related drug development. At the same time, this work also has the following shortcomings: (ⅰ) in the feature selection part, only two algorithms are compared, and subsequent research can try other algorithms to further improve the feature subset; (ⅱ) the model uses TWSVM as a classifier. In the original problem of TWSVM, only empirical risk is minimized, but structural risk is not minimized. Moreover, the algorithm can only obtain approximate solutions. Subsequent research can consider improving TWSVM or try other classification algorithms as the classifier of the model to improve the prediction performance. Future work will study emerging methods ^{[66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81]} to further improve the accuracy of the model.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (NSFC 62172076, 62072385), and the Municipal Government of Quzhou under Grant Number 2020D003 and 2021D004.

Conflict of interest

The authors declare there is no conflict of interest.