We utilized the dataset developed in the iSulf-Cys framework (Xu et al., 2016), which is based on the work of Yang et al. (2014). The dataset comprises 7,124 non-sulfenylation cysteine sites and 1,045 sulfenylation cysteine sites. In total, the data are derived from 778 protein sequences and include 1,105 sulfenylation cysteine sites. The corresponding protein sequences are available at https://ndownloader.figstatic.com/files/5004001.

For each cysteine site subjected to analysis, Xu et al. constructed peptide samples by extracting the 10 amino acids upstream and downstream, resulting in a sequence window of size 21 centered on the target cysteine. If the surrounding sequence contained fewer than 10 residues, the placeholder residue X was used to pad the sequence. The representation of a sample is as Equation 1: Peptide = A − 10 A − 9 … C A 1 A 2 A 3 … A 10 (1)

where A i ∈ { 20  amino acids } ∪ { X } , and i range from −10 to 10, representing the one-dimensional coordinate relative to the central cysteine. A peptide segment with a central cysteine that is a sulfenylation site is defined as a positive sample, while all others are classified as negative samples. To minimize redundancy and reduce homology bias, peptide segments sharing more than 40% sequence similarity were removed from the dataset. Ultimately, 145 positive and 268 negative samples were randomly selected as the test set, while the remaining 900 positive and 6,858 negative samples were allocated to the training set. A summary of the dataset is presented in Table 1. This dataset has been widely adopted in sulfenylation site prediction studies (Ju and Wang, 2018; Khan and Pi, 2020; Do et al., 2020), and it exhibits a pronounced imbalance between positive and negative classes.

In this study, two strategies were considered for processing peptide segment features. First, bioinformatics-based features were introduced by encoding sequences using the Position-Specific Scoring Matrix (PSSM). Second, natural language processing techniques were applied to generate word embeddings from protein sequences, thereby enabling the model to fully capture the latent semantic information embedded within the sequences. This approach has been validated as effective in previous study (Do et al., 2020).

Language is one of the most intuitive forms of representation, and humans have long sought to describe the natural world through linguistic expressions. By assigning specific letters to amino acids, protein sequences can be directly encoded as character strings. These sequences contain rich semantic content that can be mined using appropriate natural language processing (NLP) techniques. This section focuses on local protein-related tasks, particularly on predicting residue-level features from sequence data. Window-based models have proven effective in various PTM site prediction tasks. For instance, Brandes et al. (2015) introduced a framework that predicts PTM sites by segmenting protein sequences into fixed-size windows, highlighting the potential of semantic features in protein feature engineering.

With the continuous advancement of NLP technologies (Vaswani et al., 2017; Devlin et al., 2018), Word2Vec has emerged as a major breakthrough in the field (Mikolov et al., 2013a). By leveraging contextual information, it generates continuous low-dimensional vector representations for words (or characters), capturing latent semantic relationships embedded in language. This advancement has inspired novel approaches to protein sequence analysis. FastText, an enhanced version of Word2Vec, further expands upon this idea (Joulin et al., 2016). Compared to Word2Vec, fastText introduces subword modeling, which effectively captures internal morphological and structural features of words. This is particularly suitable for protein sequence modeling, as such sequences frequently contain variable and localized motifs. FastText not only captures local n-gram-like information during training but also addresses the Out-of-Vocabulary (OOV) problem—a critical advantage when analyzing complex biological sequences.

The core idea of fastText is to enhance traditional word embeddings by incorporating subword modeling. In contrast to the Word2Vec model, which assigns a distinct vector to each word, fastText represents a word as a combination of its constituent subwords, thereby capturing subtle internal variations within words. For example, the word protein can be split into subwords such as pro, ote, and ein, all of which participate in the model’s training. By encoding and embedding these subwords, fastText learns more nuanced semantic representations, making it especially effective in handling diverse and irregular protein sequences.

From a mathematical standpoint, fastText is built upon the Skip-gram architecture, where the objective is to maximize the conditional probability of predicting a target word given its surrounding context (Mikolov et al., 2013b). Given a training corpus composed of word sequences w 1 , w 2 , … , w T , the objective function of fastText is to learn embeddings for both words and their subwords by minimizing the following loss function (Equation 3): L = − ∑ t = 1 T ∑ − c ≤ j ≤ c , j ≠ 0 log ⁡ p w t + j | w t (3)

where w t is the current word, w t + j is the context word, c is the size of the context window, and p ( w t + j | w t ) is the conditional probability of predicting the target word w t + j . In fastText, p ( w t + j | w t ) is no longer a simple word’s conditional probability but a conditional probability of the subwords. Specifically, for each word w , it is split into several subwords s 1 , s 2 , … , s k , and the probability of the target word is calculated as Equation 4: p w t + j | w t = ∏ s ∈ S w t + j p s | w t (4)

where S ( w ) represents the set of subwords for the word w , and p ( s ∣ w t ) denotes the conditional probability of subword s given the context word w t .

In practical applications, fastText represents the semantic content of a word through both its own embedding and the embeddings of its constituent subwords. For protein sequences, dipeptides or tripeptides derived from amino acid residues can be regarded as subwords. Using this representation, fastText can convert short fragments of protein sequences—such as amino acid pairs or triplets—into vector embeddings, capturing the underlying contextual relationships between substructures. In this study, each peptide sequence is decomposed into 3-g as subwords. For a sequence consisting of n amino acids, this yields n − 2 trigrams, and hence, generates n − 2 corresponding word vectors that embed the sequence’s implicit order and semantic content. This process constitutes the fastText-based protein sequence embedding used in this study.

As described in Section 2.1, the dataset used in our study exhibits a pronounced imbalance between positive and negative samples, necessitating the application of resampling techniques. Common resampling strategies include down-sampling and up-sampling. Down-sampling involves reducing the number of majority class instances to achieve a balanced class distribution, whereas up-sampling increases the number of minority class samples to the same effect. To minimize information loss from the sample space and retain protein residue-level features, up-sampling is prioritized. Synthetic Minority Over-sampling Technique (SMOTE) is one of the most widely used up-sampling methods (Chawla et al., 2002). Numerous variants of SMOTE have been proposed in subsequent studies, and this study adopts Support Vector Machine Synthetic Minority Over-sampling Technique (SVMSMOTE) for up-sampling (Demidova and Klyueva, 2017). The resampling procedure is outlined as follows.

1. Train the SVM model: First, a support vector machine (SVM) classifier is trained using minority class samples to determine the decision boundary. By maximizing the margin between classes, SVM effectively learns to delineate the optimal separation boundary for the training data.

where x new is the generated synthetic sample, x sv is the current support vector, x neigh is its nearest minority class neighbor, and λ is a random number in the range [0, 1] that determines the interpolation proportion.

4. Generate the required number of synthetic samples: Based on the method above, SVMSMOTE generates enough synthetic samples to achieve a balanced class distribution in the training dataset.

For a schematic illustration of the SMOTE and SVMSMOTE procedures, please refer to Supplementary Figure S1.

Considering the sequential nature of the two feature spaces used for protein representation in this study, a Bi-LSTM-based approach was adopted. The strength of LSTM in handling protein sequences lies in its capability to effectively capture long-range dependencies. In protein prediction tasks, especially in the identification of post-translational modification (PTM) sites, long-range contextual dependencies are critical (Shi et al., 2025; Vaswani et al., 2017) as relevant information in protein sequences often spans beyond local neighborhoods. Accurate functional prediction therefore, relies on a model’s ability to extract such extended context. Traditional recurrent neural networks (RNNs) struggle in these scenarios due to challenges such as vanishing and exploding gradients, which hinder learning across long temporal sequences. In contrast, Bi-LSTM incorporates both forward and backward contextual information along the sequence, allowing it to capture bidirectional dependencies that are essential for accurate PTM site prediction.

A standard Bi-LSTM consists of a forward LSTM and a backward LSTM. Specifically, the LSTM architecture introduces memory cells that determine whether to retain or discard information from past inputs. These memory cells are regulated by three gates: the input gate, forget gate, and output gate. The input gate controls how much of the current input is stored in the memory cell, the forget gate decides which parts of the previous memory to discard, and the output gate determines the information to be passed to the next layer in the network.

In this study, Bi-LSTM is employed as the core model to process protein sequence data. By combining the forward and backward passes, Bi-LSTM effectively captures bidirectional contextual information within the protein sequence, thereby enhancing the accuracy of sequence-based predictions.

Let f denote the forward LSTM process and b denote the backward LSTM process. For each time step t , the forward LSTM computes the following Equations 6‐10: i t = σ W i ⋅ h t − 1 ( f ) , x t + b i (6) f t = σ W f ⋅ h t − 1 ( f ) , x t + b f (7) o t = σ W o ⋅ h t − 1 ( f ) , x t + b o (8) c t ( f ) = f t ⋅ c t − 1 ( f ) + i t ⋅ tanh W c ⋅ h t − 1 ( f ) , x t + b c (9) h t ( f ) = o t ⋅ tanh c t ( f ) (10)

Where σ is the sigmoid activation function, tanh is the hyperbolic tangent function, x t is the current input, h t − 1 ( f ) is the previous time step’s hidden state, and c t − 1 ( f ) is the previous time step’s memory state.

The backward LSTM operates similarly to the forward LSTM but processes the sequence in reverse order. In the standard Bi-LSTM model, the outputs of both the forward and backward LSTMs are combined to calculate the final output for each time step. However, in this study, the initial hidden state of the backward LSTM is modified. Specifically, the initial hidden state of the backward LSTM is set as a weighted sum of the forward LSTM’s initial hidden state, with the weights being trainable. This modification is expressed as Equation 11: h 0 ( b ) = σ W i h i ( f ) (11)

Where w i is the weight for the i -th step of the forward LSTM, and h i is the hidden state at the i -th step. This improvement enables the backward pass to integrate information learned from the forward pass, enhancing backward propagation.

Additionally, the model incorporates residual connections, which can be expressed as Equation 12: h t + 1 ( f ) = α ⋅ h t + 1 ( f ) + 1 − α ⋅ h t ( f ) (12)

Where α is a hyperparameter that adjusts the weighting between the current hidden state and the previous hidden state. The residual connection retains information from the previous time step’s hidden state, enhancing the model’s ability to capture long-range dependencies.

Feature extraction in this study is based on both bioinformatics features and semantic features derived from protein sequence fragments. These two types of features emphasize different aspects of information representation: bioinformatics features highlight the structural and functional features of the sequence, whereas semantic features capture contextual relationships at a natural language level. However, due to potential redundancy or inconsistency between these two modalities, an Adaptive Feature Fusion Module is proposed to integrate and utilize them more effectively. This module is designed with two primary objectives: first, to accurately identify and eliminate redundant information to prevent irrelevant or noisy features from impairing model performance; and second, to enhance discriminative features that are crucial to the target task, thereby improving both the learning efficiency and predictive accuracy of the model. The Adaptive Feature Fusion Module not only manages the allocation of importance across different feature types, but also dynamically adjusts the fusion strategy to accommodate variations across different tasks or data distributions.

This study proposed a fusion module leveraging both self-attention and cross-attention mechanisms. The goal of this module is to dynamically assign feature weights using these attention mechanisms, thereby extracting salient information and achieving efficient fusion. A detailed introduction to this module is provided below.

Let the features to be fused be H P (size N × d 1 ) and H W (size N × d 2 ). Since the dimensions of H P and H W may differ, projection layers were initially established to map both features into a common dimensionality d . This is done via linear projection followed by a non-linear activation function σ (Equations 13, 14): H 1 P = σ W 1 H P (13) H 1 W = σ W 2 H W (14) where W 1 ∈ R d × d 1 and W 2 ∈ R d × d 2 .

Next, self-attention is applied to both H 1 P and H 1 W . Self-attention is a mechanism used to compute the relevance of each element in the sequence to every other element, allowing it to capture global context. It dynamically assigns higher weights to important elements in the sequence, thus improving the model’s ability to understand complex dependencies. The calculation of Self-attention is as Equation 15: Self - Attention Q , K , V = softmax Q K T d V (15)

where Q , K , and V are linear transformations of H 1 P and H 1 W . Specifically, self-attention is applied to extract their respective key feature patterns (Equations 16, 17): H 2 P = Self - Attention H 1 P , (16) H 2 W = Self - Attention H 1 W , (17)

Next, H 2 P and H 2 W are combined using the cross-attention mechanism to model the interaction between the two types of features. Cross-attention extracts important information by focusing on the most relevant parts of one feature set relative to another. By doing so, it effectively fuses features from different sources and highlights key associations between them. The computation of Cross-attention is similar to Self-attention, but the design of Q , K , and V differs (Equation 18). Cross - Attention Q , K , V = softmax Q K T d V (18)

Here, Q comes from the linear transformation of one feature set, while K and V come from the linear transformations of the other feature set. Therefore, relevant semantic feature information can be extracted from the bioinformatics features, as expressed by Equations 19, 20: H P W = Cross - Attention H 2 P , H 2 W , (19) H W P = Cross - Attention H 2 W , H 2 P , (20)

Finally, we perform a weighted fusion of the outputs from both Self-attention and Cross-attention (Equation 21): H F = α H P W + H 2 W + 1 − α H W P + H 2 P (21)

where α is a weight parameter. The fused feature representation, H F , is then passed to the subsequent dense layers.

Because positional biases around modified cysteines are supported by redox chemistry and chemoproteomics (Gupta and Carroll, 2014; Huang et al., 2019), the attention weights in the Adaptive Feature Fusion module serve as a proxy for residue-level importance within the ± 10 window, encouraging the model to up-weight positions consistent with these biologically plausible contexts.

To evaluate the robustness and generalization capability of our proposed method, we employed a repeated 10-fold cross-validation strategy. Specifically, the entire dataset was randomly partitioned into 10 folds, of which 9 folds were used for training and one fold for testing. Similar to other studies, we adopted commonly used evaluation metrics, including the Area Under the Receiver Operating Characteristic Curve (auROC), Area Under the Precision-Recall Curve (auPRC), Sensitivity (Sn), Specificity (Sp), Accuracy (ACC), and Matthews Correlation Coefficient (MCC), to comprehensively assess the model’s predictive performance. The formulas for the latter four metrics are as Equations 22–25: Sn = T P T P + F N (22) Sp = T N T N + F P (23) ACC = T P + T N T P + T N + F P + F N (24) MCC = T P ⋅ T N − F P ⋅ F N T P + F P T P + F N T N + F P T N + F N (25)

Where T P , T N , F P , and F N represent the true positives, true negatives, false positives, and false negatives, respectively. Since we use an imbalanced dataset, it is highly appropriate to include auROC and auPRC as evaluation metrics. auPRC better reflects the model’s performance on the positive class, especially in situations of class imbalance (Hancock et al., 2023).

Beyond model performance assessment, we quantitatively evaluate feature adequacy using a sequence compression–based approximation of Kolmogorov complexity (KC) (Li and Vitányi, 1997). In plain terms, Kolmogorov complexity asks “what is the shortest description that can reproduce a given sequence?“. If the feature representation preserves most of the regularities and patterns in the raw sequence, then compressing the features (or the features concatenated with the raw sequence) will be nearly as effective as compressing the raw sequence itself. Therefore, a small estimated information-loss rate indicates that the features retain most sequence-level biochemical signals relevant to prediction, whereas a large loss suggests that important patterns may have been discarded. The information-loss rate η loss is defined as Equation 26: η loss = K X ∣ f X K X = K X , f X − K f X K X . (26) Here K ( X ) denotes the (prefix) Kolmogorov complexity of a string X , formally (Equation 27) K X = min {   | p | :  a program  p  outputs  X on a fixed universal  Turing machine  U   } . (27)

In words, K ( X ) is the length (in bits) of the shortest program that generates X on U . This definition is invariant up to an additive constant (independent of X ) due to the invariance theorem (Li and Vitányi, 1997). Likewise, K ( X , f ( X ) ) is the joint complexity of the concatenated description of X and f ( X ) , and K ( X ∣ f ( X ) ) is the conditional complexity (the additional information needed to reconstruct X given f ( X ) ).

In this study, X denotes the raw amino-acid sequence fragment and f ( X ) the feature-extraction mapping (fastText subword embeddings and the PSSM). Since exact Kolmogorov complexity is non-computable (due to the undecidability of the halting problem) (Li and Vitányi, 1997), we approximate it using off-the-shelf lossless compressors, which provide practical upper bounds to K . Concretely, let C ( ⋅ ) be the compressed length (bytes) under a fixed ZIP/DEFLATE compressor; we estimate (Equation 28) K ̂ X = C X , K ̂ f X = C f X , K ̂ X ∣ f X ≈ C f X ‖ X − C f X , (28) where ‖ denotes concatenation with a delimiter to avoid boundary artifacts. This estimator is closely related to the normalized compression distance framework (Cilibrasi and Vitányi, 2005). For biological symbol sequences, compression-based complexity estimates—often based on Lempel-Ziv variants—have a long history and have been successfully applied to genetic texts (Gusev et al., 2000).

This metric is particularly crucial for our sulfenylation site prediction task where we process 21-residue sequence fragments without structural information. Given that our features convert symbolic sequences to numerical representations for BiLSTM classification, quantifying information preservation ensures the extracted features maintain biochemically relevant patterns while reducing noise (Zvonkin and Levin, 2007; Kolmogorov, 1983).

To better illustrate the differences between sulfenylation site samples and non-sulfenylation site samples at the residue level, this study employed the Two-Sample Logo technique. This method utilizes statistical analysis to identify significant distinctions between central cysteine residues and their surrounding amino acids in positive and negative samples (Schneider and Stephens, 1990; Crooks et al., 2004). The Two-Sample Logo is a commonly used visualization tool in sequence content analysis for various PTMs, offering clear insights into amino acid distribution patterns around the modification site and serving as explanatory support for site prediction.

Figure 2 depicts the Two-Sample Logo comparison between the positive and negative sample sets used in this study, with a significance threshold (P-value) set to 0.5. The figure reveals notable differences between multiple positive and negative samples. For instance, non-central cysteine residues are more prevalent in negative samples, particularly at positions +1, +4, +5, +6, +7, +8, +9, +10, +14, +15, +17, and +21. This observation is consistent with the findings of the Bi-directional Gated Recurrent Unit network with Self-Attention (BiGRUD-SA) model, despite differences in datasets, highlighting the strong generalization capability of the BioSemAF-BiLSTM model proposed in this study (Zhang et al., 2023). Moreover, other amino acids also exhibit statistically significant differences between the two sample groups, such as.

• Lysine (K): Frequently appears in positive samples at positions +1, +3, +4, +5, +6, +7, +9, +15, +18, +19, +20;

Importantly, the enrichment of charged residues such as Lys (K) and Glu (E) flanking the central cysteine is consistent with established redox chemistry and chemoproteomic observations: local electrostatics and hydrogen-bonding networks modulate cysteine thiol pKa and the lifetime of sulfenic acid, biasing oxidation toward specific sequence microenvironments (Gupta and Carroll, 2014). Large-scale surveys likewise report redox-sensitive cysteines occurring in characteristic local contexts across proteomes (Huang et al., 2019; Fu et al., 2023; Akter et al., 2018). This concordance provides biological support for the sequence patterns highlighted by our Two-Sample Logo and for subsequent sequence-based prediction.

These features capture significant distinctions between positive and negative samples at the corresponding P-values, thereby confirming the reliability and effectiveness of the sample data and laying a robust foundation for further research.

Additionally, according to Do et al. (2020), no definitive correlation has been established between sulfenylation sites and their surrounding amino acid compositions. However, investigating its biological significance remains crucial. Mu et al. (Mu et al., 2024) further elucidated the mechanisms underlying protein sulfenylation. They highlighted that variations in the generation and distribution of intracellular oxidants (e.g., hydrogen peroxide) lead to differential oxidant levels across organelles and subcellular compartments, thereby influencing sulfenylation occurrence. The pKa of the cysteine thiol group is significantly modulated by protein structure, including its secondary and tertiary structures, hydrogen bonding, macrodipoles, and electrostatic interactions. Lower pKa values enhance the nucleophilicity of the thiol group, rendering it more susceptible to sulfenylation (Kortemme and Creighton, 1995; Ferrer-Sueta et al., 2011). Additionally, the polarization of the thiol group and the non-bonding electrons of adjacent atoms directly modulate its nucleophilicity, thereby affecting sulfenic acid formation.

These studies offer critical biological insights into the sequence-based relationship between sulfenylation sites and their surrounding amino acids. They also provide theoretical support for model development and optimization, further advancing our understanding of sulfenylation mechanisms.

To ensure optimal model performance, we performed an extensive tuning process for the key hyperparameters of the BioSemAF-BiLSTM model. These hyperparameters include n-gram subword length (Length), fastText word vector dimensionality (Dim), and the number of attention heads in the adaptive feature fusion module (Head). For each type of parameter, values were progressively adjusted within reasonable ranges, with appropriate step sizes set to precisely identify local optima. Subsequently, by integrating the optimal settings across multiple parameters, the global best configuration for overall model performance was determined.

Table 2 illustrates the effects of varying parameter configurations on model performance, including n-gram subword lengths (2, 3, 4), fastText embedding dimensions (30, 60, 90), and the number of attention heads in the adaptive feature fusion module (3, 6, 9). The evaluation metrics used in the analysis include ACC, MCC, Sn, and Sp. After individually tuning each hyperparameter and evaluating the overall performance, the optimal settings for these three key parameters were determined as follows.

• n-grams Subword Length: Selected as 3, dividing the protein sequence into subwords with three amino acids as one unit;

Through analyzing the impact of different parameters on model performance during the debugging process, the following conclusions were drawn.

1. The length of n-grams subwords significantly affects the model’s feature representation ability. Too short may lead to insufficient feature information, while too long may introduce irrelevant contextual information, thus affecting model performance.

The final experimental results validated the effectiveness of these hyperparameter configurations, laying the foundation for future model applications.

To comprehensively assess the predictive performance of BioSemAF-BiLSTM, we conducted experiments from two perspectives: data resampling strategies and classifier performance evaluation. To mitigate the issue of data imbalance, multiple up-sampling and down-sampling techniques were explored to enhance the model’s classification performance. Results indicate that SVMSMOTE significantly improves the model’s generalization ability while preserving the original feature distribution, making it the preferred resampling method for this study. Table 3 summarizes the experimental outcomes of various sampling techniques.

It can be seen that there exists a strong gap between the results of the up-sampling and down-sampling methods. To further interpret the significant performance differences observed between them, we conducted a supplementary experiment using a compression-based approximation of Kolmogorov complexity to quantify the information content of the training data under different sampling strategies. Specifically, the original, upsampled, and downsampled training sets were each saved in CSV format and compressed using the gzip algorithm. For practical estimation, we used gzip as a standard lossless compressor to approximate Kolmogorov complexity. Alternative strategies, such as reduced-alphabet encodings or specialized compressors designed for biological sequences (Gusev et al., 2000), could yield slightly different absolute compression sizes, but the relative differences between sampling strategies remain robust. Since Kolmogorov complexity is incomputable in general, the size of the compressed file provides a practical estimate of the dataset’s algorithmic information content. Figure 3 shows the results.

It can be observed that the up-sampling methods produce gzip-compressed files with sizes similar to that of the original dataset. This is because the new samples generated during up-sampling are typically derived from the existing data distribution through interpolation or duplication. Although the number of samples increases, the redundancy in the data also rises, which limits any substantial increase in information complexity. In contrast, down-sampling directly reduces the number of samples, potentially discarding minority class or boundary instances, which results in a significant loss of information. Consequently, the compressed file sizes for down-sampling methods are noticeably smaller than those of the up-sampling methods, indicating a reduction in information complexity.

This finding provides an information-theoretic justification for the performance differences reported in the main experiments: down-sampling reduces the dataset’s descriptive richness and diversity, thereby impairing the model’s generalization capability, whereas up-sampling maintains the complexity while improving class balance (Gao et al., 2024; Thabtah et al., 2020).

In the classifier performance comparison experiment, we selected traditional machine learning algorithms (such as KNN, SVM) and other deep learning models (such as DNN, 2DCNN, BiLSTM) for comparative testing. To ensure fairness, the hyperparameters of each model were adjusted to optimal values. Specifically, the k value of KNN was set to 10, and the penalty coefficient C for SVM was set to 2, consistent with Do’s study (Do et al., 2020). DNN used a four-layer fully connected network, with dimensions of 256, 128, 64, and 32 for each layer. BiLSTM consisted of two layers of bidirectional long short-term memory units, followed by three fully connected layers. 2DCNN included two convolutional layers, two average pooling layers, and one fully connected output layer, with Sigmoid used as the activation function. In these settings, the models were compared according to the auROC and auPRC metrics, with the results shown in Figure 4. The experimental results demonstrate that BioSemAF-BiLSTM performs excellently on both of these key metrics, with an auROC of 0.96, higher than DNN (0.90), BiLSTM (0.93), and 2DCNN (0.87). Its auPRC also reached 0.89, further validating the model’s ability to accurately predict approximately 89% sulfenylation sites.

Comprehensive analysis indicates that the outstanding performance of BioSemAF-BiLSTM can be primarily ascribed to two key factors. Firstly, SVMSMOTE effectively balances the ratio of positive to negative samples, thereby ensuring more robust input data. Secondly, integrating bidirectional long short-term memory units with the adaptive feature fusion module enables the model to comprehensively capture contextual information within protein sequences, leading to enhanced feature learning and classification performance in comparison with conventional machine learning algorithms and other deep learning architectures. These findings underscore the high reliability and generalization capacity of BioSemAF-BiLSTM in sulfenylation site prediction, offering valuable insights for future research.

We employed a combination of word embedding techniques and sequence evolutionary conservation information to generate features, aiming to enhance the model’s predictive performance. Regarding word embedding techniques, fastText was compared with other feature representation methods, such as one-hot encoding and Global Vectors (GloVe) (Pennington et al., 2014). When all methods were evaluated under their respective optimal parameter settings, fastText demonstrated a clear advantage. This result suggests that fastText’s ability to utilize subwords plays a crucial role in improving the performance of neural networks, which is consistent with the findings reported by Do et al. (2020).

In the generation of features based on sequence evolutionary conservation information, this study compared PSSM with BLOSUM62 and amino acid composition (AAC). For BLOSUM62, the substitution scores corresponding to each residue were extracted from the matrix and arranged into per-residue 20-dimensional vectors, resulting in a feature representation of the same shape as the PSSM, which was then fed into the BiLSTM in the same way as the PSSM features. Similarly, under the condition that all parameters were set to optimal values, the PSSM features exhibited superior prediction performance.This finding further suggests a significant correlation between sequence evolutionary conservation information and sulfenylation sites.

The results are illustrated in Table 4. In conclusion, fastText’s subword technique and PSSM features have demonstrated their superiority in semantic feature extraction and biological information representation, respectively, providing a solid foundation for the model’s efficiency and accuracy. This not only validates the importance of multimodal feature generation but also offers strong reference for future research.

To further explore the potential complementarity between semantic and biological features, several fusion strategies were designed and evaluated. Among them, the combination of fastText and PSSM achieved promising performance, indicating that integrating diverse modalities enhances the model’s representation capability.

Notably, further improvements were observed when additional features such as GloVe and BLOSUM62 were included. For instance, the combination of fastText, PSSM, and GloVe achieved the Sn of 92.10, and integrating all six features yielded the best overall performance (ACC: 90.65, Sn: 92.34). These results demonstrate that while deeper feature fusion can be beneficial, the marginal gains may decrease due to redundancy. In addition, it is worth noting that feature fusion increases not only the input dimensionality but also the time cost of feature extraction and model training. In practice, the combination of fastText and PSSM offers a favorable trade-off between performance and computational efficiency.

In order to evaluate the impact of different feature information on model performance, we designed ablation experiments to observe the role of individual features in the prediction task. Table 5 presents the performance metrics when using only PSSM features and only word embedding features. The experimental results show that when using only PSSM features or only word vector features as inputs, the model’s prediction performance decreases. However, when both PSSM and word vector features are used together, the model’s metrics reach their optimal values.

This conclusion indicates that information from a single aspect (such as sequence evolutionary conservation or semantic features) cannot fully represent the global statistical patterns of the data. Relying solely on one feature may overlook some key information. PSSM features capture the biological evolutionary conservation of protein sequences, revealing their potential functional importance, while word vector features, through subword modeling, effectively represent local semantic relationships in sequences. Therefore, combining bioinformatics features (PSSM) with semantic features (word vectors) allows for a more comprehensive integration of the biological significance and statistical properties of the sequence, significantly improving the model’s prediction ability. For the optimal combination of fastText and PSSM features, the estimated information loss during feature extraction is approximately 11%, which is considered acceptable for maintaining sequence-level biochemical relevance (Gusev et al., 2000). This supports the reliability of the extracted features for downstream sulfenylation site prediction. Through the ablation experiment, we further validated the necessity and effectiveness of multidimensional feature fusion, providing theoretical support for building more robust prediction models.

We compared the performance of the BioSemAF-BiLSTM model with existing prediction tools. To ensure the accuracy and fairness of the results, three state-of-the-art prediction tools were selected: iSulf-Cys, Sulf_FSVM, and fastSulf-DNN. These tools were trained and tested using the same dataset as in this study. Additionally, to maintain consistency with previous studies, we followed the standard evaluation protocol of the iSulf-Cys benchmark (Xu et al., 2016), reporting results from both repeated 10-fold cross-validation on the training set and evaluation on the independent gold-standard test set provided with iSulf-Cys. Within each training set, 10% of the data was further set aside as a validation subset for hyperparameter tuning and early stopping, ensuring that model optimization was not influenced by the final test fold. To reduce randomness due to single partitioning, the 10-fold cross-validation process was repeated 10 times with different random seeds.

These three tools use different feature generation methods and model architectures.

• fastSulf-DNN: Uses only biological subwords as features and adopts a relatively simple deep learning network architecture.

Table 6 reports the mean and standard deviation (±) of evaluation metrics in 10-fold CV and illustrates the results in independent test. The experimental results show that BioSemAF-BiLSTM outperforms all of these tools in all evaluation metrics. Specifically, on the MCC metric, it exceeds the next highest by about 17% on the independent test; on the ACC metric, it outperforms the next highest method by 12%; on the Sn metric, it surpasses the next highest method by 5%; and on the Sp metric, it exceeds the next highest method by 18%. Furthermore, despite integrating multidimensional features, the network architecture of BioSemAF-BiLSTM is not overly complex, indicating that the model design achieves improved performance while maintaining a certain level of simplicity.