The cornerstone of a robust and reliable model is to generate a reliable and strict benchmark dataset. In previous literature, scholars used 50°C as a cutoff to construct a benchmark dataset. However, this criterion did not seem objective because proteins might be stable even above the OGT of microorganisms. For instance, a protein produced by microorganisms living at 45°C is likely not to denature at 60°C. According to the 50°C cutoff criterion, this protein is included in the negative dataset, but it should be included in the positive dataset as it is still stable above the 50°C. To eradicate this effect as much as possible, we used Lin and Chen’s (2011) strict and objective standard to generate a benchmark dataset. According to Lin and Chen’s (2011) criterion, the proteins in the microorganism with OGT > 60°C and <30°C were regarded as thermophilic and non-thermophilic proteins, respectively. Of course, even after using Lin and Chen’s (2011) criterion, the effect mentioned still exists but not as strongly as when compared to the 50°C cutoff criterion. All protein sequences were extracted from a universal protein resource (UniProt). Subsequently, the following steps were used to ensure the quality of protein data: (I) the proteins which have been manually reviewed remained; (II) proteins containing ambiguous residues were excluded; (III) sequences which are a fragment of other proteins were excluded; (IV) proteins which infer from prediction or homology were excluded; (V) to remove redundancy and homology bias, CD-HIT program (Huang et al., 2010) was used by setting a cutoff of sequence identity to 30%. As a result, the final benchmark dataset contained 1,443 non-thermophilic and 1,366 thermophilic proteins. Our final dataset contains only a few thousand proteins because the growth temperature of some microorganisms is known (Li G. et al., 2019) and UniProt contains few confirmed proteins. We only included experimental data. Moreover, noise and redundancy were removed, which also caused a reduction in the number of proteins. For training model, the dataset was divided into 80:20 ratios; model was trained on 80% dataset and validated on 20% dataset.

Protein sequences were transformed into numerical vectors to identify thermophilic proteins by machine learning methods (Liu et al., 2019, 2020; Li et al., 2021; Zhang et al., 2021a,b). To accomplish this task, we used the iFeature program (Chen et al., 2018) to generate seven kinds of protein features, namely amino acid composition (AAC), traditional pseudo amino acid composition (tPseAAC), amphiphilic pseudo amino acid composition (aPseAAC), the composition of k-spaced amino acid pairs (CKSAAP), dipeptide composition (DC), dipeptide deviation from the expected mean (DDE), and composition, transition, and distribution (CTD). These features will be described in detail in the following sections.

Redundant features and noise affect the prediction performance of the model. In order to get the best prediction performance, it is necessary to remove redundant features and noise using feature selection methods (Tang et al., 2020; Zhang Z. M. et al., 2020; Dao et al., 2021c). In this study, the analysis of variance (ANOVA; Tang et al., 2018) was applied for feature ranking, and a sequential backward selection strategy was used to pick out optimal features. The following section will introduce the method briefly.

Analysis of variance (ANOVA) can be used to select the best feature subsets based on F-value. F-value is the ratio of the variance between the sample types and the variance within the samples. A feature’s greater F-value implies that the feature can contribute more to discriminating between positive and negative samples.

F-value for a feature m can be calculated as:

where s²_b is the variance between the features and s²_w is the variance with each feature’s sample. These variances can be represented as:

where K denotes the total features, N denotes the total samples, fij(m) denotes the m-th feature of the j-th sample in the i-th group, and n_i denotes sample in the i-th group. The degree of freedom for between features df_b and within features df_w was K-1 and N-1, respectively. Detailed descriptions about ANOVA can be referred to as reference (Tang et al., 2018).

For classification, we examined a number of classifiers, including Support Vector Machine (SVM; Tang et al., 2017), K Nearest Neighbor (KNN; Zuo et al., 2013; Zulfiqar et al., 2021a), Random Forest (RF), and Multi-layer Perceptron (MLP) for training the model. The following sections will introduce these classifiers briefly.

In order to evaluate the overall model performance, the following parameters were used (Lv et al., 2020b,d; Shao et al., 2021).

where Sn, Sp, Acc, and MCC denote sensitivity, specificity, accuracy, and Matthews’s correlation coefficient. Thermophilic proteins classified as thermophilic were denoted TP (true positive), Non-thermophilic proteins classified as non-thermophilic were denoted TN (true negative), Non-thermophilic proteins classified as thermophilic were denoted by FP (false positive), and thermophilic proteins classified as non-thermophilic were denoted by FN (false negative).

For performance evaluation, seven descriptors including AAC, tPseAAC, aPseAAC, DC, DDE, CKSAAP, and CTD were used to create numerical vectors from protein sequences. In order to use these numerical vectors, MLP-based models were trained to evaluate their performances. Results showed that the AUC are 0.9723, 0.9551, 0.9519, 0.8812, 0.9081, 0.9081, and 0.9786 for AAC, tPseAAC, aPseAAC, DC, DDE, CKSAAP, and CTD, respectively (as shown in Table 2). In order to remove the redundant features and improve the prediction performance of the model, a feature selection method should be used to pick out the optimal features from each descriptor. In this work, ANOVA was used to rank features for selecting the best feature subsets from the seven types of descriptors. Table 2 also recorded the performance of each descriptor after feature selection. It showed that AAC, tPseAAC, aPseAAC, DC, DDE, CKSAAP, and CTD produced the best AUC of 0.9735, 0.9580, 0.9610, 0.9143, 0.9165, 0.8349, and 0.9644, respectively. Obviously, the performance of each descriptor increased after the feature selection except the CTD descriptor; therefore, we considered all features of CTD in our study.

The above results and analysis have demonstrated that each descriptor has useful information to discriminate thermophilic proteins from non-thermophilic proteins. We adopted a feature fusion strategy to include the valuable information of all selected features from each descriptor in model training. In feature fusion, the selected optimal feature subsets of seven descriptors were fused and inputted into the MLP classifier to distinguish thermophilic proteins from non-thermophilic proteins. Table 2 shows that the AUC increased to 0.9864, suggesting that feature fusion is very effective and has made an outstanding contribution to improving the model’s prediction performance.

In order to demonstrate that the MLP classifier has better prediction performance, we also investigated the performance of other machine learning methods, including SVM, Random forest, and KNN. These methods were trained and tested using the same fused features. The results are recorded in Figure 2. As shown in Figure 2, the performance of MLP classifiers was better than other classifiers. Therefore, we considered using a MLP-based model to establish a software package.

Many models have been proposed for thermophilic protein identification (Gromiha and Suresh, 2008; Lin and Chen, 2011; Wang et al., 2011; Nakariyakul et al., 2012; Fan et al., 2016; Tang et al., 2017; Li J. et al., 2019; Feng et al., 2020; Guo et al., 2020). All proposed models were established based on machine learning methods and were evaluated by cross-validation. However, our model was examined on independent data. Moreover, the benchmark dataset used in the present study was rigorous and objective. Moreover, most of these published works did not establish available tools that are not only non-practical but also prevent us from making a fair comparison. The only available web-server for the identification of thermophilic proteins was established by Lin and Chen (2011). We performed a comparison with the web server using the same validation dataset. Their model (Lin and Chen, 2011) displayed 95.30% accuracy, while our model produced an accuracy of 96.26%.

Our model produces good prediction performance and shows that the characteristics used can effectively characterize thermophilic proteins. Thus, we performed an analysis on features based on their contribution to model performance. In order to find feature contribution, we used permutation feature importance. The contribution of features to the performance of the model is represented in Supplementary Table 1. The following section will analyze the feature of each descriptor briefly.

The composition and arrangement of amino acids determine the unique function of a protein. At present, the research on thermophilic proteins uses the composition characteristics of amino acids. The current study involves a detailed analysis of AAC. We found that the frequencies of alanine (A), lysine (K), valine (V), isoleucine (I), glutamine (Q), aspartic acid (D), tyrosine (Y), serine (S), glutamic acid (E), and threonine (T) were significantly different between the two classes. It is speculated that these amino acids have crucial information in providing either thermophilicity or non-thermophilicity to proteins. Tyrosine contributed the most to model performance among these amino acids and showed the weight 0.0249 ± 0.0080. Moreover, lysine, glutamic acid, glutamine, and aspartic acid also contributed well to model performance and showed the weights 0.0033 ± 0.0026, 0.0041 ± 0.0017, 0.0036 ± 0.0049, and 0.0162 ± 0.0036, respectively. Glutamate, lysine, tyrosine, glutamic acid, and aspartic acid residues were more common in thermophilic proteins than non-thermophilic proteins. Thermophilic proteins contain highly charged amino acids, which contribute to the thermal stability of proteins. Lysine, glutamine, aspartic acid, and glutamic acid residues belong to charged amino acids, while tyrosine belongs to polar amino acids. These amino acids participate in forming salt bridges and hydrogen bonds, which provide thermal stability to proteins. These results are consistent with previous studies (Liu et al., 2011; Wang et al., 2011; Panja et al., 2020).

Valine and isoleucine showed good ability for thermophilic protein identification. In permutation feature importance, valine and isoleucine showed the weights 0.0201 ± 0.0056 and 0.0101 ± 0.0028, respectively. Isoleucine and valine are hydrophobic amino acids. It has been reported that hydrophobicity contributes to the thermal stability of proteins, as during protein folding, hydrophobic amino acids get buried inside the protein to form a hydrophobic core; this hydrophobic core contributes to the thermal stability of proteins (Baldwin, 2007; Gromiha and Suresh, 2008).

Amino acid alanine, threonine, and serine indicated an important role in model performance and showed the weights 0.0087 ± 0.0018, 0.0122 ± 0.0045, and 0.0031 ± 0.0039, respectively. Figure 3 illustrates the contribution of AAC features to model performance. Non-thermophilic proteins contain more alanine, threonine, and serine residues than thermophilic proteins, consistent with a previous study by Cambillau and Claverie (2000). Alanine carries less charge, while threonine and serine are neutral amino acids, so these amino acids are rarely involved in forming hydrogen bonds and salt bridges, indicating that the proteins enriched with these amino acids can be prone to thermal denaturation (Lin and Chen, 2011).

Amino acid composition is an excellent descriptor to discriminate thermophilic proteins from non-thermophilic proteins. Previous studies have also confirmed the contribution of AAC to protein classification tasks (Gromiha and Suresh, 2008; Mahmoudi et al., 2016). Although AAC plays a good role in protein classification, it also lacks sequence information. The traditional tPseAAC and aPseAAC (Chou, 2001, 2005) are good options for the lack of sequence information in AAC. Wang et al. (2011) and Chen et al. (2016) also confirmed the critical role of these descriptors in protein classification.

Both tPseAAC and aPseAAC are used to describe the sequence information of amino acid residues in protein sequence. In tPseAAC, Xc1.K, Xc1.E, Xc1.D, Xc1.Q, Xc1.T, Xc1.A, Xc1.G, and Xc1.S were valuable features with the weights of 0.0045 ± 0.0019, 0.0041 ± 0.0017, 0.0040 ± 0.0029, 0.0036 ± 0.0049, 0.0035 ± 0.0015, 0.0030 ± 0.0029, 0.0026 ± 0.0027, and 0.0018 ± 0.0026, respectively (Figure 3). The features Pc1.Q, Pc1.E, Pc1.I, Pc1.T, Pc1.A, Pc1.S, Pc1.G, and Pc1.K in aPseAAC presented important contribution to model performance. They showed the weights 0.0086 ± 0.0030, 0.0081 ± 0.0046, 0.0060 ± 0.0057, 0.0047 ± 0.0030, 0.0047 ± 0.0026, 0.0042 ± 0.0030, 0.0034 ± 0.0023, and 0.0018 ± 0.0014, respectively (Figure 3). Our in-depth analysis showed that hydrophobic amino acid and polar amino acid based features were more frequent in thermophilic protein, while uncharged and neutral amino acid based features were more frequent in non-thermophilic proteins.

Dipeptides are also an important feature to distinguish thermophilic proteins from non-thermophilic proteins. Our statistical analysis showed that the occurrence frequencies of KE, LK, EE, EK, AA, LA, KI, IK, KK, and EI have a considerable variance between the two classes of proteins. The ranking of features also confirmed the role of these dipeptides in model performance. Dipeptide KE, LK, EE, EK, AA, LA, KI, IK, KK, and EI showed the weights 0.0113 ± 0.0029, 0.0086 ± 0.0045, 0.0072 ± 0.0019, 0.0070 ± 0.0040, 0.0069 ± 0.0043, 0.0058 ± 0.0010, 0.0043 ± 0.0013, 0.0040 ± 0.0017, 0.0029 ± 0.0026, and 0.0023 ± 0.0017, respectively (Figure 3). Dipeptide KE, LK, EE, EK, KI, IK, KK, and EI have charged at biological pH, showing a great trend of forming salt bridges and hydrogen bonds, which contributes to the thermal stability of proteins. AA and LA have poor charge capability and were found more in non-thermophilic proteins (Nakariyakul et al., 2012; Panja et al., 2020). Previous studies have also confirmed the role of dipeptide composition in identifying thermophilic proteins (Gromiha et al., 2005; Lin and Chen, 2011). MLP model trained on these selected features reveals that these features have good capability to distinguish thermophilic proteins.

The dipeptide deviation from the expected mean also showed meaningful information for the identification of thermophilic proteins. Features including EE, AA, and KE deviation from expected mean showed good ability to identify thermophilic proteins (Table 2). The dipeptide deviation for EE, AA, and KE showed the weights 0.0098 ± 0.0032, 0.0039 ± 0.0028, and 0.0025 ± 0.0012, respectively (Figure 3). Previous studies have also reported the effective contribution of dipeptide deviating from the expected mean in protein classification tasks (Saravanan and Gautham, 2015; Ho Thanh Lam et al., 2020). In addition to these dipeptide-related descriptors, we also considered the composition of k-spaced amino acid pairs, representing the paired amino acid frequency separated by any other amino acid. It is a valuable descriptor and has been widely used in previous studies for protein classification (Jang et al., 2020; Ju and Wang, 2020; Zhang L. et al., 2020; Zulfiqar et al., 2021b). In the present study, E^∗∗K, E^∗∗∗K, A^∗∗A, and A^∗A were found to be containing meaningful information for thermophilic protein identification and showed the weight 0.0077 ± 0.0041, 0.0057 ± 0.0014, 0.0031 ± 0.0034, and 0.0028 ± 0.0015, respectively (Figure 3).

Composition, transition, and distribution involves the composition, transition, and distribution of hydrophobic, polar, and neutral residues. Like other descriptors, the hydrophobic and polar residue-based features of CTD were more frequent in thermophilic proteins while neutral residues-based features were more frequent in non-thermophilic proteins. Permutation feature importance of descriptor CTD is represented in Supplementary Table 2. In previous studies, the CTD has been extensively used for protein classification purposes. Wang et al. (2011) and Zulfiqar et al. (2021b) also reported CTD as a valuable descriptor for thermophilic protein identification. In the present study, the CTD showed an excellent capability to identify thermophilic proteins (Table 2). For CTD, all features performed better than the selected features, so we used CTD features without selection. MLP model trained on CTD features performed good results (Table 2).