MDDeep-Ace: species-specific acetylation site prediction based on multi-domain adaptation

Dataset

We collected ten species of lysine acetylation sites from PLMD dataset (Xu et al., 2017), including 6,078 H. sapiens proteins, 3,645 M. musculus proteins, 2,960 S. cerevisiae proteins, 4,359 R. norvegicus proteins, 1,251 S. japonicum proteins, 231 A. thaliana proteins, 1,860 E. coli proteins, 1,146 B. velezensis proteins, 1,214 P. falciparum proteins, 336 O. sativa proteins. We utilized the CD-HIT tool (Huang et al., 2010) to cluster protein sequences for each species, filtering out homologous proteins exceeding a 40% similarity threshold. Negative samples were defined as lysine residues that were not experimentally verified as acetylation sites from the non-homologous proteins. Subsequently, 10% of both positive and negative samples per species were randomly allocated to an independent test dataset, with the remaining data used for training and validation. To address sample imbalance and prevent over-optimization, we balanced the dataset by equalizing the number of positive and negative samples (Table 1). Finally, protein sequences were encoded into numerical vectors using one-hot encoding, where each amino acid is represented by a vector with a single ‘1’ and all other elements as ‘0’. For example, for lysine (K), the one-hot encoding is “000000100000000000000”, while for serine (S), it is “000000000000000100000”. By following previous studies (Jia et al., 2016; Liu et al., 2021), symmetrical 31-residue windows, spanning from −15 to +15 with lysine at the center, were extracted from protein fragments to serve as training and testing samples. Meanwhile, we also conducted comparative experiments on different window sizes, demonstrating the rationality of choosing a length of 31. The relevant settings and results of the experiments are listed in the Table S1. To standardize fragment lengths, we used the placeholder amino acid ‘X’ for padding missing residues or replacing non-standard amino acids. The fragments were transformed into numerical vectors using one-hot encoding. Each fragment was encoded as a 31 × 21 feature, where 31 represents the sequence length and 21 the amino acid categories.

MDDeep-Ace architecture

MDDeep-Ace is a novel multi-domain adaptation learning based species-specific lysine acetylation site prediction method, with its workflow illustrated in Fig. 1. We curated data from ten species, using nine as source domains when training a predictive model for a specific target species.

The training dataset consists of two components: N labeled source dataset D^S = {D^S₁, D^S₂, …, D^S_N}. The categories across all source domains are identical. $D^{S_j}={\left\{x_i^{S_j},y_i^{S_j}\right\}}_{i=1}^{N_{S_j}}$, where x denotes the one-hot encoding sequence feature of sample and y denotes the label of sample, $x_i^{S_j}$ and $y_i^{S_j}$ denote sequence feature of jth source domain and corresponding label, N_Sj denotes the number of sample of jth source domain; a labeled target dataset $D^T={\left\{x_i^T,y_i^T\right\}}_{i=1}^{N_T}$, where $x_i^T$ and $y_i^T$ denote sequence feature and label of target domain, N_T denotes the number of labeled target domain.

MDDeep-Ace employs a CNN-LSTM hybrid network g to extract sequence features and map protein sequences into an embedding space, followed by a classifier h for feature classification. The CNN layer, with 21 input channels and 128 output channels, uses a kernel size of 3 with a padding of 1 to capture local sequence features. The subsequent LSTM layer, with a hidden size of 128, processes the 128-channel CNN output to capture long-range dependencies. An average pooling layer (kernel size 2, stride 2) reduces dimensionality, enabling efficient extraction of both local and long-range features. The classifier applies a softmax operation to generate predictions. (1)${\displaystyle p\left(y|x\right)=h\left(g\left(x:\theta \right);W\right).}$

The entire framework employs CNN-LSTM hybrid network as shared feature extractors across multiple domains. This shared structure allows data from different domains to be mapped to the same feature space, ensuring the extraction of domain-invariant features common to all fields. The domain alignment component minimizes feature distribution discrepancies between source and target domains through a domain discrepancy loss.

We employ Kullback–Leibler (KL) divergence to measure feature distribution distances between domains, calculating a dynamic adjustment factor to prioritize source domain samples closer to the target domain while suppressing less similar ones. By differentially reducing inter-domain discrepancies, the framework better learns domain-invariant feature representations, thereby achieving more accurate recognition of target domain samples. To support the domain adaptation, the total loss comprises two parts: the classification loss of all source domains and the inter-domain discrepancy loss. By adjusting the weight proportions of each source domain within the total loss function, the influence of each source domain on the target domain is determined by the specific weight allocation. The total loss function is shown below: (2)${\displaystyle L_{total}=L_C^T\left(g,h\right)+{\alpha }_j\left(L_C^{S_j}\left(g,h\right)+L_{KL}^{S_j}\left(g\right)\right)}$ where $L_C^T\left(\cdot \right)$ and $L_C^{S_j}\left(\cdot \right)$ represent the classification loss function applied to the target domain and the jth source domain, respectively, $L_{KL}^{S_j}\left(\cdot \right)$ denotes domain difference loss function for the jth source domain whileα_j is the weight coefficient of the jth source domain in the loss function. The classification loss function is defined as follows: (3)${\displaystyle L_c\left(g,h\right)=-\frac{1}{N}\sum _{i=1}^Q{y}_{i}lnp\left(y_i|x_i\right)+\left(1-y_i\right)ln\left(1-p\left(y_i|x_i\right)\right),}$ where Q is the number of samples. The domain difference loss function is defined as follows: (4)${\displaystyle L_{KL}^{S_j}=\mathrm{E}\left(\sum _{k=1}^{M}g{\left(x^T\right)}_k\cdot log\frac{g{\left(x^T\right)}_k}{g{\left(x^{S_j}\right)}_k}\right)}$ where E(⋅) denotes the mathematical expectation, g(x)_k denotes the kth dimension feature of g(x). The weight coefficient is computed by: (5)${\displaystyle {\alpha }_j=\frac{exp\left({\mathrm{D}}_j\right)}{\sum _{c=1}^{N}exp\left({\mathrm{D}}_c\right)}}$ where D_j denotes the KL divergence between the features of jth source domain and the target domain. Supplementary File S1 details the corresponding code for loss function calculations.

Performance evaluation indicators

To evaluate MDDeep-Ace, we used five performance metrics: area under the receiver operating characteristic (ROC) curves (AUC), sensitivity (Sn), accuracy (Acc), precision (Pre) and F1 score. The formula of these metrics are as follows. (6)${\displaystyle Acc=\frac{TP+TN}{TP+TN+FP+FN}}$ (7)${\displaystyle Sn=\frac{TP}{TP+FN}}$ (8)${\displaystyle Sp=\frac{TN}{TN+FP}}$ (9)${\displaystyle Pre=\frac{TP}{TP+FP}}$ (10)${\displaystyle F1=\frac{2\times Pre\times Sn}{Pre+Sn}}$

where TP, TN, FP, FN represent true positives, true negatives, false positives and false negatives, respectively.

Sequence analysis

We analyzed the patterns of lysine acetylation using two-sample logo tool (Fig. 2). The analysis revealed distinct sequence differences across species. For example, Glycine/G residue is found enriched in position −1, −4 and +2 for H. sapiens, but not in O. sativa and A. thaliana. However, compared with R. norvegicus and H. sapiens, phenylalanine/F and tyrosine/Y residue prefer to enrich in position −1, −2 and +2 on B. velezensis. These findings underscore the need for species-specific prediction models. However, there still exist similarities between different species. For example, in most species, lysine (K) shows a significant depletion at diverse positions surrounding the acetylation site, both upstream and downstream. Meanwhile, the asparagine/D is enriched at −1 and +1 position on several species. Additionally, we also found that the degree of similarity varies among different species. For example, the types of amino acids enriched and depleted near acetylation sites in S. cerevisiae and E. coli are almost identical, whereas there are significant differences in the distribution of amino acids near acetylation sites between S. cerevisiae and B. velezensis. These observations highlight the potential of domain adaptation methods to leverage inter-species similarities, enhancing species-specific acetylation site prediction.

Comparison with existing methods

We compare MDDeep-Ace with several existing acetylation site prediction methods including DeepDA-Ace (Liu, Wang & Xi, 2022), CapsNet (Wang, Liang & Xu, 2019), and PAIL (Deng et al., 2016) based on test dataset, and the construction approach of the test data set is shown as 2.1. In order to visualize the predictive performance, we show the ROC curves of different methods on multiple species in Fig. 3 and Fig. S1. From the Figure, we find that MDDeep-Ace, DeepDA-Ace and CapsNet achieved better performance than PAIL on all species. Take S. cerevisiae as an example, the AUC value of MDDeep-Ace, DeepDA-Ace and CapsNet are 0.808, 0.780 and 0.736 respectively, while PAIL only obtains 0.541 AUC value, demonstrating the superiority of deep learning-based methods. Besides, MDDeep-Ace and DeepDA-Ace outperform the other methods, suggesting that species-specific models have more advantages than general models. Meanwhile, it indicates that the transfer learning technique has an advantage over other deep learning methods in species-specific acetylation site prediction. For instance, DeepDA-Ace achieves 3.1%, 4.9%, 4.4% and 4.9% improvement compared to CapsNet on S. japonicum, R. norvegicus, S. cerevisiae and M. musculus. Additionally, our proposed MDDeep-Ace combines the benefits of transfer learning and multi-domain adaptation technique, and obtains the highest AUC value among existing methods on almost all species in test dataset. Taking B. velezensis as an example, compared to DeepDA-Ace, CapsNet and PAIL, the AUC of MDDeep-Ace has improved by 5.3%, 7.7% and 29.4%, respectively. As for R. norvegicus, MDDeep-Ace obtains 0.764 AUC value, while the corresponding AUC value of DeepDA-Ace, CapsNet and PAIL are 0.732, 0.683 and 0.538, respectively.

Additionally, to evaluate the performance of MDDeep-Ace more comprehensively, these four metrics are also considered by us: sensitivity (Sn), accuracy (Acc), precision (Pre) and F1 score at high specificity (Sp). We compared these metrics of the four methods at a specificity of 0.900, and the results are illustrated in Fig. 4 and Fig. S2. Clearly, MDDeep-Ace demonstrates superior performance across most metrics for all species. For example, on S. cerevisiae, MDDeep-Ace shows notable improvements: accuracy (Acc) increases by 5%, sensitivity (Sn) by 9.8%, precision (Pre) by 3.6%, and F1 score by 9.3% compared to PAIL. When compared to DeepDA-Ace, the improvements are also significant, with Acc, Sn, Pre, and F1 score increasing by 10.7%, 21%, 10.3%, and 21.8%, respectively. Moreover, compared with CapsNet, the gains are even more pronounced, with Sn, Acc, Pre, and F1 score increasing by 34.9%, 17.8%, 28.5%, and 40.4%, respectively. MDDeep-Ace not only outperforms these models in S.cerevisiae but also achieves competitive results across other species.

The superior performance of MDDeep-Ace can be attributed to two main factors. Firstly, MDDeep-Ace leverages multi-domain adaptation, effectively transferring acetylation knowledge from multiple species to enhance the model’s predictive power. Secondly, MDDeep-Ace employs a CNN-LSTM hybrid network, which captures both local features and long-range dependencies, further boosting the model’s predictive performance. Overall, MDDeep-Ace not only surpasses existing methods with respect to accuracy, sensitivity, precision, and F1 score, but also represents a significant advancement in species-specific lysine acetylation site prediction, offering a robust and versatile tool for bioinformatics research.

Ablation study

To validate the efficacy of the introduced multi-domain adaptation method, we compared MDDeep-Ace with three baseline models. (1) Species-specific direct training (SSDT) model: for each species, a supervised model is trained using only the labeled data specific to that species. (2) Species-specific domain adaptation (SSDA) model: this model leverages domain adaptation techniques, but uses only human data as source domain. It trains species-specific prediction models based on labeled data from both human and the target species. (3) Uniform-weight MDDeep-Ace model: to assess the impact of the domain weighting strategy, this model uses the same architecture as MDDeep-Ace but applies equal weighting to all source domains. The above three methods are all based on the CNN-LSTM hybrid network. The differences among them lie in either the distinct training data they utilize or the different training strategies they adopt. The results clearly demonstrate that the adaptative weighting mechanism in MDDeep-Ace plays a critical role in enhancing prediction accuracy.

Table 1 illustrates the performance of MDDeep-Ace compared to baseline methods, demonstrating that MDDeep-Ace generally outperforms the uniform-weight MDDeep-Ace in predicting PTM sites across various species. For instance, MDDeep-Ace achieves the highest AUC values across all species, reflecting its superior capability. However, it is also notable that the uniform-weight MDDeep-Ace performs comparably to the SSDA method in some species. For example, the AUC values of the uniform-weight MDDeep-Ace are only 0.5% and 0.7% higher than those of the SSDA method in M. musculus and R. norvegicus, respectively. This can be attributed to the fact that human acetylation sites constitute a significant portion of the training data for the uniform-weight MDDeep-Ace, and the acetylation patterns of M. musculus and R. norvegicus are quite similar to those of humans. Conversely, the uniform-weight MDDeep-Ace performs poorly on O. sativa and A. thaliana, likely due to the considerable differences in acetylation patterns between these species and humans, making it challenging for the uniform-weight MDDeep-Ace to effectively capture species-specific patterns. These findings highlight the limitations of the uniform-weight MDDeep-Ace in accurately predicting PTM sites across all species, underscoring the necessity of establishing adaptative weight species-specific prediction models for enhanced performance.

Although SSDT is a species-specific method, it exhibits poor performance in some species such as A. thaliana, O. sativa, and B. velezensis, indicating the difficulty of building effective species-specific prediction models through direct training when training data is insufficient. To address this issue, domain adaptation technology, as employed in the SSDA method, transfers knowledge from extensive human PTM data to improve the prediction accuracy for other species. For instance, the SSDA method achieves AUC values of 0.793, 0.804, and 0.805 in A. thaliana, O. sativa, and B. velezensis, representing improvements of 8%, 11.5%, and 6.4% over SSDT. Moreover, with the integration of multiple species data, MDDeep-Ace further enhances prediction accuracy, achieving AUC values of 0.813, 0.823, and 0.847 in these species, respectively. In addition to these species, MDDeep-Ace also achieves the optimal performance across all other species, confirming its effectiveness in species-specific acetylation site prediction. The multi-domain adaptation strategy, combined with the incorporation of diverse species data, allows MDDeep-Ace to outperform existing methods.

Besides, we also compared the F1 scores, accuracy (Acc), sensitivity (Sn), and precision (Pre) for these methods, and the results are listed in Table 2. From the table, it is evident that the methods utilizing multi-domain adaptation technology achieved higher performance across several indicators. For example, in R. norvegicus, the Pre values of MDDeep-Ace and SSDA methods are 0.783 and 0.752, respectively, with MDDeep-Ace showing a 3.1% improvement over SSDA. Additionally, the Acc values of MDDeep-Ace and SSDA methods are 0.620 and 0.592, representing a 2.8% increase for MDDeep-Ace over SSDA. Moreover, compared to the uniform-weight MDDeep-Ace model, MDDeep-Ace demonstrates notable improvements. For instance, in S. cerevisiae, the Sn, Acc, Pre, and F1 score of MDDeep-Ace are 0.463, 0.677, 0.827, and 0.593, respectively, while the corresponding values for the uniform-weight MDDeep-Ace model are 0.431, 0.661, 0.817, and 0.565. This reflects a 3.2% increase in Sn, a 1.6% increase in Acc, a 1% increase in Pre, and a 2.8% improvement in F1 score. Among the domain adaptation-based methods, our proposed MDDeep-Ace shows clear advantages across most species. For example, in B. velezensis, the Pre, Acc, F1, and Sn of MDDeep-Ace are 0.846, 0.694, 0.634, and 0.507, respectively, whereas the SSDA method only achieves 0.774, 0.594, 0.449, and 0.316. This represents improvements of 7.2%, 10%, 18.5%, and 19.1% in Pre, Acc, F1, and Sn, respectively, when using MDDeep-Ace. These results further support the conclusions drawn from the comparison using AUC values.