SAGEPhos: Sage Bio-Coupled and Augmented Fusion for Phosphorylation Site Detection

Dear reviewer 1xkZ,

Thank you for raising these valuable points. Our responses are as follows:

1.Do you consider directly comparing SAGEPhos with a model that simply combines sequence and structure data? It seems the improvements in SAGEPhos may result more from the utilization of structural information itself rather than from the hierarchical fusion strategy. A direct fusion comparison would clarify the benefits of the selective approach.

We have indeed conducted systematic comparisons between SAGEPhos and baseline models with different feature combination approaches. As presented in our ablation study (Section 4.4, Table 3), we evaluated three distinct fusion strategies:

(1) Direct concatenation of substrate sequence, substrate structure and kinase sequence features (w/o fusion);

(2) Direct concatenation of substrate sequence and structure, then use inter-modal fusion between substrate feature and kinase sequence (w/- inter fusion);

(3) Our proposed hierarchical fusion approach in SAGEPhos (w/- fusion & w/- empha).

For your convenience, we have extracted the relevant experimental results for these three strategies:

The experimental results reveal that the simple concatenation baseline (w/o fusion) achieves notably lower performance across all metrics compared to SAGEPhos (e.g., Accuracy: 65.8% vs 80.6%). Even using inter-modal fusion strategy only (w/- inter fusion) shows limited improvement. This performance gap suggests that our Bio-Augmented fusion strategy (fusion strategy between substrate sequence and substrate structure) not only effectively leverages structural information but also successfully filters noise while preserving essential structural features. These results demonstrate that both the incorporation of structural information and our specialized fusion mechanism contribute synergistically to the model's optimal performance.

2.Do you consider other structure encoders besides R-GCN, there are multiple encoding models available for protein structure.

We selected R-GCN as our graph encoder based on its unique ability to model diverse residue-residue relationships in phosphorylation site prediction. Through its relation-specific weight matrices, R-GCN can capture and learn the distinct importance of different types of residue relationships (sequential, spatial, and k-nearest neighbor connections).

To validate this choice, we conducted comprehensive comparisons with representative graph encoders:

R-GCN consistently outperforms alternatives across all metrics, with notable improvements in accuracy (+2% compared to GIN), AUC-ROC (+2.2% compared to GCN), and a significant reduction in false positive rate (21.7% vs ≥23.7%). While GCN, GIN, and GAT are effective in many graph learning tasks, their uniform treatment of residue relationships limits their ability to capture the nuanced structural information critical for phosphorylation site prediction.

3.How to explain the performance differences in FPR across models in Table 2? The substantial FPR variance raises questions about the stability and consistency of SAGEPhos compared to other baselines. A more detailed discussion of the factors behind these differences would be helpful.

SAGEPhos, Phosformer, and Phosformer-ST are designed to predict phosphorylation by considering both kinase and substrate information, while other baselines only utilize substrate information for prediction.

The FPR variations can be primarily attributed to the incorporation of kinase information in SAGEPhos and Phosformer-ST. Compared to substrate sequences, kinase information is relatively sparse, leading to performance fluctuations, including FPR, particularly in kinase cold-start scenarios. This effect is evident even in Phosformer-ST despite its high FPR values. While Phosformer also considers kinase information, its strong bias toward negative predictions results in extremely low FPR with minimal variation, compromising its predictive value.

As for the substrate-only baselines (MusiteDeep, DeepPhos, PhosIDN, PhosIDNSeq), only MusiteDeep and PhosIDN series show relative stability since they only consider substrate information, which is more abundant and thus less affected by cold-start partitioning. However, their higher FPR values, especially PhosIDN, indicate inferior performance. Based on this analysis, the observed FPR variations are natural fluctuations caused by kinase absence, and SAGEPhos demonstrates reliable stability and consistency in its overall performance.

Dear reviewer dmnR,

Thank you for raising these valuable points. Our responses are as follows:

1.What steps did you take to ensure the model’s generalizability, given that it was only evaluated on a limited set of datasets?

As reported in our manuscript Appendix C, Table A1, we validated our model's generalizability on CDK17, a kinase absent from our training set, using substrate data from an independent dataset. We have now extended this evaluation to include four additional novel kinases (MYLK4, CDKL1, PHKG2, SRPK3) [1]. Comparing with MusiteDeep, the best-performing baseline, our model demonstrates strong generalization capability across different kinases:

[1] Johnson et al., 2023, An atlas of substrate specificities for the human serine/threonine kinome.

2.Could you clarify why you chose complex models like R-GCN for structural representation over simpler alternatives, and what benefits they provide?

Our choice of R-GCN is motivated by the complex nature of phosphorylation site prediction. Unlike simpler graph models, R-GCN uses relation-specific transformation matrices to effectively capture and learn the distinct importance of different types of residue relationships (sequential, spatial, and k-nearest neighbor connections) in phosphorylation site prediction. Comprehensive comparisons with representative graph models (GCN, GIN, and GAT) demonstrate R-GCN's superior performance across all evaluation metrics:

3.How does SAGEPhos handle incomplete or low-confidence structural data, given its reliance on AlphaFold2 predictions?

To handle varying confidence levels in AlphaFold2's structural predictions, we propose a Dynamic fusiOn meChanism (DOC) that employs a gated module to selectively filter structural information based on relevance. Instead of using hard thresholds to discard low-confidence structures, DOC preserves valuable features while reducing the impact of noisy predictions. Comparative experiments with commonly used pLDDT thresholds (0.5 and 0.7) demonstrate SAGEPhos's superior performance across all metrics:

4.Why does the model only use both sequence and structure data for substrates, while only sequence data is used for kinases?

We employ different input features for kinases and substrates based on their distinct biological properties. As shown in previous studies [1], kinases are highly conserved enzymes with similar core structures and functional sites, making their sequence information sufficient for prediction purposes. Furthermore, not all kinase sequences have corresponding structural information available - incorporating both modalities would significantly reduce our training dataset size (from 29,376 to 24,893 samples), which consequently leads to decreased model performance, as shown below:

In contrast, the high sequence and structural diversity of substrates, which critically influences phosphorylation events, necessitates the integration of both features for accurate prediction.

[1] Taylor et al., 2011, Protein kinases: evolution of dynamic regulatory proteins.

5.Why was a GCN applied to the AlphaFold structures instead of directly using AlphaFold’s embeddings?

Given the resource requirements of AlphaFold, we took a more lightweight approach by utilizing pre-computed structures from the AlphaFold Database and processing them with R-GCN. The R-GCN's relation-specific architecture enables explicit modeling of various molecular interactions, which is essential for capturing diverse structural patterns in kinase-substrate interactions.

Specifically, due to resource constraints, AlphaFold2's intensive computational requirements make it impractical to directly generate embeddings for our dataset of 30,000+ sequences, especially within the limited rebuttal period. Our current approach is significantly more efficient while still effectively capturing structural information. Thank you for suggesting this potential direction, we plan to explore it in future work when sufficient computational resources become available.

Dear reviewer dmnR,

We sincerely appreciate your positive feedback and are delighted that our revisions have addressed your concerns. We believe our work brings novel perspectives to post-translational modification analysis, particularly in phosphorylation studies, and we are committed to continuing our research endeavors in this important field. If possible, we would be grateful if you could consider increasing the score.

Thank you again for your valuable guidance throughout this review process.

Dear reviewer aM45,

Thank you for raising these valuable points. Our responses are as follows:

1.It is unclear how the cold-start and warm-start experiments are designed. The cited reference does not include this information either.

For warm-start setting, the dataset was randomly split into training, validation, and testing sets with a ratio of 8:1:1;

For kinase cold-start, data points with the same kinase were grouped into the same set to prevent kinase overlap across train/valid/test sets;

For substrate cold-start, data points with identical substrate sequences were grouped together to ensure no sequence overlap between sets.

These experimental designs evaluate model generalization across random, unseen kinase, and unseen substrate scenarios.

2.Are negative sequences chosen from other positions within the same substrate of the kinase, or the entire substrate set across all kinases? Are these positions known phosphosite locations on other substrates, or are they potential sites that could accept a phosphate group without a reported kinase association?

Our negative sequences are drawn from the entire substrate set across all kinases, comprising two types:

(1) Known phosphosites that are validated substrates of other kinases;

(2) Experimentally verified phosphorylatable sites without reported kinase associations (majority)

For each kinase, we maintained a 1:1 ratio between positive examples (confirmed phosphorylation sequences) and randomly selected negative examples from these two sources. We have added detailed description of this negative example selection procedure to Appendix A "DATASETS AND PARTITION METHOD".

3.When splitting the train/test folds, is substrate similarity or kinase similarity taken into account?

Rather than using sequence similarity thresholds, we enforce strict data separation by grouping identical kinases or substrates into the same split named cold-start. To test model generalization, we conducted additional experiments using sequence similarity-based splits with MMseqs2 (50% similarity threshold, 8:1:1 ratio) for both kinase sequences and substrate sequences respectively.

Our results demonstrate that sequence similarity-based clustering achieves relatively better performance (kinase clustering: Acc 78.8%, AUC-ROC 86.4%; substrate clustering: Acc 79.2%, AUC-ROC 87.3%) compared to the more challenging cold-start settings (kinase cold-start: Acc 68.6%, AUC-ROC 76.9%; substrate cold-start: Acc 79.1%, AUC-ROC 87.2%). This performance gap validates that our cold-start splitting strategy poses a more rigorous evaluation scenario than sequence similarity-based splits. As expected, these results are lower than those obtained under random split settings (Acc: 80.6%, AUC-ROC: 88.3%), which aligns with the intuition that performance naturally decreases as the splitting criteria become more stringent.

4.ESM2 is used for feature extraction, but the specific ESM2 model size is not specified.

We used ESM-2-650M for feature extraction in our experiments. This has been added to Appendix B “IMPLEMENTATION DETAILS” in revised manuscript.

5.How are the conservation scores computed?

The conservation scores in our method are implemented as learnable embeddings specifically for phosphorylation sites. This design is motivated by the biological observation that different positions in protein sequences exhibit varying degrees of conservation, with functionally critical sites typically being more conserved. In particular, the central residue (S/Y/T) within the 11-mer peptide sequences represents a highly conserved motif for phosphorylation events. These scores are dynamically optimized during model training to capture the relative importance of these conserved sites.

The effectiveness of this approach is demonstrated in Section 4.4 "ABLATION STUDY", where we compare model performance with and without the conservation score emphasis ("w/- fusion & w/- empha" vs. "w/- fusion & w/o empha").

6.Which physicochemical properties are used and with what kind of representation? Continuous scale, binarized categories?

We employ four fundamental physicochemical properties commonly used in protein structure-function analysis: (1) Aliphatic, (2) Aromatic, (3) Acidic charged, and (4) Basic charged. Each property is represented as a binary category. These basic properties capture essential amino acid characteristics that influence protein structure and interactions. We have added this description to Appendix B "IMPLEMENTATION DETAILS AND HYPERPARAMETERS" in revised manuscript.

7.Figure A1 needs a legend for the color scale.

We have added a legend of the color scale in Figure A1.

8.In CDK17, how many substrates are taken from the bottom and the top?

After structural matching and balancing the dataset to maintain a 1:1 ratio between positive and negative samples, our final test set consisted of 46 positive samples (ranked first) and 46 negative samples (ranked at the bottom).