Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

Genomics Proteomics Bioinformatics

[Link]/locate/gpb
[Link]

REVIEW

A Review of Machine Learning and Algorithmic

Methods for Protein Phosphorylation Site
Prediction
Farzaneh Esmaili 1,#, Mahdi Pourmirzaei 1,#, Shahin Ramazi 2,*,
Seyedehsamaneh Shojaeilangari 3, Elham Yavari 1

1
Department of Information Technology, Tarbiat Modares University, Tehran 14115-111, Iran
2
Department of Biophysics, Faculty of Biological Sciences, Tarbiat Modares University, Tehran 14115-111, Iran
3
Biomedical Engineering Group, Department of Electrical Engineering and Information Technology, Iranian Research
Organization for Science and Technology (IROST), Tehran 33535-111, Iran

Received 16 March 2022; revised 16 January 2023; accepted 23 March 2023

Available online 19 October 2023

Handled editor Yu Xue

KEYWORDS Abstract Post-translational modifications (PTMs) have key roles in extending the functional diver-
Phosphorylation; sity of proteins and, as a result, regulating diverse cellular processes in prokaryotic and eukaryotic
Machine learning; organisms. Phosphorylation modification is a vital PTM that occurs in most proteins and plays a
Deep learning; significant role in many biological processes. Disorders in the phosphorylation process lead to mul-
Post-translational tiple diseases, including neurological disorders and cancers. The purpose of this review is to orga-
modification; nize this body of knowledge associated with phosphorylation site (p-site) prediction to facilitate
Database future research in this field. At first, we comprehensively review all related databases and introduce
all steps regarding dataset creation, data preprocessing, and method evaluation in p-site prediction.
Next, we investigate p-site prediction methods, which are divided into two computational groups:
algorithmic and machine learning (ML). Additionally, it is shown that there are basically two main
approaches for p-site prediction by ML: conventional and end-to-end deep learning methods, both
of which are given an overview. Moreover, this review introduces the most important feature
extraction techniques, which have mostly been used in p-site prediction. Finally, we create three test
sets from new proteins related to the released version of the database of protein post-translational
modifications (dbPTM) in 2022 based on general and human species. Evaluating online p-site pre-
diction tools on newly added proteins introduced in the dbPTM 2022 release, distinct from those in
the dbPTM 2019 release, reveals their limitations. In other words, the actual performance of these

* Corresponding author.
E-mail: [Link]@[Link] (Ramazi S).
#
Equal contribution.
Peer review under responsibility of Beijing Institute of Genomics, Chinese Academy of Sciences / China National Center for Bioinformation and
Genetics Society of China.
[Link]
1672-0229 Ó 2023 The Authors. Published by Elsevier B.V. and Science Press on behalf of Beijing Institute of Genomics, Chinese Academy of Sciences /
China National Center for Bioinformation and Genetics Society of China.
This is an open access article under the CC BY license ([Link]
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1267

online p-site prediction tools on unseen proteins is notably lower than the results reported in their
respective research papers.

Introduction Phosphorylation is present in more than one-third of

human proteins, and this modification is regulated by approx-
imately 568 human PKs and 156 protein phosphatases [29]. In
Post-translational modifications (PTMs) are biochemical reac-
this sense, phosphorylation is one of the widest spread and
tions occurring on a protein after its translation [1,2], which
most extensively studied protein PTMs, which has a significant
change the regulated physicochemical properties, maturity,
role in the control of biological processes like proliferation,
and activity of most proteins [3,4]. It includes cutting, folding,
differentiation, and apoptosis [29,34]. Site mutations or dys-
ligand-binding, adding a modifying group to one or more
regulation of kinase activity, their hyperactivity, malfunction,
amino acids, and finally changing the chemical nature of
or overexpression, and also the hyperphosphorylation of
amino acids [5,6]. In recent years, an increasing volume of
human proteins are associated with certain disease states such
PTM data has become available because of improvements in
as cancer, Alzheimer’s disease (AD), Parkinson’s disease (PD),
mass spectrometry (MS) based on high-throughput proteomics
frontotemporal dementia (FTD), and various pathways
[7]. There are more than 600 types of PTMs [8] that affect
involving the immune system [27–29,35]. Therefore, identifying
many aspects of cellular functionalities, such as metabolism,
kinase-specific p-sites is essential for understanding the regula-
signal transduction, activity, stability, and localization of var-
tory mechanisms of phosphorylation.
ious proteins [9,10]. Recent studies have shown that each mod-
Multiple experimental methods are used for the detection
ification leads to a multitude of effects on the structure and,
assays of protein phosphorylation, like liquid
therefore, the function of the proteins [11]. PTMs include
chromatography-tandem mass spectrometry (LC-MS/MS),
phosphorylation, glycosylation, ubiquitination, sumoylation,
radioactive chemical labeling, and immunological detection
acetylation, succinylation, and nitrosylation, as well as numer-
[such as proximity ligation assay (PLA), chromatin immuno-
ous others involved in most cellular activities [9,12–16]. More-
precipitation, and Western blotting], although the combina-
over, PTMs play key roles in a variety of biological regulatory
tion of the LC-MS/MS method with the
pathways, including metabolic pathways, DNA damage
immunoprecipitation (IP) strategy is suitable for the detection
response, transcriptional regulation, signaling pathways, pro-
of p-sites in proteins [36,37]. Cellular processes are regulated
tein–protein interactions (PPIs), apoptosis, cell death, insulin
by phosphorylation, which is highly conserved and signifi-
signaling, immune response, and aging [17,18]. The dysregula-
cantly affects protein stability. Figure 1 shows the schematic
tion of PTMs is linked to diseases such as cancer, diabetes, car-
of protein phosphorylation process.
diovascular disease, and neurological disorders [19–24].
The number of known p-sites has grown since 2003, and it
Phosphorylation is one of the most important reversible
rose from 2000 to more than 500,000 in experimental data-
PTMs. Phoebus Levene discovered phosphorylation in 1906
bases. The p-sites are involved in the regulation of at least
in the protein vitellin (phosvitin) [25]. In phosphorylation, a
30% of human proteomes [38,39].

2 phosphate group is covalently added to serine (S), thre-
Indeed, experimental approaches are generally difficult,
onine (T), tyrosine (Y), and histidine (H) residues and then
slow, and costly, and need specialized equipment and knowl-
removed by protein phosphatases. It is known that protein
edge. Over the last two decades, PTM research has made
phosphorylation regulates the activity of various enzymes
remarkable progress due to technological advancements and
and receptors, including signal pathways [26], and can greatly
the emergence of new computational methods.
impact the folding, function, stability, and subcellular localiza-
A previous study [36] has reviewed PTM tools, resources,
tion of the protein [25,27,28]. In eukaryotes, this modification
and related databases and also investigated the challenges of
plays a vital role in signal transduction and other biological
algorithmic methods. The authors divided 10 types of PTMs
functions, including protein synthesis, cell division, signal
into small chemical groups, lipids, and small proteins, and they
transduction, DNA repair, environmental stress response,
investigated databases and algorithmic approaches for differ-
transcriptional regulation, apoptosis, cellular motility, immune
ent PTM sites [36]. Shi et al. [40] reviewed 19 available tools
response, metabolism, cell growth, development, cellular dif-
for phosphorylation networks. They reported different analy-
ferentiation, and aging [29,30].
ses for their functionality, data sources, performance, network
In eukaryotes, the phosphorylation process is catalyzed via
visualization, and implementation. Rashid et al. [41] reviewed
protein kinases (PKs) differentially and specifically in which
specified machine learning (ML) methods, main feature selec-
each PK only modifies a subset of substrates to ensure signal-
tion methods, databases, and current online tools for microbial
ing fidelity [31]. The phosphorylation process typically occurs
p-sites. They only investigated microbial p-sites and did not
on T, S, and Y residues. However, in prokaryotes, plants,
mention other p-sites in organisms nonetheless. Also, their
and fungi, this modification also occurs on H and aspartate
work was limited to classical ML methods.
(D) residues in protein sequences, which play important roles
In this study, unlike other previous studies, we investigated
in two-component signaling systems [32,33]. Indeed, phospho-
all features, databases, and methods concerning p-site predic-
histidine (PhosH) is involved in many biological processes,
tion. First, valid PTM databases that contain phosphorylation
such as central metabolism and signaling in eukaryotes and
experimental data were introduced. Then, the two most impor-
bacteria, making it important to develop and potentially
tant phosphorylation databases were reviewed, in which the
improve computational tools for accurately predicting his-
numbers of organisms and p-sites were covered in detail.
tidine phosphorylation sites (p-sites) [32,33].
1268 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

Figure 1 Schematic of protein phosphorylation process

The red balls represent amino acids. ATP, adenosine triphosphate; ADP, adenosine diphosphate.

Furthermore, two main data preparation steps for p-site data- Considering different types of PTMs, databases are
sets, including data collection and data preprocessing, were arranged into specific and general terms, in which general
reviewed. In other words, this study investigated methods for PTM databases cover a wide domain of data for different types
data collection and also introduced the most important and of PTM, but specific databases are constructed based on spe-
functional approaches for data preprocessing. Additionally, cial types of PTMs like phosphorylation.
all evaluation metrics that have been used for p-site prediction Databases such as dbPTM [43], SysPTM [3], Swiss-Prot
were introduced. Then, the most common and important fea- [44], and HPRD [42] are general databases that cover different
ture extraction methods were described based on the physico- types of PTMs, and phosphorylation is one of them. On the
chemical, sequence, evolutionary, and structural properties of other hand, the eukaryotic phosphorylation site database
amino acids. It was found that there are generally two ML- (EPSD) [45], LymPHOS2 [46], Phospho3D [47], Phospho.
based approaches for p-site prediction, which are divided into ELM [48], and Regulatory Network in Protein Phosphoryla-
conventional ML and end-to-end deep learning (DL) methods. tion (RegPhos) [49] are specifically gathered for p-sites.
In the present study, the methods of both approaches were In the following, two important databases, dbPTM and
reviewed, and the available online tools for p-site prediction EPSD, which are known as general and specific databases
were briefly introduced. for p-site data, are going to be introduced. Furthermore,
Finally, we created three test sets from new proteins related Table 1 summarizes both general and specific databases
to the released version of the database of protein post- according to their statistical information for p-sites.
translational modifications (dbPTM) in 2022, and then evalu-
ated and compared the available online tools together in differ-
EPSD
ent metrics on the three specific test sets.

EPSD is one of the most specific and comprehensive databases

Databases
for p-sites which has been updated in 2020. EPSD has been
updated from two databases, the database of p-sites in plants
Developing a prediction model requires a dataset of (dbPPT) [58] and the database of p-sites in animals and fungi
experimentally-validated p-sites. Therefore, the availability of (dbPAF) [59], which includes roughly  82,000 p-sites from
general and specific databases for p-sites is the first step toward 20 plants and more than 483,000 p-sites from 7 different types
this end [36]. Databases are constantly evolving due to the of animals and fungi. Moreover, EPSD collects p-sites from 13
advent of technology and may be updated by providing accu- additional databases, including PhosphoSitePlus [50], Phos-
rate details. These databases contain information about differ- [Link] [48], Universal Protein Resource (UniProt) [60],
ent organisms, such as viruses, animals, and plants, that has PhosphoPep [61], Biological General Repository for Interac-
been collected manually and experimentally. For instance, all tion Datasets (BioGRID) [62], dbPTM, Fungi Phosphorylation
the information in the Human Protein Reference Database Database (FPD) [63], HPRD, Medicago PhosphoProtein Data-
(HPRD) has been collected manually, and it contains more base (MPPD) [64], Plant Protein Phosphorylation Database
than 95,000 p-sites extracted from  13,000 proteins [42]. (P3DB) [57], phosphorylation site database (PHOSIDA) [53],
Table 1 Summary of phosphorylation databases
Type Acronym General statistic Type of data and database URL
No. of covered organisms No. of p-sites and phosphoproteins
General dbPTM [43] More than 1000 organisms p-site:  1,770,000 Experimentally-validated and predicted; [Link]
database P:  557,700 secondary
PhosphoSitePlus [50] 26 organisms p-site:  240,000 Experimentally-validated; [Link]
P:  20,200 primary
p-site:  316,500

Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction

PTMcode v2 [51] 19 organisms Experimentally-validated; [Link]
P:  45,300 secondary
qPTM [52] Human p-site:  199,000 Experimentally-validated; [Link]
P:  18,402 secondary
HPRD [42] Human p-site:  1100 Experimentally-validated; [Link]
P:  30,000 primary
PHOSIDA [53] 9 organisms p-site:  70,000 Experimentally-validated; –
P:  28,700 secondary
PTM-SD [1] 7 model organisms p-site:  1600 Experimentally-validated; [Link]
P:  842 secondary dsimb_tools/PTM-SD
SysPTM 2.0 [3] 6 organisms p-site:  353,000 Experimentally-validated; –
P:  53,200 secondary
EPSD [45] 68 organisms p-site:  1,616,800 Experimentally-validated; [Link]
P:  209,300 secondary
Specific PhosphoNET [54] Human p-site:  966,000 Experimentally-validated and predicted; [Link]
database for P:  20,000 secondary
phosphorylation RegPhos 2.0 [49] Human, mouse, and rat p-site:  113,000 Experimentally-validated and predicted; [Link]
P:  18,700 secondary
dbPSP [55] 200 prokaryotic organisms p-site:  19,300 Experimentally-validated; [Link]
P:  8600 secondary
pTestis [56] Mouse p-site:  17,800 Experimentally-validated and predicted; –
P:  3900 secondary
LymPHOS [46] Human and mouse p-site:  18,300 Experimentally-validated and predicted; [Link]
P:  4900 primary
P3DB [57] 45 plant organisms p-site:  220,000 Experimentally-validated and predicted; [Link]
P:  57,000 secondary
Note: This table was adapted from Table 1 [36] with permission. The type of database can be secondary or primary: secondary databases are the integration of other databases, and primary databases
are independent. p-site, phosphorylation site; P, phosphoprotein; URL, uniform resource locator; PTM, post-translational modification; dbPTM, database of protein post-translational modifications;
HPRD, Human Protein Reference Database; PHOSIDA, phosphorylation site database; PTM-SD, post-translational modification structural database; EPSD, eukaryotic phosphorylation site
database; RegPhos, Regulatory Network in Protein Phosphorylation; P3DB, Plant Protein Phosphorylation Database.

1269
1270 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

A No. of p-sites in animals

B No. of p-sites in animals, plants, and fungi

Animal
Fungus
Plant

Figure 2 EPSD database analysis

A. Number of p-sites in the animal proteins distributed by different animal species. B. Number of p-sites in proteins from animals, fungi,
and plants. All figures are based on EPSD. p-site, phosphorylation site; EPSD, eukaryotic phosphorylation site database.

Arabidopsis Protein Phosphorylation Site Database (Phos- 130 types of PTMs in more than 1000 organisms [36]. The
PhAt) [65], and SysPTM [38]. Totally, this database contains new version of dbPTM [66] in 2022 has curated more than
 1,616,800 experimentally known p-sites in  209,300 phos- 2,777,000 PTM sites from 41 published databases and
phoproteins of 68 eukaryotes (18 animals, 24 plants, 19 fungi,  82,000 research articles.
and 7 protists) [45]. Figure 2 shows the number of p-sites in dif-
ferent animal species (Figure 2A), and also depicts the distribu- Identifying driver mutations and their effects on p-sites
tion of p-sites in animals, fungi, and plants (Figure 2B) in the
EPSD database. Figure 3 shows the number and the percentage Phosphorylation is involved in many aspects of cellular organi-
of S, T, and Y p-sites in the EPSD database for animal, fungus, zation and signaling pathways associated with diseases. Vari-
and plant species. Figure 4 shows the number of p-sites in dif- ous studies have demonstrated that p-sites are evolutionarily
ferent plant and fungus species in the EPSD database. constrained in human genomes, as well as prevalent in cancer
driver mutations and causal variants of inherited diseases.
dbPTM Therefore, phosphorylation information and knowledge of
its function are useful for interpreting genetic variations, geno-
A general database called dbPTM integrates PTM’s data from type–phenotype associations, and molecular diseases and their
30 databases and  92,600 research articles. dbPTM covers treatment [67].
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1271

Figure 3 Distribution of phosphorylated S, T, and Y in EPSD

Number (A) and percentage (B) of phosphorylated S, T, and Y in proteins from animals, plants, and fungi in EPSD. S, serine;
T, threonine; Y, tyrosine.

A p-sites in fungi B p-sites in plants

Log10 No. of p-sites

Log10 No. of p-sites

subsp.

subsp.
subsp.

subsp.

Figure 4 p-site distribution in fungi and plants in EPSD

Number of p-sites in figui (A) and plants (B) in EPSD.

The most common type of sequence change, DNA single and updated tools are required to interpret rapidly increasing
nucleotide variants (SNVs), is caused by a single nucleotide genomic and phosphoproteomic data to explain the signaling
change. Genetic variation of p-sites via SNVs can have an networks. We are briefly going to describe the ActiveDriverDB
effect directly by modifying target residues or indirectly by database as well as mutation impact on phosphorylation
modifying the consensus binding sequences (i.e., short linear (MIMP) and PTMsnp tools in this field [67–69].
motifs) located in the flanking sequences of phosphorylated ActiveDriverDB is a web database that was designed to
residues. As a result, this can change signaling networks by understand how protein coding varies in the human genomes.
making, changing, and disrupting the p-sites [68]. There have The ActiveDriverDB database contains more than 260,000
been reports of phosphorylation-related SNVs that disrupt experimentally identified PTM sites in human proteome using
existing sites, create new sites, disturb kinase–substrate interac- public databases like PhosphoSitePlus, UniProt, Phospho.
tions, and cause disease phenotypes. A major challenge faced ELM, and HPRD, which contains  149,300 p-sites
by biomedical research is the identification of genotype–phe- [42,48,50,60]. As evidenced in the ActiveDriverDB database,
notype associations, molecular mechanisms, and cancer driver changes in target amino acid substitutions in p-sites influence
mutations [67]. the creation of pathogenic disease mutations, somatic muta-
There are various databases with a useful list of genome tions in cancer genomes, and germline variants in humans.
variants in p-sites and other PTM sites. However, they provide Additionally, the ActiveDriverDB database contains phospho-
no perspective on how mutations on p-sites and other protein proteomics data reflecting the cellular response to severe acute
sites will affect kinase binding [67–69]. Therefore, databases respiratory syndrome coronavirus 2 (SARS-CoV-2) infection,
1272 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

which can be used to predict the impact of human genetic vari- Data preprocessing
ation on SARS-CoV-2 infection and coronavirus disease 2019
(COVID-19) disease course [68]. After constructing the primary positive and negative datasets,
An online tool called MIMP ([Link] one important task is removing inconsistent or redundant sam-
can be used for predicting kinase–substrate interactions based ples to gain a more reliable dataset.
on missense SNVs. MIMP analyzes kinase sequence specifici- Cluster database at high identity with tolerance (cd-hit)
ties and predicts whether SNVs disrupt the existing p-sites or program is a protein clustering program widely used to reduce
create new ones. This helps discover mutations that modify the sequence homology and filter out the similar ones. Accord-
protein function by altering kinase networks and provides ing to different phosphorylation prediction studies [70,72–74],
insights into disease biology and therapy development [67]. a threshold of identity is considered to range from 30% to
PTMsnp is another online tool for identifying driver genetic 60% in many phosphorylation prediction studies [75].
mutations aiming at PTM sites in proteins across different There are three main steps in the literature for removing
cohorts of TCGA by using a Bayesian hierarchical model. inconsistent or redundant proteins [36,76]. First, redundant
There are more than 411,500 modification sites in PTMsnp phosphoproteins should be removed by using the cd-hit pro-
from 33 different types of PTMs and 1,776,800 mutation sites gram. In the second step, identical subsequences are removed
from 33 types of cancer. The web server detects proteins with a within positive and negative sets by selecting the optimal win-
higher frequency of PTM-specific mutations in the motif dow size. Finally, identical subsequences between the positive
region, considered to be the key targets in human disease and negative datasets are removed by choosing the size of
development [69]. the optimal window size.

Data gathering and preprocessing Class-imbalance problem

In this section, we are going to describe steps concerning cre- It is a common problem in ML when there is an imbalance
ating and preprocessing datasets before p-site prediction. In between the distribution ratios of data classes. In other words,
the last decade, due to the importance of phosphorylation in a dataset that has unequal samples in classes is imbalanced.
understanding the biological systems of proteins and in guid- This is not an issue when the difference is not that much. Nev-
ing basic biomedical drug design, research on phosphorylation ertheless, when one or more classes are infrequent, many mod-
has boomed. Several experimental methods are used to identify els do not work well at identifying the minority classes. For
p-sites in a large number of phosphorylation examples with example, in p-site prediction, preprocessed datasets are mostly
high accuracy (ACC), but many of them are labor-intensive imbalanced because the number of negative samples is much
and time-consuming. Therefore, low-cost and fast algorithmic greater than the positive samples. Figure 5 shows the prepro-
and ML methods have become popular to overcome the prob- cessing framework for balancing data.
lems associated with experimental methods [70]. In order to There are three most commonly used approaches to deal
build a dataset for p-site prediction, all verified data from mul- with class-imbalance problems: upsampling, downsampling,
tiple databases are considered. Mainly, there are two main and customized loss function.
steps to prepare a dataset: data collection and data preprocess-
ing [36,70]. Upsampling
It generates additional data for minority classes either by mak-
Data collection ing copies of the minimum class or by creating synthetic data
which can represent samples of minimum classes.
Data collection includes negative and positive data collection
steps. Downsampling
It removes data from the majority class either randomly or
Negative data collection using intelligent approaches of sample selection to handle the
S, T, and Y residues existing in experimentally-validated pep- issue.
tides without any phospho-groups are considered as non-
p-sites or negative samples. There are two major strategies Customized loss function
available to choose the negative samples. Firstly, from phos- This is a technique to deal with imbalance problems in ML
phoproteins, the negative random samples of the target residue that tries to customize the model loss function by assigning lar-
that did not undergo the phosphorylation modifications are ger weights to minority. Customized losses have demonstrated
selected. Secondly, from non-phosphoproteins with none of better performance and attracted more attention than upsam-
their target residues (S, T, and Y) that have undergone specific pling and downsampling approaches [77].
phosphorylation (based on experimental evidence) are selected
as the negative set [36,71].
Evaluation
Positive data collection
The well-known evaluation metrics for protein p-sites are clas-
S, T, and Y residues as p-sites or the positive samples are usu-
sified into five methods: ACC, sensitivity (SN), specificity (SP),
ally compiled from the aforementioned databases (e.g., EPSD
Matthews coefficients of correlation (MCC), and the area
and dbPTM). These samples are usually known from experi-
under the receiver operating characteristic (ROC) curve
ments [36].
(AUROC). These metrics are evaluated with a confusion
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1273

Figure 5 Work flow of data collection and preprocessing

Data preprocessing includes the balancing step. dbPTM, database of protein post-translational modifications.

matrix that summarizes the performance of models; it com- into. That is why the procedure is called k-fold cross validation,
pares the real target values with those predicted by a model. in which specific values for k can be chosen. Considering the
The number of rows and columns in this matrix is based on scenario of 5-fold cross validation (k = 5), a dataset is divided
the number of classes. From the confusion matrix, we will into 5 bunches. Within the first iteration, the primary fold is uti-
end up with four values [36,76]. True positive (TP) represents lized to assess the model, and the rest are utilized to train the
the number of positive samples classified correctly. False pos- model. Within the second iteration, the subsequent fold is uti-
itive (FP) represents the number of negative samples classified lized as the validation set, whereas the rest serve as the training
incorrectly. True negative (TN) represents the number of neg- set. This process is repeated until each fold has been used as the
ative samples classified correctly. False negative (FN) repre- validation set. Each sample is given the opportunity to be uti-
sents the number of positive samples classified incorrectly. lized within the validation set one time as well as utilized to
train the model k
1 times. The k-fold cross validation is usu-
Model evaluation ally used when the amount of train–valid data is limited. On
the contrary, when dealing with huge amounts of data, we do
Basically, there are three methods for model evaluation for not need to have a big valid set. In other words, the proportion
p-site prediction: independent test (train–test), k-fold cross val- of the train–valid split can sometimes go below 1% for the valid
idation, and jackknife cross validation (or leave-one-out cross set. This approach is mostly used when massive amounts of
validation). In the first one, a dataset is split into two sets: a data are accessible. However, in low-data regimes, they usually
train set and a test set. Then, the train set is divided into two split with proportions of 30%–70%.
subsets again: a train subset and a valid subset. The basic pro- Note that there is also another evaluation strategy named
cedure is that the train subset is used to train models, and the the jackknife cross validation test [78], rarely used for p-site
valid subset is used for the evaluation of the trained models. prediction. As the most objective method, the jackknife cross
After selecting the best model with respect to the valid subset validation (or the leave-one-out cross validation) delivers
result, we need to evaluate it on the test set. At the end, we unique results for a dataset in which one sample is selected
should report the test set, and there shouldn’t be much differ- to serve as the test data, whereas the rest are used as the train-
ence between the valid subset and the test set results (Figure 6). ing data. This procedure is repeated N times for a dataset with
On the other hand, there is another assessment strategy uti- N samples, which can be expensive for large datasets [79].
lized to assess ML models on restricted data samples. The In summary, k-fold should be used in low-data regimes, and
method contains a single parameter called k that alludes to an independent method with a small percentage of a test set
the number of bunches that data samples should be divided should be used when we have access to lots of data.
1274 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

Figure 6 Evaluation
The evaluation step can be done by two methods: k-fold cross validation and independent test. Independent test method sometimes is
called ‘‘train–test” or ‘‘train–valid–test” as well. MCC, Matthews coefficients of correlation.

Methods for predicting p-sites some of them, they finally reported a series of patterns as the
output during an iterative cycle.
He et al. [82] showed that the number of patterns to be
In the following sections, we are going to review methods of
examined around each position is growing exponentially based
p-site prediction by dividing them into two main categories:
on the length of the window. They referred to two developed
algorithmic methods and ML. Likewise, ML methods are also algorithms to find phosphorylation patterns, named Motif-X
divided into two approaches: conventional ML methods and
and model-based DL (MoDL) algorithms. They supposed that
end-to-end DL methods.
these algorithms do not detect all patterns, and some patterns
remain hidden from biologists. Therefore, they introduced a
new algorithm called Motif-ALL to discover and report all
Algorithmic methods possible patterns based on previous algorithms.
There has been a family of algorithms called Group-based
Innovative algorithms based on statistical approaches have Prediction System (GPS) for many years as algorithmic meth-
been used in many studies. Here, we need to define algorithmic ods [83–89]. In 2004, an algorithm was developed for group-
methods as computational methods in which there are no based p-site prediction and scoring 1.0, based on the hypothe-
learning algorithms to gain information directly from data. sis that similar short peptides exhibit similar biological func-
Schwartz and Gygi [80] proposed a statistically repetitive tions. Likewise, the algorithm was refined and created an
method, using a set of phosphorylated peptide sequences to online service called GPS 1.1, which could predict p-sites for
extract the patterns and a set of peptide sequences to evaluate 71 PK clusters. Then, GPS 2.0 and 2.1 were presented with
the predictions. They mapped two sets of sequences to the the same scoring strategy using two methods named matrix
position weight matrix so that in the matrices, the number of mutation (MaM) and motif length selection (MLS), which
repetitions of each residue was determined from 6 positions were designed to improve the ACC. Consequently, GPS 2.2,
higher to 6 positions lower than each p-site (which means their 3.0, 4.0, and 5.0 algorithms were developed, which are used
window size for each peptide is 13 amino acids long). Then, for the prediction of other PTM sites rather than p-sites [31].
they formed a binary matrix based on these two matrices. This
final matrix indicates the probability of observing a specific ML methods
residue around a p-site by examining this matrix and compar-
ing it with other p-sites. Most algorithms used for phosphorylation prediction are
Chen et al. [81] presented a new method for predicting based on ML. Moreover, with explosions of the DL method
p-sites by collecting four background datasets, including phos- in the early 2010s, ML has become even more popular than
phorylated and non-phosphorylated sequences. They chose a before. ML is generally the ability of machines to do actions
given length of 13 amino acids for windows around p-sites. Ini- based on prior knowledge and experience [90]. There are more
tially, they formed the position weight matrices and then than 40 different methods for predicting p-sites, and many of
extracted the patterns. By scoring those patterns and deleting them are based on ML techniques, including logistic regression
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1275

(LR), support vector machine (SVM), random forest (RF), Average accumulated hydrophobicity
and k-nearest neighbor (KNN) [70]. Average accumulated hydrophobicity (ACH) quantifies the
In general, there are two main strategies in ML to predict tendency of amino acids surrounding S, T, or Y residues to
phosphorylation: conventional ML methods and end-to-end be exposed to solvent [96]. For different window sizes, ACH
DL methods. The conventional approach stands for using is calculated by averaging the cumulative hydrophobicity
ML algorithms as a part of solving a solution besides other indices around the p-site. Note that every site is located in
steps in a pipeline design such as feature extraction and the center of the sliding windows [97,98].
hand-feature engineering. In other words, usually, in a conven-
tional ML-based system, there are multiple stages of process- Encoding scheme based on attribute grouping
ing which need to be designed individually. However, the The encoding scheme based on attribute grouping (EBAG)
end-to-end DL approaches can replace all those steps with a represents the hydrophobicity attribute of the amino acids
single neural network. This type of learning tries to eliminate and divides the residues into 4 classes based on their physico-
the need for explicit feature engineering steps inside the learn- chemical properties: hydrophobic class c1 = {A, F, G, I, L, M,
ing system by feeding raw data as the input to it. P, V, W}, polar class c2 = {C, N, Q, S, T, Y}, acidic class c3 =
{D, E}, and basic class c4 = {H, K, R} [99,100].
Feature extraction
Overlapping property
Overlapping property (OP) clusters each protein based on its
In protein phosphorylation prediction, various types of con-
chemical attributes. Each amino acid is classified into 10
ventional approaches have been studied. Feature extraction
physicochemical properties: polar, positive, negative, charged,
is an important step in those approaches [91]. In this review,
hydrophobic, aliphatic, aromatic, small, tiny, and proline [98].
we summarized 20 feature extraction techniques suggested
according to the physicochemical, sequence, evolutionary,
Pseudo amino acid composition
and structural properties of amino acids. We have tried to
Pseudo amino acid composition (PseAAC) is first defined by
introduce the most important and practical methods of feature
Chou et al. [101] for coding proteins. They proposed sequence
extraction in the following.
order and physicochemical information in protein sequences.
For more details, refer to [102–105].
Physicochemical property-based features
Encoding based on grouped weight
Sequence-based features
Encoding based on grouped weight (EBGW) divides 20 amino
Quasi-sequence order
acids into 7 categories based on their hydrophobicity and
Quasi-sequence order (QSO) describes the physicochemical
charge characteristics [92,93]. For each group Hi (i = 1, 2,
distance between amino acids [92]. Most physicochemical
3), a 25-dimensional array Si (i = 1, 2, 3) of the same element
properties are hydrophobicity, hydrophilicity, polarity, and
in the segment should be generated. If the amino acid at that
side-chain volume. This feature was originally proposed by
position belongs to the Hi group, the element in the array will
Chou et al. in [101]. For more detail, refer to [101,106].
be set to 1; otherwise, it will be set to 0. Each array will be
divided into sub-arrays (J-ones), which are represented as
Numerical representation for amino acids
D(j). This value can be taken from cutting the main Si from
It converts each character of amino acids into numerical num-
the first window with len(D(j)) defined as Equation (1):

 bers by mapping them in alphabetic order from 1 to 20, and
jL
lenðDðjÞÞ ¼ int the dummy amino acid X represents 21 [92].
J;
j ¼ 1; 2;    ; J; L ¼ length of segments ð1Þ Binary encoding of amino acids
Binary encoding of amino acids (BINA) represents each
For each group of Hi, a vector with a length of J based on
amino acid as 21-dimensional binary vectors, which encodes
its sub-arrays should be defined, in which the j-th element of
ðjÞ
1 for the target amino acid and 0 for the residues (other
Xi , is calculated based on Equation (2): 20 amino acids). For example, alanine (‘‘A”) is shown as
ðjÞ SumðDðjÞÞ 100000000000000000000 [92].
Xi ¼ ð2Þ
lenðDðjÞÞ
Logo
This feature is defined by calculating the occurrence of amino
Amino acid index acid frequencies and encoding them in a sequence with the
The features based on amino acid index (AAINDEX) are Two Sample Logo program [92].
extracted from the AAINDEX database. This database is used
to represent various physicochemical and biochemical proper- Position weight amino acid composition
ties of each amino acid alone and also in pairs of them in every Position information of each amino acid is another key point
PTM [94]. The feature encodes 14 properties: hydrophobicity, that shall be considered in feature extraction. Position weight
polarity, polarizability, solvent/hydration potential, accessibil- amino acid composition (PWAA) can reveal sequence order
ity reduction ratio, net charge index of side chains, molecular information around P, S, and Y residues [107]. PWAA can
weight, ionization equilibrium constant pKa (-COOH), ioniza- be declared from Equation (3), in which L represents the num-
tion equilibrium constant pKa (-NH3), melting point, optical ber of upstream or downstream amino acids from p-sites in
rotation, entropy of formation, heat capacity, and absolute specific windows. If xi;j = 1, it means that each amino acid
entropy [92,95]. belongs to the j-th position in the window, otherwise xi;j = 0.
1276 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

X

1 L
jjj Shannon entropy
Ci ¼ xi;j j þ Shannon entropy (H) in information theory quantifies the
LðL þ 1Þ j¼
L L
amount of uncertainty of a random variable. To be more pre-
j ¼
L;    ; L ð3Þ cise, it is the average (expected value) amount of information
obtained from observing a random variable. It means that
Composition of k-spaced amino acid pairs when the entropy of a random variable is high, we have more
The encoding of the composition of k-spaced amino acid pairs ambiguity about that random variable [115]. In science and
(CKSAAP) is pretty easy and can be directly calculated from engineering in general, entropy is a measure of the degree of
the sequence pieces of p-sites and non-p-sites. CKSAAP is ambiguity or disorder [116].
one of the important feature encoding schemes in lots of pre-
diction tasks, especially in representing short sequence residues Relative entropy
in protein sequences or subsequences. All 21 amino acids con- Relative entropy (RE) is known as Kullback–Leibler which is
tain 441 different possible pairs. For scanning pieces to count aggregated entropy for more than 20 sites in proteins [74].
all pairs of amino acids with k-space, we can use different win-
dow sizes. For example, window AXXV is a 2-space amino Information gain
acid pair in k = 2 [73,108]. The CKSAAP equation is pro- Information gain (IG) can be computed by subtracting RE
posed as Equation (4) [73]. In this equation, L denotes the from H [Equation (6)] [98].
length of the window, and Ai Aj is an amino acid pair.
  IG ¼ H
RE ð6Þ
Num Ai Aj
fi;j ¼
L
K
1
Accessible surface area
i; j ¼ 1; 2;    ; 21 ð4Þ Accessible surface area (ASA) or solvent-ASA is a biomolecule
surface that can access the solvent. This is an essential struc-
Amino acid composition tural feature determining the protein’s folding and stability
Amino acid composition (AAC) is the most commonly used [98].
feature, which simply calculates the frequency of each amino
acid in subsequences of a protein while encoding the informa- Conventional ML approach
tion into 20 bits [109]. This feature is also represented as amino
acid frequency (AF) in some research. Both AF and AAC
Once the features have been extracted, classification models
reflect the frequency of each amino acid or amino acid pair’s
should be adopted to predict the p-sites. One of the most pop-
occurrence. Lin et al. [109] proposed the AAC equation as
ular classifiers is SVM [97,109,117].
Equation (5), in which ci is the number of amino acid i in
SVM is a linear model for classification and regression
the sequence and vi refers to AAC.
problems that uses a line or hyperplane to separate data. In
ci
vi ¼ other words, SVMs calculate the maximum margin boundary
lenðseqÞ
that leads to the equivalent division of all data points. First,
i ¼ 1;    ; 20 ð5Þ SVM uses a line to classify each data point based on their dis-
tance. If data points are not linearly separable in low-
Evolutionary-based features dimensional space, there may be multiple transformations
KNN enabling the data to be linearly separable in higher dimensions.
The most popular feature selection method that is used in var- Therefore, SVMs can find a hyperplane in higher dimensions
ious ML problems, especially in PTM and phosphorylation between different classes of data such that the distance
classification, is KNN. It classifies sequences based on their between data points falling on either side of that hyperplane
distance. The algorithm classifies sequences by looking at k is maximized [118,119]. Nowadays, SVMs have been widely
of nearest neighbor sequences and finding out the majority used in bioinformatics, especially in PTM problems
of votes from nearest neighbors that have similar attributes [97,107,120].
and the shortest distance as those used to map the items [110]. RF is another well-known and important classifier fre-
quently used in this field. This algorithm can randomly build
Position-specific scoring matrix-based transformation a forest that contains a large number of decision trees. Each
Position-specific scoring matrix-based transformation (PSSM) tree constructs a class prediction, and the class with the most
encodes the evolutionary data of a protein, which is very infor- votes will become the model prediction [121]. Figure 7 demon-
mative and useful for some biological classification problems. strates the procedure of feature extraction, and Figure 8 shows
The PSSM matrix in a protein with a sequence of length L is the process of conventional ML methods.
a matrix with L  20 dimensions. In the matrix, each row rep- As of recent years, kinase-specific methods have been used
resents an amino acid in the protein sequence, and the columns since, in general, some protein prediction sites have not yet
represent the 20 amino acids in proteins [111]. been explored and kinases can assist in locating these sites.
NetPhos [122] and NetPhosK [123] both used deep neural net-
Structural-based features works (DNNs) based on consensus sequences and MS experi-
Protein disorder features mental methods. These algorithms are specific to the kinase’s
All PTMs include p-sites located within disorder positions family. In the Quokka framework [30], the LR approach was
[112]. Protein disorder features (DFs) were used as features suggested to classify 43 S/T and 22 Y kinase family sites.
in many studies [97,113,114]. Kim et al. [117] proposed to use the consensus sequence struc-
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1277

Figure 7 A common procedure of feature extraction

Figure 8 Conventional ML procedure

ML, machine learning; ANN, artificial neural network; SVM, support vector machine; RF, random forest.

ture as features and a SVM classifier to predict four kinase Furthermore, Lin et al. [109] used KNN, AF, and
groups and families. The best ACCs achieved by their model CKSAAP as features and combined different features together
were reported around 83%–95% at the kinase family level to feed it into their model to investigate the best features. The
and 76%–91% at the kinase group level. Liu et al. [124] pro- combination of AF and CKSAAP provided the best ACC for
posed a method for prediction of four kinase families based their SVM model. They believed SVM could classify rice pro-
on RF, which extracted features with an auto covariance teins as universal p-sites. Their work was named Rice_Phos-
(AC) transform and seven physicochemical properties and pho 1.0, which achieved 82% ACC.
achieved over 90% ACC. Cheng et al. [97] proposed a granular SVM (GSVM) for
To recognize protein p-sites in universal proteins, Huang predicting universal p-sites. They used KNN, AF, and DF fea-
et al. [107] proposed a method based on SVM in viruses. They tures in every p-site position to make the train set. To split data
used EBAG and PWAA features for extracting the physico- into high-dimensional feature spaces, they used kernel fuzzy C-
chemical and sequence information of viral proteins around means clustering as a feature extraction method. The method
p-sites. They used 10-fold cross validation and an independent was applied to plant and animal dataset types and could
test set for different window sizes ranging from 15 to 27 amino achieve 80% and 85% ACC scores, respectively.
acids. They got the best results for window size of 23 amino By using the PhosPred-RF method, Banerjee et al. [71] used
acids with ACC scores of 88.8%, 95.2%, and 97.1% for information extracted from PSSM and trained individuals with
S, T, and Y sites, respectively. They also showed the influence RF with odd window sizes ranging from 9 to 25 amino acids.
of using different features. Their model improved almost 15% They got approximately 70% ACC for 26 protein sequences.
when they used the combination of two EBAG and PWAA RF-Phos 1.0 transformed each amino acid into vectors by
features. using eight algorithms of feature selection (H, RE, ASA, OP,
1278 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

ACC, QSO, and the sequence order coupling number of each End-to-end DL approach
sequence) based on a window size of 9 amino acids. They
specifically showed which features are the most important End-to-end learning has become a hot topic in the ML field by
and have more effects on ACC. It was mentioned that AAC taking advantage of DL. DNN is almost the same as tradi-
was the best feature for S and T sites. Then, these features were tional artificial neural networks (ANNs), which is composed
used as RF input with 10-fold cross validation. The ACC of of many connected neurons that work together to solve specific
the model was approximately 80% for S, T, and Y sites [74]. issues, inspired by the functionality of biological neural net-
Moreover, in the RF-Phos 2.0, their RF model was improved works in the human brain. Inspired by the human brain, each
by using window sizes of 5 to 21 amino acids and using differ- DNN’s layer (or group of layers) could be used for learning the
ent features. QSO was the best feature for S and T sites [98]. hierarchical abstraction for downstream tasks. In other words,
RF-Phos 1.0 and RF-Phos 2.0 specifically predicted universal usually raw input sequences are just fed to a DNN, and the
p-sites. It should be mentioned that feature selection methods process of feature selection automatically happens between
helped to improve the ACC of various approaches. layers. Since it refers to training a possibly complex learning
Microbial Phosphorylation Site predictor (MPsites) was system by applying gradient-based learning to the system as
proposed by Hasan et al. [125] to recognize universal microbial a whole, it is called end-to-end DL. These systems are specially
p-sites with different sequence features. In order to convert designed so that all components are created to be differen-
each sequence to numerical vectors, they used various tiable, and consequently, learnable. That is to say, it is a pro-
sequence encoding strategies, including AF, BINA, AAIN- cedure in which a model learns all the steps, including feature
DEX, and PWAA. They used naı̈ve bayes, SVM, neural net- selection and extraction, between the first and last layers [128].
works, decision trees, and RF algorithms to recognize S and Figure 9 shows the common procedure of end-to-end DL
T p-sites. Results showed that RF has better performance than methods.
the other algorithms. It got 68% ACC for S sites and 75% Need to mention that in order to prepare a sequence of
ACC for T sites [125]. amino acids for the end-to-end DL system, there are two pre-
Cao et al. [126] proposed a method named PreSSFP to pre- requisite steps [129]: (1) sequence encoding, and (2) converting
dict p-sites in seven species-specific fungi proteins. They used a the encoded sequence to numerical vectors. The second step
strategy including two steps for feature optimization to could be done via either one-hot encoding [72] or another pop-
improve the SVM prediction performance. KNN, AAC, di- ular technique named word embedding [130]. Therefore, one-
amino acid composition, and physiochemical properties were hot encoding is not considered a feature extraction method
used as features. First, with the RF model, they sorted each and is simply used to represent categorical inputs (e.g., amino
input feature based on the mean ACC. In the second step, acid codes) into numerical vectors in order to feed to DL mod-
the top ten features from the previous step were merged to els. However, in an end-to-end DL network, word embedding
train the SVM model. Finally, they achieved over 80% ACC. is often used in PTM due to the similarity between PTM and
Chen et al. [127] proposed a feature selection method natural language processing (NLP) domains as well as the
named ga-aided ant colony system (GAS), based on ant colony effectiveness of the technique.
and genetic algorithms, for the classification of six kinase We have shown a great success of DL in solving problems
types. in different domains of science, especially in biological prob-
Qiu et al. [105] developed an approach called iPhos-PseEvo. lems with finding non-obvious patterns or making predictions
Protein sequence evolutionary and PseAAC were selected as in datasets [130–136]. In recent years, DL has been applied to
features for an ensemble RF model. The ACC for their model PTM classification of proteins, such as p-site prediction. As
was 71% with the jackknife test evaluation approach. mentioned earlier, the main aspect of this approach compared
As a final example in this section, multi-iPPseEvo [104] is with the conventional ML approach is that the feature extrac-
similar to iPhos-PseEvo but with a different implementation tion step is not designed by human engineers or manually.
strategy while using k-fold cross validation. This method con- These layers are acquired from input data to extract the best
tains a multi-ensemble RF classifier for each S, T, and Y site patterns accurately and quickly. Though the most important
and proposes multi-label p-site prediction for each site. point about DL is that it needs huge amounts of data, by

Figure 9 End-to-end DL procedure

DL, deep learning.
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1279

increasing the size of the dataset, it can perform better. This performance interpretable deep tabular learning network
can be counted as a drawback; when the dataset is not big (TabNet) provides an extremely powerful framework for solv-
enough, it quickly falls behind other ML methods in terms ing more challenging learning problems [144]. For example,
of performance. Khalili et al. [76] developed a TabNet model to predict p-
Among all DL architectures, convolutional neural net- sites in soybean with a high ACC rate that outperformed other
works (CNNs), recurrent neural networks (RNNs), and long common ML methods (LR-L1, LR-L2, RF, SVM, and
short-term memory (LSTM) are the most famous [70,72,137]. XGBoost). They assessed and compared the strength and reli-
Wang et al. [72] provided a DL architecture called Musite- ability of all models using 10-fold cross validation. Experi-
Deep to predict general and kinase-specific families’ positions ments assessed the performance of AAC, dipeptide
in a sequence. The window size used for the input sequence composition (DPC), tripeptide composition (TPC), PSSM,
was 33 amino acids. Then, they presented their network with and physicochemical properties as individual features. To
a multi-layer CNN and attention layer architecture. In con- extract training sequences for model development, various
trast to multi-layer models of MusiteDeep, DeepPhos [70] used window sizes ranging from 7 to 35 amino acids were used.
dense CNN blocks that could show different and multiple rep- They got the best results for window size of 13 amino acids
resentations of proteins for p-site predictions by using the con- with an ACC of 87.34% based on PSSM features.
catenation of intra-block layers and inter-block layers. The Naseer et al. [145] compared human-based feature represen-
method could improve the performance of MusiteDeep by tation with DL-based representation for the reorganization of
using different window sizes with lengths of 15, 33, and 51 phosphoserine p-sites. The combination of the RNN–LSTM
amino acids. Both of these methods, DeepPhos and Musite- model got 81.1% ACC, and the CNN-based model achieved
Deep, have been developed for kinase families and universal 78.3% ACC. In contrast to human engineering with 77%
p-sites. Moreover, PhosTransfer [137] is a DL-based frame- ACC, DL methods have performed better for phosphoserine
work that constructs a pre-train architecture with CNNs based p-site prediction.
on hierarchy kinases systems and transfer learning. It was spe- Even though most DL approaches worked well with large
cialized for improving kinase p-site prediction. The method volumes of data, a study [26] with a small amount of data from
was to accumulate the information of a hierarchical kinase’s only two kinase families proposed a simple DNN architecture
classification tree at family, subfamily, and group levels. It and achieved around 80% ACC. It means that end-to-end
could achieve AUROC of 0.89 on average. The DeepPPSite learning can also perform successfully in low-data regions.
is another DL model based on universal p-site prediction with This algorithm was designed for both the kinase families and
consideration of sequence information [73]. Ahmed et al. used universal p-sites.
one-hot encoding sequence as input, PSSM, EBGW, Guo et al. [146] collected phosphoprotein-binding domains
CKSAAP, and AAINDEX as features, and stacked LSTM (PPBDs) that interact with PPBD-containing proteins (PPCPs)
architecture as a predictive model. The MCC values reported from 12 eukaryotic species and developed a DNN framework
for S, T, and Y are 0.358, 0.356, and 0.350, respectively. based on transfer learning to classify the protein binding
Awais et al. [32] developed a computational model named domains into a hierarchical structure with three levels, includ-
iPhosH-PseAAC using an ANN algorithm to predict PhosH ing group, family, and single PPBD cluster.
sites in protein sequences. The model was based on features Despite most end-to-end approaches using raw sequences
such as PseAAC, statistical moments, and position-relative (one-hot encoding) as input, PhosIDN [147] trained a DNN
features. To validate the iPhosH-PseAAC predictor, they per- by combining raw sequences and PPI information together.
formed self-consistency testing, 10-fold cross validation, and This architecture contains three sub-networks: (1) sequence
jackknife test, which resulted in ACC scores of 100%, feature encoding sub-network (SFENet), (2) PPI feature
94.26%, and 97.07%, respectively. Wang et al. [93] proposed encoding sub-network (IFENet), and (3) heterogeneous fea-
a hybrid model named MaloPred for predicting PhosH sites ture combination sub-network (HFCNet).
in the proteome. The model was composed of two CNN-
based classifiers and a RF-based classifier and was trained p-site prediction tools
on three types of features: one-of-K coding, enhanced grouped
amino acid content (EGAAC), and composition of k-spaced Due to the high cost and low speed of using experimental
amino acid group pairs (CKSAAGP) encoding. They found methods to recognize p-sites, in recent years, many computa-
that MaloPred was able to accurately predict PhosH sites from tional online tools have been developed to help increase the
sequence information through both 10-fold cross validation quality of p-site prediction. Table 2 introduces famous publicly
and independent tests. accessible online tools or GitHub repositories for p-site
Furthermore, there have been some researches such as the prediction.
work of Lv et al. [138] that used hybrid architectures. They
presented a specific hybrid end-to-end architecture that com-
bined both CNN and LSTM together, called DeepIPs, to pre-
Current limitations
dict universal p-sites in host cells infected with SARS-CoV-2
[139,140]. Lv et al. utilized three approaches in NLP as word In general, it is unfair to compare different ML algorithms
embedding layers to represent amino acids as vectors: GloVe applied to p-site prediction tasks to choose the best technique
[141], FastText [142,143], and Word2vec [129] pre-training due to variation in preprocessing steps, evaluation methods,
word embedding methods. The final ACC for this method and, more importantly, database diversity in the literature.
was reported as 80.45% for S/T and 75.22% for Y. Therefore, we tried to evaluate several tools by creating three
DL provides a highly effective framework for dealing with new test datasets. For this purpose, we selected the newly
modern-day learning challenges. The modern high- released version of the dbPTM [66] database in 2022 and
Table 2 Summary of p-site prediction tools

1280
Tool Type/description Method Feature extraction method Dataset size Window size (amino acid) Negative dataset Unbalance strategy Redundancy threshold Evaluation strategy URL U/K

NetPhos [122] Conventional ANN Sequence composition features 902 p-sites 21 (Y, S) – – – 5-fold [Link] K
25 (T)
Kim et al. [117] Conventional SVM – 855 p-sites on S, 3–25 Phosphoproteins Downsampling 70% 7-fold [Link] K
216 p-sites on T
Liu et al. [124] Conventional RF Auto covariance transform, 1911 p-sites – Phosphoproteins Downsampling 40% 5-fold, independent test – K
7 physicochemical properties
Huang et al. [107] Conventional SVM EBAG, PWAA 230 p-sites on S, 23 Phosphoproteins Downsampling – 10-fold, independent test – U
61 p-sites on T,
14 p-sites on Y
Rice_Phospho 1.0 Conventional RF AF, CKSAAP, KNN 4220 p-sites on S, 25 Phosphoproteins Downsampling – 10-fold, independent test [Link] U
[109] 605 p-sites on T,
141 p-sites on Y
GSVM [97] Conventional SVM KNN, AF, DF  50,000 p-sites 13 Phosphoproteins Downsampling 30% 10-fold – U
RF-Phos 1.0 [74] Conventional RF H, RE, ASA, OP, AAC, QSO  28,000 p-sites 5–21 Phosphoproteins Downsampling 30% 10-fold, independent test – U
RF-Phos 2.0 [98] Conventional RF H, RE, IG, ASA, OP, AAC, QSO  28,000 p-sites 5–21 Phosphoproteins Downsampling 30% 10-fold, independent test [Link] Phos/ U
PhosTransfer [137] Conventional CNN H, RE, DF, OP  10,000 p-sites on S, – – Downsampling 40% Independent test [Link] K
 34,000 p-sites on T,
 3000 p-sites on Y

Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

DeepPPsite [73] Conventional LSTM CKSAAP, EBGW, PSSM  7000 p-sites on S, 15, 19, 21 Phosphoproteins Downsampling 30% 10-fold, independent test [Link] U
 2000 p-sites on T,
 700 p-sites on Y
GPS 5.0 [31] Conventional LR Structural features  15,000 p-sites 20 Phosphoproteins Downsampling – 10-fold [Link] K
MPSite [125] Conventional RF AF, OP, PSSM, PWAA  2700 p-sites on S, 7–25 Phosphoproteins Downsampling 30% 10-fold, independent test [Link] U
 2100 p-sites on T
Quokka [30] Conventional LR KNN, AF, BLOUSM64  2400 p-sites on S, 15, 19, 21 Phosphoproteins Downsampling 30% 5-fold, independent test [Link] K
 370 p-sites on T
PhosContext2vec [148] Conventional SVM H, BLOUSM64, DF, OP, ACH, Universal: 25 Phosphoproteins Downsampling – 10-fold, independent test [Link] K/U
secondary structure  20,000 p-sites on S,
 5600 p-sites on T,
 2100 p-sites on Y
Kinases:
 4100
PhosphoSVM [120] Conventionlal SVM H, RE, secondary structure, DF, ASA,  25,000 p-sites on S, 15, 19, 21 Phosphoproteins Downsampling 30% 10-fold, independent test [Link] U
OP, KNN  7200 p-sites on T,
 2700 p-sites on Y
PhosPredRF [71] Conventional RF H, RE, OP  4300 p-sites on S, 15, 19, 21 Phosphoproteins Downsampling 30% 10-fold, independent test [Link] U
 2700 p-sites on T
PreSSFP [126] Conventional SVM Sequence information, evolutionary Various for organisms – Phosphoproteins Downsampling 30% 10-fold, independent test [Link] U
information, physicochemical properties
GasPhos [127] Conventional Multiple GAS  3400 p-sites – Phosphoproteins Downsampling – 5-fold [Link] K
classifiers
iPhos-PseEvo [105] Conventional Ensemble RF KNN, PseAAC 845 p-sites on S, – Phosphoproteins Downsampling 50% Jackknife test [Link] U
386 p-sites on T,
249 p-sites on Y
Multi-iPPseEvo [104] Conventional RF KNN, PseAAC 845 p-sites on S, – Phosphoproteins Downsampling 50% 5-fold [Link] U
386 p-sites on T,
249 p-sites on Y
DeepIPs [138] End-to-end CNN-LSTM – 5387 p-sites on S/T, 33 Phosphoproteins Downsampling 30% Independent test [Link] U
102 p-sites on Y [Link]
DeepPhos [70] End-to-end CNN – 140,000 p-sites on S/T, 15, 33, 51 Phosphoproteins Downsampling 40% 10-fold, independent test [Link] U/K
27,000 p-sites on Y
MusiteDeep [72] End-to-end CNN + –  35,000 p-sites on S/T, 33 Phosphoproteins Downsampling 50% 5-fold, independent test [Link] U/K
attention  2000 p-sites on Y [Link]
Lumbanraja et al. [26] End-to-end DNN –  1800 p-sites on S, 9 Phosphoproteins Downsampling 20% 10-fold – U/K
 700 p-sites on T,
 200 p-sites on Y
PhosIDN [147] End-to-end SFENet + PPI graph embedding  160,000 p-sites 15, 33, 71 Phosphoproteins Downsampling 40% Independent test [Link] U/K
IFENet +
HFCNet
Khalili et al. [76] Neural network + feature TabNet AAC, DPC, PSSM, physicochemical  4500 p-sites 7–35 Phosphoproteins Downsampling 40% 10-fold [Link] U
properties

iPhosH-PseAAC [32] Neural network + feature ANN PseAAC  1300 histidine p-sites – Phosphoproteins Downsampling – 10-fold, – U/K
jackknife test
PROSPECT [33] End-to-end + CNN + RF One-hot encoding,  1600 histidine p-sites 27 – – 40% 10-fold, [Link] U/K
Conventional EGAAC, CKSAAGP independent test

Note: U stands for universal which includes all types of p-sties; K stands for kinase which includes only kinase-specific p-sites. ANN, artificial neural network; LR, logistic regression; SVM, support
vector machine; RF, random forest; KNN, k-nearest neighbor; CNN, convolutional neural network; LSTM, long short-term memory; EBAG, encoding scheme based on attribute grouping; PWAA,
position weight amino acid composition; AF, amino acid frequency; CKSAAP, composition of k-spaced amino acid pairs; DF, protein disorder feature; H, Shannon entropy; RE, relative entropy;
ASA, solvent accessible surface; OP, overlapping properties; AAC, amino acid composition; QSO, quasi-sequence order; PseAAC, pseudo amino acid composition; PPI, protein–protein interaction; –,
not available.
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1281

Table 3 Online tool evaluation

Test set MusiteDeep [72] PhosIDN [147] NetPhos [122]
161–all 161–humans 100–top 161–all 161–humans 100–top 161–all 161–humans 100–top
TP 168 194 150 249 308 297 – 447 339
FP 1656 745 1044 4356 1140 2597 – 3701 5349
TN 11,378 5927 7781 8678 5532 6228 – 2971 3476
FN 201 517 346 120 403 199 – 264 157
ACC (%) 86.14 82.91 85.09 66.60 79.10 70.00 – 46.30 40.93
Precision 0.09 0.21 0.13 0.05 0.21 0.10 – 0.11 0.06
Recall 0.46 0.27 0.30 0.67 0.43 0.60 – 0.63 0.68
F1 0.15 0.24 0.18 0.10 0.29 0.18 – 0.18 0.11
Specificity 0.87 0.89 0.88 0.67 0.83 0.71 – 0.45 0.39
Note: TP, true positive; FP, false positive; TN, true negative; FN, false negative; ACC, accuracy; –, not available.

picked up all new phosphoproteins in all organisms that did Conclusion

not exist in the previous versions. Subsequently, we built
161-all, 161-humans, and 100-top test sets.
Almost all proteins contain phosphorylation, which is respon-
In the 161-all set, 161 new proteins with p-sites were ran-
domly selected from 161 newly released organisms’ proteins sible for critical functions in the cell. Various diseases can be
caused by disruptions of this modification. The discovery of
(one protein per organism). This test set consists of 13,403
phosphorylation as one of the most important PTMs by
sites, 402 of which are p-sites. The maximum and minimum
length of sequences are 7096 and 49 amino acids, respectively. high-throughput experimental methods is labor-intensive and
time-consuming. Therefore, it is urgent to develop a tool or
In the 161-humans set, 161 proteins with p-sites were ran-
method to automatically predict the p-sites. As we investigated
domly selected from newly released proteins of Homo sapiens.
the literature, there is not a complete review for p-site predic-
This test set consists of 7383 sites, 714 of which are p-sites. The
maximum and minimum length of sequences are 921 and 714 tions based on ML algorithms. Due to the importance of the
issue, this review briefly introduces some popular PTM data-
amino acids, respectively.
bases (including phosphorylation), methods, and online tools
In the 100-top set, 100 new proteins with p-sites were ran-
domly selected from the top 10 organisms that had the biggest for p-site prediction to provide a guide to current research.
In this review, we introduced two important databases:
new protein numbers (ten proteins per organism). This test set
EPSD and dbPTM, while comparing them in terms of p-site
consists of 9321 sites, 507 of which are p-sites. The maximum
and minimum length of sequences are 3498 and 102 amino distribution. Then, we gave a brief overview of protein p-site
prediction by ML techniques, which are mainly divided into
acids, respectively.
classical ML and end-to-end DL methods. In addition to
Next, we tried to evaluate several universal p-site prediction
ML, we slightly discussed algorithmic methods as well. Algo-
tools introduced in Table 2 on these datasets. However, there
were many hurdles in the evaluation stage. Kim et al. [117], rithmic methods have a statistical basis which are slow and
have high time complexity. On the other hand, ML algorithms,
RF-Phos 2.0, PhosPred-RF, PreSSFP [126], iPhos-PseEvo,
which are quite popular these days, have attracted a lot of
and Multi-iPPseEvo were not available. Moreover, Rice_Phos-
pho 1.0 and PhosphoSVM only take one sequence as input in attention in p-site prediction including SVM, LR, and RF.
In conventional methods, SVM has shown better performance,
order to process, and since the process was time-consuming,
although the feature extraction step would obviously have a
we could not evaluate our three test datasets on them. Further-
more, DeepIPs did not have any response to our request. significant impact on the final result. Therefore, this study
introduced 20 important and widely used feature extraction
Finally, we selected three tools — MusiteDeep [72], PhosIDN
methods based on structural, sequence, evolutionary, and
[147], and NetPhos [122] — to evaluate. By the way, since Net-
Phos [122] could not predict sequences with a length more than physicochemical properties. Additionally, CNN- and RNN-
based architectures, known as efficient end-to-end learning
4000 amino acids, and the 161-all test set had proteins more
styles, were introduced, which are able to predict p-sites
than that length, we could not evaluate the performance of
directly from the raw input sequences without any feature
NetPhos on the 161-all test set. Table 3 shows the results.
As shown in Table 3, all three tools performed weakly com- extraction steps.
In the next stage, the evaluation methods for predicting p-
pared with the performances reported in their related studies
site approaches were reported to give the standard metrics for
[72,122,147]. We interpreted from the results that there are
no valid benchmarks for p-site prediction. In other words, each comparison between performances. Finally, in order to
demonstrate the current limitations of p-site prediction meth-
study proposed a method applied to a unique test set to report
ods, we created three test sets and evaluated several available
the results, which makes it difficult to compare different meth-
ods together. Therefore, for fair and precise competition, we online tools. All those methods performed poorly compared
with the performances reported in their related studies
suggest that uniform, comprehensive, unique, and well-
[72,122,147], suggesting the importance of creating uniform
defined test benchmarks for p-site prediction will be prepared
and well-defined benchmarks for p-site prediction.
as a crucial step for future research in this field.
1282 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

[11] Pejaver V, Hsu WL, Xin F, Dunker AK, Uversky VN, Radivojac
Competing interests
P. The structural and functional signatures of proteins that
undergo multiple events of post-translational modification.
The authors have declared no competing interests. Protein Sci 2014;23:1077–93.
[12] Khoury GA, Baliban RC, Floudas CA. Proteome-wide post-
translational modification statistics: frequency analysis and
CRediT authorship contribution statement
curation of the Swiss-Prot database. Sci Rep 2011;1:90.
[13] Strumillo M, Beltrao P. Towards the computational design of
Farzaneh Esmaili: Conceptualization, Investigation, Data protein post-translational regulation. Bioorg Med Chem
curation, Visualization, Writing – original draft. Mahdi Pour- 2015;23:2877–82.
mirzaei: Investigation, Data curation, Formal analysis, Visual- [14] Johnson ES. Protein modification by SUMO. Annu Rev
ization, Software, Writing – review & editing. Shahin Ramazi: Biochem 2004;73:355–82.
Investigation, Visualization, Writing – review & editing, Pro- [15] Ahmad I, Hoessli DC, Qazi WM, Khurshid A, Mehmood A,
Walker-Nasir E, et al. MAPRes: an efficient method to analyze
ject administration. Seyedehsamaneh Shojaeilangari: Writing
protein sequence around post-translational modification sites. J
– review & editing. Elham Yavari: Writing – review & editing. Cell Biochem 2008;104:1220–31.
All authors have read and approved the final manuscript. [16] Nickchi P, Jafari M, Kalantari S. PEIMAN 1.0: post-transla-
tional modification enrichment, integration and matching anal-
ysis. Database 2015;2015:bav037.
Acknowledgments [17] Zhou B, Du Y, Xue Y, Miao G, Wei T, Zhang P. Identification
of malonylation, succinylation, and glutarylation in serum
We gratefully thanks Hadi Pourmirzaei for preparing pictures proteins of acute myocardial infarction patients. Proteomics
and Mohamad Ezati for his helps and recommendation. Clinical Appl 2020;14:e1900103.
[18] Huang H, Arighi CN, Ross KE, Ren J, Li G, Chen SC, et al.
iPTMnet: an integrated resource for protein post-translational
ORCID modification network discovery. Nucleic Acids Res 2018;46:
D542–50.
ORCID 0000-0002-0002-7305 (Farzaneh Esmaili) [19] Kamath KS, Vasavada MS, Srivastava S. Proteomic databases
ORCID 0000-0003-4621-0372 (Mahdi Pourmirzaei) and tools to decipher post-translational modifications. J Pro-
teomics 2011;75:127–44.
ORCID 0000-0002-9043-1140 (Shahin Ramazi)
[20] Karve TM, Cheema AK. Small changes huge impact: the role of
ORCID 0000-0002-7013-3330 (Seyedehsamaneh Shojaeilangari) protein posttranslational modifications in cellular homeostasis
ORCID 0000-0002-3035-0445 (Elham Yavari) and disease. J Amino Acids 2011;2011:207691.
[21] Schedin-Weiss S, Winblad B, Tjernberg LO. The role of protein
References glycosylation in Alzheimer disease. FEBS J 2014;281:46–62.
[22] Falkenberg KJ, Johnstone RW. Histone deacetylases and their
[1] Craveur P, Rebehmed J, de Brevern AG. PTM-SD: a database of inhibitors in cancer, neurological diseases and immune disorders.
structurally resolved and annotated posttranslational modifica- Nat Rev Drug Discov 2014;13:673–91.
tions in proteins. Database 2014;2014:bau041. [23] Park G, Tan J, Garcia G, Kang Y, Salvesen G, Zhang Z.
[2] Sreedhar A, Wiese EK, Hitosugi T. Enzymatic and metabolic Regulation of histone acetylation by autophagy in Parkinson
regulation of lysine succinylation. Genes Dis 2020;7:166–71. disease. J Biol Chem 2016;291:3531–40.
[3] Li J, Jia J, Li H, Yu J, Sun H, He Y, et al. SysPTM 2.0: an [24] Popovic D, Vucic D, Dikic I. Ubiquitination in disease patho-
updated systematic resource for post-translational modification. genesis and treatment. Nat Med 2014;20:1242–53.
Database 2014;2014:bau025. [25] Levene PA, Alsberg CL. The cleavage products of vitellin. J Biol
[4] Audagnotto M, Dal Peraro M. Protein post-translational mod- Chem 1906;2:127–33.
ifications: in silico prediction tools and molecular modeling. [26] Lumbanraja FR, Mahesworo B, Cenggoro TW, Budiarto A,
Comput Struct Biotechnol J 2017;15:307–19. Pardamean B. An evaluation of deep neural network perfor-
[5] Xu Y, Wang X, Wang Y, Tian Y, Shao X, Wu LY, et al. mance on limited protein phosphorylation site prediction data.
Prediction of posttranslational modification sites from amino Procedia Comput Sci 2019;157:25–30.
acid sequences with kernel methods. J Theor Biol [27] Tenreiro S, Eckermann K, Outeiro TF. Protein phosphorylation
2014;344:78–87. in neurodegeneration: friend or foe? Front Mol Neurosci
[6] Farriol-Mathis N, Garavelli JS, Boeckmann B, Duvaud S, 2014;7:42.
Gasteiger E, Gateau A, et al. Annotation of post-translational [28] Barber KW, Rinehart J. The ABCs of PTMs. Nat Chem Biol
modifications in the Swiss-Prot knowledge base. Proteomics 2018;14:188–92.
2004;4:1537–50. [29] Ardito F, Giuliani M, Perrone D, Troiano G, Lo ML. The
[7] Ramazi S, Zahiri J, Arab S, Parandian Y. Computational crucial role of protein phosphorylation in cell signalingand its
prediction of proteins sumoylation: a review on the methods and use as targeted therapy (Review). Int J Mol Med 2017;40:271–80.
databases. J Nanomed Res 2016;3:00068. [30] Li F, Li C, Marquez-Lago TT, Leier A, Akutsu T, Purcell AW,
[8] Xu H, Wang Y, Lin S, Deng W, Peng D, Cui Q, et al. PTMD: a et al. Quokka: a comprehensive tool for rapid and accurate
database of human disease-associated post-translational modifi- prediction of kinase family-specific phosphorylation sites in the
cations. Genomics Proteomics Bioinformatics 2018;16:244–51. human proteome. Bioinformatics 2018;34:4223–31.
[9] Duan G, Walther D. The roles of post-translational modifica- [31] Wang C, Xu H, Lin S, Deng W, Zhou J, Zhang Y, et al. GPS 5.0:
tions in the context of protein interaction networks. PLoS an update on the prediction of kinase-specific phosphorylation
Comput Biol 2015;11:e1004049. sites in proteins. Genomics Proteomics Bioinformatics
[10] Alleyn M, Breitzig M, Lockey R, Kolliputi N. The dawn of 2020;18:72–80.
succinylation: a posttranslational modification. Am J Physiol [32] Awais M, Hussain W, Khan YD, Rasool N, Khan SA, Chou
Cell Physiol 2018;314:C228–32. KC. iPhosH-PseAAC: identify phosphohistidine sites in proteins
by blending statistical moments and position relative features
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1283

according to the Chou’s 5-step rule and general pseudo amino [53] Gnad F, Gunawardena J, Mann M. PHOSIDA 2011: the
acid composition. IEEE/ACM Trans Comput Biol Bioinform posttranslational modification database. Nucleic Acids Res
2021;18:596–610. 2011;39:D253–60.
[33] Chen Z, Zhao P, Li F, Leier A, Marquez-Lago TT, Webb GI, [54] Safaei J, Maňuch J, Gupta A, Stacho L, Pelech S. Prediction of
et al. PROSPECT: a web server for predicting protein histidine 492 human protein kinase substrate specificities. Proteome Sci
phosphorylation sites. J Bioinform Comput Biol 2011;9:S6.
2020;18:2050018. [55] Pan Z, Wang B, Zhang Y, Wang Y, Ullah S, Jian R, et al.
[34] Peng C, Lu Z, Xie Z, Cheng Z, Chen Y, Tan M, et al. The first dbPSP: a curated database for protein phosphorylation sites in
identification of lysine malonylation substrates and its regulatory prokaryotes. Database 2015;2015:bav031.
enzyme. Mol Cell Proteomics 2011;10:M111.012658. [56] Qi L, Liu Z, Wang J, Cui Y, Guo Y, Zhou T, et al. Systematic
[35] Ferguson FM, Gray NS. Kinase inhibitors: the road ahead. Nat analysis of the phosphoproteome and kinase-substrate networks
Rev Drug Discov 2018;17:353–77. in the mouse testis. Mol Cell Proteomics 2014;13:3626–38.
[36] Ramazi S, Zahiri J. Posttranslational modifications in proteins: [57] Yao Q, Ge H, Wu S, Zhang N, Chen W, Xu C, et al. P3DB 3.0:
resources, tools and prediction methods. Database (Oxford) from plant phosphorylation sites to protein networks. Nucleic
2021;2021:baab012. Acids Res 2014;42:D1206–13.
[37] Thapa N, Chaudhari M, Iannetta AA, White C, Roy K, [58] Cheng H, Deng W, Wang Y, Ren J, Liu Z, Xue Y. dbPPT: a
Newman RH, et al. A deep learning based approach for comprehensive database of protein phosphorylation in plants.
prediction of Chlamydomonas reinhardtii phosphorylation sites. Database 2014;2014:bau121.
Sci Rep 2021;11:12550. [59] Ullah S, Lin S, Xu Y, Deng W, Ma L, Zhang Y, et al. dbPAF: an
[38] Li H, Xing X, Ding G, Li Q, Wang C, Xie L, et al. SysPTM: a integrative database of protein phosphorylation in animals and
systematic resource for proteomic research on post-translational fungi. Sci Rep 2016;6:23534.
modifications. Mol Cell Proteomics 2009;8:1839–49. [60] UniProt Consortium. UniProt: a worldwide hub of protein
[39] Newman RH, Zhang J, Zhu H. Toward a systems-level view of knowledge. Nucleic Acids Res 2019;47:D506–15.
dynamic phosphorylation networks. Front Genet 2014;5:263. [61] Bodenmiller B, Malmstrom J, Gerrits B, Campbell D, Lam H,
[40] Shi XX, Wu FX, Mei LC, Wang YL, Hao GF, Yang GF. Schmidt A, et al. PhosphoPep — a phosphoproteome resource
Bioinformatics toolbox for exploring protein phosphorylation for systems biology research in Drosophila Kc167 cells. Mol Syst
network. Brief Bioinform 2021;22:bbaa134. Biol 2007;3:139.
[41] Rashid MM, Shatabda S, Hasan MM, Kurata H. Recent [62] Oughtred R, Rust J, Chang C, Breitkreutz B, Stark C, Willems
development of machine learning methods in microbial phos- A, et al. The BioGRID database: a comprehensive biomedical
phorylation sites. Curr Genomics 2020;21:194–203. resource of curated protein, genetic, and chemical interactions.
[42] Keshava Prasad TS, Goel R, Kandasamy K, Keerthikumar S, Protein Sci 2021;30:187–200.
Kumar S, Mathivanan S, et al. Human protein reference [63] Bai Y, Chen B, Li M, Zhou Y, Ren S, Xu Q, et al. FPD: a
database — 2009 update. Nucleic Acids Res 2009;37:D767–72. comprehensive phosphorylation database in fungi. Fungal Biol
[43] Huang KY, Lee TY, Kao HJ, Ma CT, Lee CC, Lin TH, et al. 2017;121:869–75.
dbPTM in 2019: exploring disease association and cross-talk of [64] de Bruijn FJ. Medicago truncatula proteomics: introduction. In:
post-translational modifications. Nucleic Acids Res 2019;47: de Bruijn FJ, editor. Model legume Medicago truncatula. New
D298–308. York: John Wiley & Sons, Inc; 2020, p. 1069.
[44] Boeckmann B, Bairoch A, Apweiler R, Blatter MC, Estreicher [65] Durek P, Schmidt R, Heazlewood JL, Jones A, MacLean D,
A, Gasteiger E, et al. The Swiss-Prot protein knowledgebase and Nagel A, et al. PhosPhAt: the Arabidopsis thaliana phosphory-
its supplement TrEMBL in 2003. Nucleic Acids Res lation site database. An update. Nucleic Acids Res 2010;38:
2003;31:365–70. D828–34.
[45] Lin S, Wang C, Zhou J, Shi Y, Ruan C, Tu Y, et al. EPSD: a [66] Li Z, Li S, Luo M, Jhong JH, Li W, Yao L, et al. dbPTM in
well-annotated data resource of protein phosphorylation sites in 2022: an updated database for exploring regulatory networks
eukaryotes. Brief Bioinform 2021;22:298–307. and functional associations of protein post-translational modi-
[46] Nguyen TD, Vidal-Cortes O, Gallardo O, Abian J, Carrascal M. fications. Nucleic Acids Res 2022;50:D471–9.
LymPHOS 2.0: an update of a phosphosite database of primary [67] Wagih O, Reimand J, Bader GD. MIMP: predicting the impact
human T cells. Database 2015;2015:bav115. of mutations on kinase-substrate phosphorylation. Nat Methods
[47] Zanzoni A, Ausiello G, Via A, Gherardini PF, Helmer-Citterich 2015;12:531–3.
M. Phospho3D: a database of three-dimensional structures of [68] Krassowski M, Pellegrina D, Mee MW, Fradet-Turcotte A, Bhat
protein phosphorylation sites. Nucleic Acids Res 2007;35: M, Reimand J. ActiveDriverDB: interpreting genetic variation in
D229–31. human and cancer genomes using post-translational modifica-
[48] Dinkel H, Chica C, Via A, Gould CM, Jensen LJ, Gibson TJ, tion sites and signaling networks (2021 update). Front Cell Dev
et al. [Link]: a database of phosphorylation sites — Biol 2021;9:626821.
update 2011. Nucleic Acids Res 2011;39:D261–7. [69] Peng D, Li H, Hu B, Zhang H, Chen L, Lin S, et al. PTMsnp: a
[49] Huang KY, Wu HY, Chen YJ, Lu CT, Su MG, Hsieh YC, et al. web server for the identification of driver mutations that affect
RegPhos 2.0: an updated resource to explore protein kinase- protein post-translational modification. Front Cell Dev Biol
substrate phosphorylation networks in mammals. Database 2020;8:593661.
2014;2014:bau034. [70] Luo F, Wang M, Liu Y, Zhao XM, Li A. DeepPhos: prediction
[50] Hornbeck PV, Zhang B, Murray B, Kornhauser JM, Latham V, of protein phosphorylation sites with deep learning. Bioinfor-
Skrzypek E. PhosphoSitePlus, 2014: mutations, PTMs and matics 2019;35:2766–73.
recalibrations. Nucleic Acids Res 2015;43:D512–20. [71] Banerjee S, Basu S, Ghosh D, Nasipuri M. PhosPred-RF:
[51] Minguez P, Letunic I, Parca L, Garcia-Alonso L, Dopazo J, prediction of protein phosphorylation sites using a consensus of
Huerta-Cepas J, et al. PTMcode v2: a resource for functional random forest classifiers. Int Conf Work Comput Commun
associations of post-translational modifications within and 2015;2015:1–7.
between proteins. Nucleic Acids Res 2015;43:D494–502. [72] Wang D, Zeng S, Xu C, Qiu W, Liang Y, Joshi T, et al.
[52] Yu K, Wang Y, Zheng Y, Liu Z, Zhang Q, Wang S, et al. qPTM: MusiteDeep: a deep-learning framework for general and kinase-
an updated database for PTM dynamics in human, mouse, rat specific phosphorylation site prediction. Bioinformatics
and yeast. Nucleic Acids Res 2023;51:D479–87. 2017;33:3909–16.
1284 Genomics Proteomics Bioinformatics 21 (2023) 1266–1285

[73] Ahmed S, Kabir M, Arif M, Khan ZU, Yu DJ. DeepPPSite: a [94] Kawashima S, Pokarowski P, Pokarowska M, Kolinski A,
deep learning-based model for analysis and prediction of Katayama T, Kanehisa M. AAindex: amino acid index database,
phosphorylation sites using efficient sequence information. Anal progress report 2008. Nucleic Acids Res 2008;36:D202–5.
Biochem 2021;612:113955. [95] Xu Y, Ding YX, Ding J, Wu LY, Xue Y. Mal-Lys: prediction of
[74] Ismail HD, Jones A, Kim JH, Newman RH, Dukka BKC. lysine malonylation sites in proteins integrated sequence-based
Phosphorylation sites prediction using random forest. IEEE 5th features with mRMR feature selection. Sci Rep 2016;6:38318.
Int Conf Comput Adv Bio Med Sci 2015:1–6. [96] Zhang T, Zhang H, Chen K, Shen S, Ruan J, Kurgan L.
[75] Li W, Godzik A. Cd-hit: a fast program for clustering and Accurate sequence-based prediction of catalytic residues. Bioin-
comparing large sets of protein or nucleotide sequences. Bioin- formatics 2008;24:2329–38.
formatics 2006;22:1658–9. [97] Cheng G, Chen Q, Zhang R. Prediction of phosphorylation sites
[76] Khalili E, Ramazi S, Ghanati F, Kouchaki S. Predicting protein based on granular support vector machine. Granul Comput
phosphorylation sites in soybean using interpretable deep tabular 2021;6:107–17.
learning network. Brief Bioinform 2022;23:bbac015. [98] Ismail HD, Jones A, Kim JH, Newman RH, Kc DB. RF-Phos: a
[77] Huang C, Li Y, Loy CC, Tang X. Learning deep representation novel general phosphorylation site prediction tool based on
for imbalanced classification. Proc IEEE Conf Comput Vis random forest. Biomed Res Int 2016;2016:3281590.
Pattern Recognit 2016:5375–84. [99] Fan SC, Zhang XG. Characterizing the microenvironment
[78] Chou KC, Zhang CT. Prediction of protein structural classes. surrounding phosphorylated protein sites. Genomics Proteomics
Crit Rev Biochem Mol Biol 1995;30:275–349. Bioinformatics 2005;3:213–7.
[79] Chen Z, Liu X, Li F, Li C, Marquez-Lago T, Leier A, et al. [100] Zhang ZH, Wang ZH, Zhang ZR, Wang YX. A novel method
Large-scale comparative assessment of computational predictors for apoptosis protein subcellular localization prediction combin-
for lysine post-translational modification sites. Brief Bioinform ing encoding based on grouped weight and support vector
2019;20:2267–90. machine. FEBS Lett 2006;580:6169–74.
[80] Schwartz D, Gygi SP. An iterative statistical approach to the [101] Chou KC. Prediction of protein subcellular locations by incor-
identification of protein phosphorylation motifs from large-scale porating quasi-sequence-order effect. Biochem Biophys Res
data sets. Nat Biotechnol 2005;23:1391–8. Commun 2000;278:477–83.
[81] Chen YC, Aguan K, Yang CW, Wang YT, Pal NR, Chung IF. [102] Xiang Q, Feng K, Liao B, Liu Y, Huang G. Prediction of lysine
Discovery of protein phosphorylation motifs through explora- malonylation sites based on pseudo amino acid. Comb Chem
tory data analysis. PLoS One 2011;6:e20025. High Throughput Screen 2017;20:622–8.
[82] He Z, Yang C, Guo G, Li N, Yu W. Motif-All: discovering all [103] Liu LM, Xu Y, Chou KC. iPGK-PseAAC: identify lysine
phosphorylation motifs. BMC Bioinformatics 2011;12:S22. phosphoglycerylation sites in proteins by incorporating four
[83] Zhou FF, Xue Y, Chen GL, Yao X. GPS: a novel group-based different tiers of amino acid pairwise coupling information into
phosphorylation predicting and scoring method. Biochem Bio- the general PseAAC. Med Chem 2017;13:552–9.
phys Res Commun 2004;325:1443–8. [104] Qiu W, Zheng Q, Sun B, Xiao X. Multi-iPPseEvo: a multi-label
[84] Xue Y, Zhou F, Zhu M, Ahmed K, Chen G, Yao X. GPS: a classifier for identifying human phosphorylated proteins by
comprehensive www server for phosphorylation sites prediction. incorporating evolutionary information into Chou’s general
Nucleic Acids Res 2005;33:W184–7. PseAAC via grey system theory. Mol Inform 2017;36:1600085.
[85] Xue Y, Ren J, Gao X, Jin C, Wen L, Yao X. GPS 2.0, a tool to [105] Qiu W, Sun B, Xiao X, Xu D, Chou K. iPhos-PseEvo:
predict kinase-specific phosphorylation sites in hierarchy. Mol identifying human phosphorylated proteins by incorporating
Cell Proteomics 2008;7:1598–608. evolutionary information into general PseAAC via grey system
[86] Xue Y, Liu Z, Cao J, Ma Q, Gao X, Wang Q, et al. GPS 2.1: theory. Mol Inform 2017;36:1600010.
enhanced prediction of kinase-specific phosphorylation sites with [106] Wang J, Yang B, Leier A, Marquez-Lago TT, Hayashida M,
an algorithm of motif length selection. Protein Eng Des Sel Rocker A, et al. Bastion6: a bioinformatics approach for
2011;24:255–60. accurate prediction of type VI secreted effectors. Bioinformatics
[87] Liu Z, Yuan F, Ren J, Cao J, Zhou Y, Yang Q, et al. GPS-ARM: 2018;34:2546–55.
computational analysis of the APC/C recognition motif by [107] Huang SY, Shi SP, Qiu JD, Liu MC. Using support vector
predicting D-boxes and KEN-boxes. PLoS One 2012;7:e34370. machines to identify protein phosphorylation sites in viruses. J
[88] Deng W, Wang Y, Ma L, Zhang Y, Ullah S, Xue Y. Mol Graph Model 2015;56:84–90.
Computational prediction of methylation types of covalently [108] Chen Z, Chen YZ, Wang XF, Wang C, Yan RX, Zhang Z.
modified lysine and arginine residues in proteins. Brief Bioinform Prediction of ubiquitination sites by using the composition of k-
2017;18:647–58. spaced amino acid pairs. PLoS One 2011;6:e22930.
[89] Zhao Q, Xie Y, Zheng Y, Jiang S, Liu W, Mu W, et al. GPS- [109] Lin S, Song Q, Tao H, Wang W, Wan W, Huang J, et al.
SUMO: a tool for the prediction of sumoylation sites Rice_Phospho 1.0: a new rice-specific SVM predictor for protein
and SUMO-interaction motifs. Nucleic Acids Res 2014;42: phosphorylation sites. Sci Rep 2015;5:1–9.
W325–30. [110] Kramer O. K-nearest neighbors. In: Kramer O, editor. Dimen-
[90] Jordan MI, Mitchell TM. Machine learning: trends, perspectives, sionality reduction with unsupervised nearest neighbors. Berlin:
and prospects. Science 2015;349:255–60. Springer; 2013, p.13–23.
[91] Jamal S, Ali W, Nagpal P, Grover A, Grover S. Predicting [111] Wang J, Yang B, Revote J, Leier A, Marquez-Lago TT, Webb
phosphorylation sites using machine learning by integrating the G, et al. POSSUM: a bioinformatics toolkit for generating
sequence, structure, and functional information of proteins. J numerical sequence feature descriptors based on PSSM profiles.
Transl Med 2021;19:218. Bioinformatics 2017;33:2756–8.
[92] Zhang Y, Xie R, Wang J, Leier A, Marquez-Lago TT, Akutsu T, [112] Dunker AK, Oldfield CJ, Meng J, Romero P, Yang JY, Chen
et al. Computational analysis and prediction of lysine JW, et al. The unfoldomics decade: an update on intrinsically
malonylation sites by exploiting informative features in an disordered proteins. BMC Genomics 2008;9:S1.
integrative machine-learning framework. Brief Bioinform [113] Iakoucheva LM, Radivojac P, Brown CJ, O’Connor TR, Sikes
2019;20:2185–99. JG, Obradovic Z, et al. The importance of intrinsic disorder for
[93] Wang LN, Shi SP, Xu HD, Wen PP, Qiu JD. Computational protein phosphorylation. Nucleic Acids Res 2004;32:1037–49.
prediction of species-specific malonylation sites via enhanced [114] Obradovic Z, Peng K, Vucetic S, Radivojac P, Dunker AK.
characteristic strategy. Bioinformatics 2017;33:1457–63. Exploiting heterogeneous sequence properties improves predic-
 Esmaili F et al / Review of Machine Learning Methods for Phosphorylation Prediction 1285

tion of protein disorder. Proteins Struct Funct Bioinforma [132] Webb S. Deep learning for biology. Nature 2018;554:555–8.
2005;61:176–82. [133] Wainberg M, Merico D, Delong A, Frey BJ. Deep learning in
[115] Shannon CE. A mathematical theory of communication. Bell biomedicine. Nat Biotechnol 2018;36:829–38.
Syst Tech J 1948;27:379–423,623–656. [134] Angermueller C, Pärnamaa T, Parts L, Stegle O. Deep learning
[116] Capra JA, Singh M. Predicting functionally important residues for computational biology. Mol Syst Biol 2016;12:878.
from sequence conservation. Bioinformatics 2007;23:1875–82. [135] Wei GW. Protein structure prediction beyond AlphaFold. Nat
[117] Kim JH, Lee J, Oh B, Kimm K, Koh I. Prediction of Mach Intell 2019;1:336–7.
phosphorylation sites using SVMs. Bioinformatics [136] Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ron-
2004;20:3179–84. neberger O, et al. Highly accurate protein structure prediction
[118] Drucker H, Burges CJC, Kaufman L, Smola A, Vapnik V. with AlphaFold. Nature 2021;596:583–9.
Support vector regression machines. Proc 9th Int Conf Neural [137] Xu Y, Wilson C, Leier A, Marquez-Lago TT, Whisstock J, Song
Inf Process Syst 1997;9:155–61. J. PhosTransfer: a deep transfer learning framework for kinase-
[119] Noble WS. What is a support vector machine? Nat Biotechnol specific phosphorylation site prediction in hierarchy. In: Lauw
2006;24:1565–7. H, Wong RW, Ntoulas A, Lim EP, Ng SK, Pan S, editors.
[120] Dou Y, Yao B, Zhang C. PhosphoSVM: prediction of phos- Advances in knowledge discovery and data min-
phorylation sites by integrating various protein sequence ing. Cham: Springer; 2020, p.384–95.
attributes with a support vector machine. Amino Acids [138] Lv H, Dao FY, Zulfiqar H, Lin H. DeepIPs: comprehensive
2014;46:1459–69. assessment and computational identification of phosphorylation
[121] Pal M. Random forest classifier for remote sensing classification. sites of SARS-CoV-2 infection using a deep learning-based
Int J Remote Sens 2005;26:217–22. approach. Brief Bioinform 2021;22:bbab244.
[122] Blom N, Gammeltoft S, Brunak S. Sequence and structure-based [139] Barnes CO, Jette CA, Abernathy ME, Dam KMA, Esswein SR,
prediction of eukaryotic protein phosphorylation sites. J Mol Gristick HB, et al. SARS-CoV-2 neutralizing antibody structures
Biol 1999;294:1351–62. inform therapeutic strategies. Nature 2020;588:682–7.
[123] Hjerrild M, Stensballe A, Rasmussen TE, Kofoed CB, Blom N, [140] Hu B, Guo H, Zhou P, Shi ZL. Characteristics of SARS-CoV-2
Sicheritz-Ponten T, et al. Identification of phosphorylation sites and COVID-19. Nat Rev Microbiol 2021;19:141–54.
in protein kinase A substrates using artificial neural networks [141] Pennington J, Socher R, Manning CD. Glove: global vectors for
and mass spectrometry. J Proteome Res 2004;3:426–33. word representation. Proc 2014 Conf Empir Methods Nat Lang
[124] Liu W, Guo Y, Luo J, Zhong Y, Yang X, Pu X, et al. Prediction Process 2014:1532–43.
of kinase-specific phosphorylational interactions using random [142] Bojanowski P, Grave E, Joulin A, Mikolov T. Enriching word
forest. Chemom Intell Lab Syst 2013;126:117–22. vectors with subword information. Trans Assoc Comput Lin-
[125] Hasan MM, Rashid MM, Khatun MS, Kurata H. Computa- guist 2017;5:135–46.
tional identification of microbial phosphorylation sites by the [143] Joulin A, Grave E, Bojanowski P, Douze M, Jégou H, Mikolov
enhanced characteristics of sequence information. Sci Rep T. [Link]: compressing text classification models. arXiv
2019;9:8258. 2016;1612.03651.
[126] Cao M, Chen G, Yu J, Shi S. Computational prediction and [144] Arik SÖ, Pfister T. TabNet: attentive interpretable tabular
analysis of species-specific fungi phosphorylation via feature learning. Proc AAAI Conf Artif Intell 2021:6679–87.
optimization strategy. Brief Bioinform 2020;21:595–608. [145] Naseer S, Hussain W, Khan YD, Rasool N. Optimization of
[127] Chen CW, Huang LY, Liao CF, Chang KP, Chu YW. GasPhos: serine phosphorylation prediction in proteins by comparing
protein phosphorylation site prediction using a new feature human engineered features and deep representations. Anal
selection approach with a GA-aided ant colony system. Int J Mol Biochem 2021;615:114069.
Sci 2020;21:7891. [146] Guo Y, Ning W, Jiang P, Lin S, Wang C, Tan X, et al. GPS-PBS:
[128] Glasmachers T. Limits of end-to-end learning. Asian Conf Mach a deep learning framework to predict phosphorylation sites that
Learn 2017:17–32. specifically interact with phosphoprotein-binding domains. Cells
[129] Mikolov T, Chen K, Corrado G, Dean J. Efficient estimation of 2020;9:1266.
word representations in vector space. arXiv 2013;1301.3781. [147] Yang H, Wang M, Liu X, Zhao XM, Li A. PhosIDN: an
[130] Elnaggar A, Heinzinger M, Dallago C, Rihawi G, Wang Y, integrated deep neural network for improving protein phospho-
Jones L, et al. ProtTrans: towards cracking the language of life’s rylation site prediction by combining sequence and protein–
code through self-supervised deep learning and high performance protein interaction information. Bioinformatics
computing. arXiv 2020;2007.06225. 2021;37:4668–76.
[131] Nambiar A, Heflin M, Liu S, Maslov S, Hopkins M, Ritz A. [148] Xu Y, Song J, Wilson C, Whisstock JC. PhosContext2vec: a
Transforming the language of life: transformer neural networks distributed representation of residue-level sequence contexts and
for protein prediction tasks. Proc 11th ACM Int Conf Bioinfor- its application to general and kinase-specific phosphorylation site
matics Comput Biol Heal Informatics 2020:1–8. prediction. Sci Rep 2018;8:8240.