DeepIEP: a Peptide Sequence Model of Isoelectric Point (IEP/pI) using Recurrent Neural Networks (RNNs)
Esben Jannik Bjerrum
q-bio.BM cs.LG q-bio.QM


[   \begin@twocolumnfalse




DeepIEP: a Peptide Sequence Model of Isoelectric Point (IEP/pI) using
Recurrent Neural Networks (RNNs) 
    J. K. Verma, Harmanpreet Singh and P. K. Pathak
    December 30, 2023
=======================================================================================================

fancy




The isoelectric point (IEP or pI) is the pH where the net charge on the molecular ensemble of peptides and proteins is zero. This physical-chemical property is dependent on protonable/deprotonable sidechains and their pKa values. Here an pI prediction model is trained from a database of peptide sequences and pIs using a recurrent neural network (RNN) with long short-term memory (LSTM) cells. The trained model obtains an RMSE and R^2 of 0.28 and 0.95 for the external test set. The model is not based on pKa values, but prediction of constructed test sequences show similar rankings as already known pKa values. The prediction depends mostly on the existence of known acidic and basic amino acids with fine-adjusted based on the neighboring sequence and position of the charged amino acids in the peptide chain.

 \end@twocolumnfalse ] 




§ INTRODUCTION


The isoelectric point (pI/IEP) of peptides and proteins is the pH
value where the net charge on the molecular ensemble is zero. It is
of practical interest in work involving peptide and proteins with
regard to solubility/precipitation<cit.>, separation for
proteomics and elution in ion exchange chromatography<cit.>,
in peptide drug design<cit.> and for understanding protein
function<cit.>. Good models has long been
available<cit.>, as the pI can be obtained by merely
counting the number of each type of amino acid, and knowing their
pKa values, solve the Henderson-Hasselbalch equation to determine
the pH where the overall charge is zero. The prediction is dependent
on the pKa values chosen for the amino acid side chains and these
are not the same in a peptide or protein as in bulk solution.<cit.>
It has been a popular exercise to fine-tune these pKa-sets to obtain
better pI models, and benchmarks of different pKa sets has been conducted<cit.>.
However, these models do not take sequence dependent and position
effects into account. As examples, it is known that the clustering
of e.g. acidic residues will lower their apparent pKa as their close
proximity to each others charges will affect the deprotonation and
hence calculation of pI values.<cit.> Additionally, it
has been found that pKa values should be adjusted for amino acids
at the C- or N- terminal, as the close proximity to the amino- or
carboxylic acid from the ends of the main chain will affect the deprotonation
and apparent pKa value. Although it is unclear if this is due to proximity
to the terminal charge<cit.> or formation of artifacts
such as pyroglutamate and pyroaspartate<cit.>.
To potentially solve these effects, specific pKa values for the terminal
amino acids have been employed<cit.> as well as more
advanced machine learning modes using SVMs<cit.> and including
descriptors calculated with cheminformatics toolkits<cit.>.

Having previously good experience with employing recurrent neural
networks (RNNs) with long short-term memory cells (LSTM) to make regression
models of serialized small molecules<cit.>, it was obvious
to try and employ the same model architecture to the natural sequences
of amino acids in peptides. Previous work has been done with using
artificial neural networks or recurrent neural networks to model protein
and peptide properties and functions <cit.>,
but it is to my knowledge the first time this has been tried for prediction
of pI. The model is trained directly from a database of sequences
and their corresponding pI value, without any knowledge about pKa
values of amino acids or the Henderson-Hasselbalch equation.




§ METHODS





  §.§.§ Datasets


The dataset was combined from two datasets<cit.>
derived from proteomics and a dataset extracted from Reaxys<cit.>
as previously described<cit.>, using the plain peptide
sequences without modifications. The peptides in the proteomics studies<cit.>
were reacted with iodoacetamide and the cysteine characters C was
substituted with Z. The pIs from the proteomics studies were set to
the average of the gel fractions where the peptide was identified.
After combination the datasets were randomly divided into a 90% train
and 10% test set.




  §.§.§ Vectorization


The amino acid sequences were converted to one-hot encoded vectors,
where each row corresponded to an amino acid. The vectors were expanded
to similar dimensions by padding with ones in a separate row. Before
training the vectors were flipped along the sequence axis, so that
the padding characters were the first and the encoded sequences the
last. The amino acid sequence was thus flipped in the same operation,
and the RNN network read the sequence from the C-terminal to the N-terminal.




  §.§.§ Neural Network Modeling


An RNN model was constructed using Keras v. 2.1.1<cit.>
and Tensorflow v. 1.4<cit.>. The first layer consisted
of 200 LSTM cells<cit.>, with a recurrent dropout
of 0.3 and used in batch mode. The final internal C state was used
as input to a dense layer of 100 neurons with the rectified linear
unit activation function<cit.>. The output from the dense
layer was connected to a single output neuron with a linear activation
function.

The network was trained with mini-batches of 256 sequences for 300
epochs using the mean square error loss function and the Adam optimizer
<cit.> with an initial learning rate of 0.005. The loss
was followed on the test set and the learning rate lowered by a factor
of two when no improvement in the test set loss had been observed
for 15 epochs. Training took approximately a quarter of an hour.

After training in batch mode, a stateful model was constructed by
creating a model with the exact same architecture but the LSTM states
set to stateful and the input vector reduced to a size of one in the
sequence dimension. After creation of the stateful model, the weights
for the networks were copied from the batch model.

All computations and training were done on a Linux workstation (Ubuntu
Mate 16.04) with 32 GB of ram, i5-2405S CPU @ 2.50GHz and an Nvidia
Geforce GTX 1060 graphics card with 6 GB of ram.




§ RESULTS


Loading and vectorization of the pI dataset resulted in vectors with
the shapes and (7364, 50, 23) and (734, 50, 23) for the train and
test set. The shape is Samples, Sequence, Features (amino acids +
padding). An example of a one-hot encoded peptide is shown in Figure
<ref>.



The training loss and validation loss improved and plateaued out over
the course of training as shown in Figure <ref>.

 There is no sign of over training as there is no rise in the validation
loss function towards the end of training. The final losses were 0.029
and 0.078, for the training and test (validation) loss, respectively.
The spatial difference between the final loss on the test and training
sets appear enlarged due to the logarithmic Y-axis, but makes it easier
to visually see that the training have fully converged.

The pIs of the training and test sequences were predicted with the
final model and used in scatter plots and for calculating the root
mean square error (RMSE) and correlation coefficient (R2). (Figure
<ref> and Table <ref>).







  §.§.§ Challenging with constructed sequences


The model was challenged with designed peptide sequences in the hope
to elucidate the effects of amino acid type, position and neighbor
sequence. The normal acidic and basic amino acids was put in the middle
of an alanine peptide sequence and the pI predicted. The predicted
pI is compared to the side chains pKa value in Table <ref>.

 Standard one-letter amino acid codes are used. The predicted pI should
be a balance between the N- and C- terminal groups and the acidic
and basic side chains. The changes in pI observed, largely follows
the ranking observed for the pKa values, with cysteine (C) as an outlier. 

Table <ref> shows the effect of placing the amino acids
in the N and C-terminal respectively, using glutamate (E) and lysine
(K) as examples of an acidic or basic amino acid respectively.

 The expectation is that the proximity to the N terminal amino group
or C terminal carboxylic acid moiety, would affect the observed pKa
and thus charge of the side chain either up and down with subsequent
adjustments of the predicted pI. For the glutamate the pI is lowered
when placed in the N-terminal end and higher when placed in the C-terminal
end (to the right). The proximity between the carboxylic acid moiety
of the and the N-terminal ammonium groups positive charge, should
make deprotonation more likely at lower pH, and thus lower the pI
as also observed. The effect is opposite in the other end as also
expected. For the ammonium group of lysine, the effect is reversed.

Some examples of the effect of clustering charged amino acids on the
predicted pI is shown in Table <ref>.

 The argumentation is as above, that the proximity of the same charged
amino acids will diminish their acid/base strength. Clustering the
aspartate (D) side chains make the predicted pH slightly higher, where
the lowering of the predicted pI for the basic amino acids (K and
R) is more pronounced. Placing opposite charged side chains besides
each other, enhances the effect of the one closest to the predicted
pI as also observed in the two bottom rows of Table <ref>.




  §.§.§ Prediction in stateful mode


The batch trained model was converted into a stateful model, which
do not predict on a sequence, but rather predicts the pI by being
presented one amino acid at a time, with the current prediction being
affected by the previous predictions via the internal states of the
LSTM cells. This allows a prediction profile to be obtained by iterating
over the amino acids in a sequence and following the predicted pI.
An example of such a prediction profile is shown in Figure <ref>.

 An initial adjustment is observed in the padding phase, but as soon
as an amino acid is encountered (here lysine, K), the pI jumps up.
The pI is largely kept until an acidic aspartate (D) is encountered.
The next aspartate leads to a further drop in pI, but not as pronounced.
After this adjustment the pI is largely unaffected by the remainder
of the encountered acidic aspartates or amino acids without pronounced
acid/base properties.

The initial changes observed during the padding phase, suggests that
the prediction could be influences by the length of the padding (first
part of Figure <ref>). The model was used
for prediction of the same peptide sequence with varying lengths of
padding. The result is shown in Figure <ref>.

 It has a clear effect on the prediction, suggesting that the model
has fine-tuned the prediction based on the size of the peptide, although
for this example, it would have been beneficial to remove the padding
sequence. The error is however in the low range on the second decimal
of the pI prediction.




§ DISCUSSION


Benchmarks of pI calculation pKa sets and pI prediction algorithms
find RMSEs in the range 0.21 to 2.0<cit.>.
However, comparison of the datasets published used by Koslowski<cit.>
on the website<cit.>, with the original data<cit.>,
show that he used the predicted pI from the supplementary information
and not the experimentally related average of the gel fractions, which
invalidates the hard work with the benchmark and the tuning of pKa
values. An advanced SVM model obtains RMSE on a test set of 0.21<cit.>,
but the SVM model had originally<cit.> been tuned
on part of the dataset that comprised the benchmark, and the performance
is likely overestimated. The originally RMSE found was 0.32<cit.>.
It thus appears that the approach using RNNs with LSTM cells reach
a performance at least on par with previous methods.

The effect on apparent pKa of glutamate placed at the termini was
studied by Cargile and coworkers, who used a genetic algorithm to
adjust pKa values based on proximity to the terminals and other charged
residues. They found corrections in the range -0.414 for the N-terminus
and 0.014 for the C-terminus. However, this can't be directly compared
to the pI shift of -0.24 and 0.17 observed with the RNN (c.f. Table
<ref>), as one is calculated on the pKa level and the
other directly on the predicted pI. The trend is nevertheless the
same and it appears that the neural network has implicitly incorporated
the same necessary correction due to the position of the glutamate
residue.

A problem with the use of data driven approach could be that the experimental
values in databases are not of relevance. An example could be that
iodoacetamide peptide sequences are reported and stored as if they
had free cysteine side chains<cit.>.
Another common problem, seem to be that the databases gets polluted
by predicted values from existing models, as it is not always clear
in scientific literature if the pI is predicted or measured. This
could likely pose a problem for a data driven approach as employed
here. However, the model obtained from the data driven approach would
then represent a consensus model of several different models and tuned
pKa sets employed, which could lead to an overall better performance
of the new model due to meta-effects. If only one method had been
employed to derive the pIs for the data driven approach, the new model
would surely regress towards the existing model entirely. Some of
the dataset is also derived from proteomics studies<cit.>,
where the data has also been filtered with regard to the position
in the isoelectric focusing gel and its deviance from the theoretically
predicted value. This also introduces a bias towards the existing
models, which could likely improve the results, as this filtering
will affect both the training and test set. If the dataset is polluted
by predicted pI values and/or filtered based on predicted values real
world performance of the algorithm would likely be lower than here
observed.

Being a data driven approach, rather than a model based approach,
the method could likely be sensitive to peptide sequence distributions
and unusual local sequences. The database consists of peptides measured
with isoelectric focusing and results from proteomics, which are probably
mostly of eukaryotic and procariotic origin. If peptides following
a different amino acid distribution or being of very low homology
to the ones in the training set (such as from Archea or constructed
peptides) the prediction error could likely be higher. However, as
the model seem mostly influenced by the amino acids with known acid-base
properties (c.f. Figure <ref>), the rise in
error could possibly be low. This error could maybe be estimated by
division of the whole peptide dataset into a training and test set
based on clusters found by clustering by sequence homology and redoing
the training and testing.

A better curated or un-biased dataset could likely lead to an increase
in model accuracy and maybe precision. Other improvements could be
the inclusion of modifications if enough data could be collected.
For pIs of peptides with phosphorylated threonines/serines and acetylated/amidated
terminals, this data seems available<cit.>.




§ CONCLUSION


The feasibility of data driven modeling of peptide isoelectric point
with RNNs has been demonstrated. The model obtained seems on par with
literature sources, illustrating the promise of the approach. However,
the performance of the model is likely overestimated due to shortcomings
of the underlying dataset. The model is mostly influenced by the existence
of amino acids with known acidic and basic properties with minor corrections
due to position in the chain and nearby amino acids. Running the model
in a stateful mode gives a minor error when removing the padding sequence.
A release of the source code and trained weights on the authors GitHub
profile is planned.




§ CONFLICT OF INTERESTS


E. J. Bjerrum is the owner of Wildcard Pharmaceutical Consulting.
The company is usually contracted by biotechnology/pharmaceutical
companies to provide third party services. 

elsarticle-num


