An Alignment-Free Algorithm in Comparing the Similarity of Protein Sequences Based on Pseudo-Markov Transition Probabilities among Amino Acids

In this paper, we have proposed a novel alignment-free method for comparing the similarity of protein sequences. We first encode a protein sequence into a 440 dimensional feature vector consisting of a 400 dimensional Pseudo-Markov transition probability vector among the 20 amino acids, a 20 dimensional content ratio vector, and a 20 dimensional position ratio vector of the amino acids in the sequence. By evaluating the Euclidean distances among the representing vectors, we compare the similarity of protein sequences. We then apply this method into the ND5 dataset consisting of the ND5 protein sequences of 9 species, and the F10 and G11 datasets representing two of the xylanases containing glycoside hydrolase families, i.e., families 10 and 11. As a result, our method achieves a correlation coefficient of 0.962 with the canonical protein sequence aligner ClustalW in the ND5 dataset, much higher than those of other 5 popular alignment-free methods. In addition, we successfully separate the xylanases sequences in the F10 family and the G11 family and illustrate that the F10 family is more heat stable than the G11 family, consistent with a few previous studies. Moreover, we prove mathematically an identity equation involving the Pseudo-Markov transition probability vector and the amino acids content ratio vector.


Introduction
With the recent development of next-generation sequencing technologies, there has been an explosion in the numbers of available DNA and protein sequences. The numerous newly sequenced protein sequences present an urgent need for novel computational algorithms to compare their similarities with sequences from known protein families, to predict their structures, and thus to infer their functions [1][2][3][4][5][6].
As usually the first step in a bioinformatics pipeline, sequence comparison is very crucial since it affects all down-stream analyses. Popular methods for sequence comparison generally fall into two categories: those using sequence alignment and those using alignment-free methods. In a sequence alignment, a score function is used to represent insertion, deletion, and substitution of nucleotides or amino acids in the compared DNAs or proteins, and the objective is to identity the alignment with the highest overall alignment score through methods like dynamic programming and seeding [7][8][9]. However, sometimes alignment becomes misleading due to unequal lengths of sequences, gene rearrangements, inversion, transposition, and translocation at substring level [10]. In these scenarios, alignment-free methods present good alternatives to alignment methods, which usually quantify sequence similarities using K-mer frequencies and other sequence features [11].
An alignment-free method for comparing protein sequences usually consists of two steps. At first, the protein sequences are transformed into fixed-length feature vectors [12]. The feature vectors are then fed into a vector similarity comparison algorithm to perform downstream analysis like phylogenetic inference [13]. Feature extraction is a procedure to extract desired information from the query sequences, which is usually critical to the accuracy of an alignment-free method [14]. Widely accepted features include chemical and physical properties [15], distance frequency matrix [16], K-string dictionary [17], 2D and 3D amino acid adjacency matrices [18], pseudo amino acid composition [19], and sequential and structural evolution information [20,21]. Though these methods have their own advantages, they are suffering problems like computational intensive and low accuracy. Thus, more discriminatory features are still in demanding.
To further improve protein sequence comparison accuracy, we present a novel 440 dimensional feature vector, which models a few important information of a protein including the amino acids' abundance and position information, and the Pseudo-Markov transition probabilities among them. We then test the performance of our feature vector in two well studied datasets: (1) the ND5 dataset [22] and (2) the F10 and G11 dataset [23]. They have been widely used in evaluating protein comparison algorithms [22,24]. As a result, our method is more accurate than a few existing methods for similarity analysis on the ND5 dataset, and we achieve accurate phylogenetic tree and heat stability results on the F10 and G11 dataset.

Method
Amino acid composition and distribution are two most fundamental information about a protein sequence. They have been widely used and proven to be effective in protein sequence analyses [25], structural classification [26][27][28], pattern recognition receptor prediction [29], and fold recognition [30]. Thus, we proposed a novel representation for a protein sequence based on the two features, i.e. a 440-D feature vector consisting of (1) a 400-D Pseudo-Markov transition probability vector reflecting the order information of adjacent amino acids. (2) a 20-D amino acid content ratio vector describing the frequency of each amino acid in the sequence, and (3) a 20-D amino acid position ratio vector exhibiting the position distribution of each amino acid.

Construction of the 400 dimensional Pseudo-Markov transition probability vector
Let S = S 1 S 2 Á Á ÁS N be a protein sequence of length N defined on A = {A 1 , A 2 , Á Á Á, A 20 }, an ordered alphabet of 20 amino acids. For 1 i, j 20, 1 k N and 1 l N − 1, an amino acid A i is said to occur at position k if S k = A i , and an ordered amino acid pair A i A j is said to occur at position l if S l S l+1 = A i A j . Let n i be the number of occurrences of A i and n i,j be the number of occurrences of A i A j in S. We then define the 400 dimensional vector as (P 1,1 , P 1,2 , Á Á Á, Foundation for Advanced Talents (No A201400121 to YZ), the Educational Commission of Hebei Province on of Humanities and Social Sciences(No SZ16180 to YZ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. P 1,20 , P 2,1 , P 2,2 , Á Á Á, P 2,20 , Á Á Á, P 20,1 , P 20,2 , Á Á Á, P 20,20 ), where In particular, if there is no A i or A i appears only once at the end of S, then the numerator and denominator of P i,j are both 0. In this case, we define P i,j = 0.

Construction of the 20 dimensional amino acid content ratio vector
Given that the protein sequence is composed of only 20 amino acids, it is clear that , we define its content ratio C i as C i ¼ n i N and the 20 dimensional amino acid content ratio vector as (C 1 , C 2 , . . ., C 20 ). Obviously, Property. Suppose S 1 = A u and S N = A v for indices u and v with 1 u, v 20. Then for any 1 j 20, we have

Construction of the 20 dimensional amino acid position ratio vector
In particular, if A v occurs only once in S, i.e. n v = 1 then Proof. If j = u, from eqs (1) and (3) we have Finally, let n v = 1. By definition, we have n v,j = 0 and P v,j = 0 for any 1 j 20. Thus, completing the proof.
Quantifying the distances among protein sequences based on their feature vectors Let S and T be two proteins and V S and V T be their 440-D feature vectors. Then the distance between S and T is quantified by the Euclidean distance between V S and V T , that is,

Results and Discussions
To evaluate the performance of our method, we applied it into two datasets: (1) the ND5 dataset [22] and (2) the F10 and G11 dataset [23].

Datasets
The ND5 dataset consists of the ND5 protein sequences of 9 species including human, gorilla, pigmy chimpanzee, common chimpanzee, fin whale, blue whale, rat, mouse, and opossum ( Table 1). The sequences have lengths 602~610 base pairs (bps). It is a popular benchmark data for testing the performances of computational methods in comparing the similarity of protein sequences [15,[31][32][33][34]. The F10 and G11 datasets represent two of the xylanases containing glycoside hydrolase families, i.e., families 10 and 11 respectively. Specifically, the F10 dataset contains ten xylanases with NCBI accession IDs O59859, P56588, P33559, Q00177, P07986, P07528, P40943, P23556, P45703, and Q60041 respectively. The G11 dataset also consists of ten xylanases with NCBI IDs P33557, P55328, P55331, P45705, P26220, P55334, Q06562, P55332, P55333, and P17137 respectively. Application to the ND5 dataset We first encoded the nine protein sequences into 440-D feature vectors. In Figs 1 and 2, we showed the content ratios and position ratios of the twenty amino acids over the sequences. As can be seen, the content ratios and position ratios exhibit similar trends over the twenty amino acids. The amino acid L has the highest content ratio and position ratio over all 9 sequences whereas amino acid C has the lowest content ratio and position ratio. In addition, the 9 species are quite similar according to the amino acids distributions of both the content ratio and position ratio in the ND5 protein.
We then calculated the pairwise Euclidean distances among the nine 440-D feature vectors and showed the results in Table 2. As we can see, human, P.chim, C.chim, and gorilla are closer to each other and they are relatively far from rat, mouse and opossum. For a better view, we also plotted a heat-map based on the distances in Fig 3. In order to estimate the contribution of each part in the 440-D feature vector to the final performance in sequence similarity analysis, we plotted heat maps for the ND5 dataset based on the 20-D amino acid position ratio vector (see Fig 4), the 20-D amino acid content ratio vector (see Fig 5), and the 40-D amino acid position ratio and content ratio vector (see Fig 6), respectively. Clearly, Fig 3 is  A common strategy to evaluate an alignment-free method is to compare it with a popular alignment method like ClustalW [31], which has a much higher time and space complexity than alignment-free methods. Table 3 showed the pair-wise distances of the 9 protein sequences using ClustalW (i.e. Table 4 in [31]). We calculated the correlation coefficient between the distances from our method and those from ClustalW and compared our method with a few popular alignment-free methods [15,[31][32][33][34] using this coefficient as a criterion (see Table 4). Comparing Protein Sequences Based on Pseudo-Markov Transition Probabilities among Amino Acids As Table 4 shows, the correlation coefficient between our method and ClustalW is 0.962, which is the highest among the 6 methods. As a result, our method is more consistent with ClustalW than the other 5 methods, which indicates that our method is more accurate.

Application to the F10 and G11 dataset
We also tested our method on the F10 and G11 datasets and plotted the heat map based on the pair-wise Euclidean distances in Fig 7. As can be seen, our method accurately separated the sequences in family F10 with those in G11 with the F10 xylanases locating in the top right quarter and G11 xylanases in the lower left quarter. We also observed that the F10 dataset is more heat stable than the G11 dataset, which is consistent with other studies, e.g., [15].
It is of note that we applied the Euclidean distance in quantifying the distances among the feature vectors for different proteins. Euclidean distance is one of the simplest and most intuitive distance measures, which has been adopted in many fields, such as gene identification [35], protein 3D structure reconstruction [36], robust automatic pectoral muscle segmentation [37] and classification of normal and epileptic seizure EEG signals [38], etc. However, there are many other distance measures, which could affect protein similarity analysis. As an example, we compared the Euclidean distance with the Hamming distance for the ND5 dataset and      For the ND5 dataset, we further computed the agreement (i.e., the Pearson correlation coefficients between the protein similarity matrices) between our method (based on the Hamming distance) and ClustalW, which is 0.937, a little bit less than that for the Euclidean distance (0.962). Thus, we believe that the Euclidean distance is more effective than Hamming distance for these two datasets.

Conclusion
In this paper, we have proposed a novel alignment-free method to compare protein sequences. The method is more accurate than 5 other popular alignment-free methods in the ND5 dataset and is capable of distinguishing the F10 xylanases family from the G11 family. The comparison results of this method are quite consistent with protein sequence aligners like ClustalW. It presents an alternative of these aligners when time and space complexities become an issue.
In the future, a few machine learning methods [39] could be applied to further improve the performance of our method. For example, in contrast to phylogenetic analysis, methods like K-means analysis [40] and random forest [41] could also be applied to classify the proteins and perform taxonomy. However, it is out of the scope of this study. In addition, our novel features could also be applied into applications like essential gene identification [42] and similar problems related to DNAs or RNAs.