## 2 Materials and Methods
### CNN Deep Learning  
The structure of the used CNN model is shown in Figure 1. The CNN model consists of multiple blocks/layers that map the embedded features to the target PTM site class (0 for negative or 1 for positive):  
(i) Input: \begin{table}
\begin{tabular}{c c c} \hline \hline
**PLM** & **\# of Features in Ngly Dataset** & **\# of Features in Kgly Dataset** \\ \hline ProtBERT-BFD & 25,600 & 31,744 \\ ProtBERT & 25,600 & 31,744 \\ ProtALBERT & 102,400 & 126,976 \\ ProtXLNet & 25,600 & 31,744 \\ ESM-1b & 32,000 & 39,680 \\ TAPE & 19,200 & 23,808 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: The size of embedding features that are extracted for a sample of both Ngly and Kgly datasets based on each used PLM. The input to the CNN model is the extracted features from PLMs embedding. It is a 3D array with size \((m,n,f)\) where \(m\) is the number of samples, \(n\) is the length of the sample and \(f\) is the number of features for each token in the sample. (ii) Convolutional block: A convolutional block consists of three sublayers: the first is the CNN1D layer, batch normalization, and activation function. Unlike fully connected networks, CNNs includes conditional connections to calculate local features called feature map. The filters and kernel size (\(k\)) that represent the output of CNN1D and the length of the CNN1D window, respectively, are used as parameters by CNN1D to handle convolution operations. Depending on the kernel (\(k\)), multiple feature extractors \(K\_i\) where \(i=\{1,2,\dots,k\}\), represented as follows [55]:  
\[Z_{i}^{(l)}=H_{i}^{(l-1)}*K_{i}^{l} \tag{3}\]  
where * represents the convolution operation, \(H_{i}^{(l-1)}\) represents the grid representation for the previous layer. The convolution outputs \(Z_{1}^{(l)},Z_{2}^{(l)},\dots,Z_{k}^{(l)}\) are treated as slides and compounded to form a large tensor of convolution layer. The previous equation represents the convolution operation that can be presented as the transformation similar to DNN:  
\[z^{(l)}=W^{l}h^{(l-1)}+b^{(l)} \tag{4}\]  
where \(z^{(l)}\) is the original vector representation of \(Z^{(l)}\), which is gained as a result of the convolution of the input data from the previous layer (\(h^{(l-1)}\)) with the tensor \(W\) (containing filters), and adding bias \(b\). The second sublayer of the convolution block is the batch normalization that typically normalizes the input data for the next layers, reduces the risk of overfitting, allows higher learning rates, and makes the training of the CNN faster and more stable [55]. The third sublayer of the convolution block is the activation function which is a non-linear mathematical function used to generate a high-level feature map of CNN. We use the present rectified linear unit (ReLU) activation function which is represented as:  
\[f(x)=max(0,x) \tag{5}\]  
* Flatten layer: The output feature map of the CNN layer is typically flattened, i.e., converted into low-dimension to be fitted to the next fully connected network layers. (iv) Last layer activation function: Finally, a different activation function is applied to the output of the last FNN. This function is a decision function for the final output of the CNN model. Since our data class label is binary, the Sigmoid function is used which gives a value between 0 and 1. The Sigmoid function is represented as: \[\phi(x)=\frac{1}{1-e^{-x}}\] (6)  
CNN model is trained to classify glycosylation/glycation sites into two classes: positive (glycosylated/glycated) or negative (non-glycosylated/non-glycated) using the binary cross-entropy as a loss function and Adam optimizer. This optimizer depends on the learning rate, batch size, and maximum epoch hyperparameters which help the model to adjust its kernels or weights during the training. The early stopping strategy is utilized on the cost of the validation set to terminate the training process when the accuracy does not increase with the patience value \"8\".