PhD7Faster 2.0: predicting clones propagating faster from the Ph.D.-7 phage display library by coupling PseAAC and tripeptide composition

Benchmark datasets

The dataset used to develop the predictor was acquired from (Matochko et al., 2014). Derda et al. employed high-throughput sequencing technology to characterize both the naïve Ph.D.-7 phage display library and the same library after one round of amplification. By comparing the abundance of each peptide before and after amplification using Bioconductor package edgeR, 770 unique peptides were identified with significantly higher growth rate (parasitic sequences) (Matochko et al., 2014), which were collected into the positive training dataset. The negative dataset was composed of those peptides with the copy number of one in the amplified Ph.D.-7 phage display library. The datasets were then processed as follows: (i) peptide sequences containing ambiguous residues (such as “X”, “B,” and “Z”) were excluded; (ii) sequences within 2 Hamming distance (h = 2, the Hamming distance between two strings of equal length is the minimum number of substitutions required to change one string into the other.) were removed. Finally, 749 peptides were retained in the positive dataset. To match the size of the positive dataset, we randomly selected 749 peptides from the negative dataset. No overlapping was found between the negative and positive datasets. Finally, the benchmark dataset was composed of 749 fast-growing peptides and 749 regular-growing peptides (See positive.fasta and negative.fasta in Supplementary Data).

PseAAC and tripeptide composition

Extracting a set of informative features is a standard and important procedure for developing predictors. Chou initially formulated the PseAAC (Chou, 2001, 2005), which consists of more than 20 discrete numbers, where the top 20 represent the classical amino acid composition (AAC) of a protein sequence whereas the additional parameters incorporate some sequence-order information. PseAAC and tripeptide composition have been widely used in protein prediction related research (Chou, 2011; Lin et al., 2013). Here, they were employed to encode each peptide in the benchmark dataset.

Given a peptide P with L amino acid residues: (1) $\text{P}=\left(R_1{R}_2{R}_3{R}_4{R}_5{R}_6{R}_7\dots R_L\right)$where R_i (i = 1, 2, 3 . . . L) is the residue at the ith sequence position. Accordingly, any sequence like the peptide P of Eq. (1) can be presented using a set of feature vectors with 8,000 + nλ dimensions.

(2) $$\text{P}=\left(P_1,P_2,\cdots ,P_{8,000},P_{8,000}{}_{+1},\cdots ,P_{8,000+n\lambda }\right)$$where the first 8,000 numbers P₁, P₂, . . ., P_8,000 reflect the effect of the conventional tripeptide composition; the remaining nλ elements P_8,000+1, P_8,000+2, . . ., P_8,000+nλ reflect the amphipathic sequence-order pattern. These features are calculated through the following equations: (3) $$P_u=\{\begin{array}{l}\frac{f_u}{{\displaystyle {\sum }_{i=1}^{8,000}{f}_i+w}{\displaystyle {\sum }_{j=1}^{n\lambda }{\tau }_j}}(1\le u\le 8,000)\hfill \\ \frac{w{\tau }_u}{{\displaystyle {\sum }_{i=1}^{8,000}{f}_i+w}{\displaystyle {\sum }_{j=1}^{n\lambda }{\tau }_j}}(8,000+1\le u\le 8,000+n\lambda )\hfill \end{array}$$where f_i (i = 1, 2, 3, . . ., 8,000) are the normalized occurrence frequencies of the 8,000 tripeptides in peptide P; w is the weight factor for the sequence-order effect; (τ_j (j = 1, 2, . . ., nλ) is the j-tier sequence-correlated factor as formulated by: (4) $$\{\begin{array}{l}{\tau }_1= \frac{1}{L-1}{\displaystyle \sum _{i=1}^{L-1}}{H}_{i,i+1}^1\\ {\tau }_2= \frac{1}{L-1}{\displaystyle \sum _{i=1}^{L-1}}{H}_{i,i+1}^2\\ \hfill \dots \hfill \\ {\tau }_n=\frac{1}{L-1}{\displaystyle \sum _{i=1}^{L-1}}{H}_{i,i+1}^n\\ {\tau }_{n+1}=\frac{1}{L-2}{\displaystyle \sum _{i=1}^{L-2}}{H}_{i,i+2}^1\\ {\tau }_{n+2}=\frac{1}{L-2}{\displaystyle \sum _{i=1}^{L-2}}{H}_{i,i+2}^2\\ \hfill \dots \hfill \\ {\tau }_{2n}=\frac{1}{L-2}{\displaystyle \sum _{i=1}^{L-2}}{H}_{i,i+2}^n\\ \hfill \dots \hfill \\ {\tau }_{n\lambda -1}=\frac{1}{L-\lambda }{\displaystyle \sum _{i=1}^{L-\lambda }}{H}_{i,i+\lambda }^{n-1}\\ {\tau }_{n\lambda }=\frac{1}{L-\lambda }{\displaystyle \sum _{i=1}^{L-\lambda }}{H}_{i,i+\lambda }^n\end{array}$$where $H_{i,j}^n$ is the physicochemical property correlation function and can be computed according to the following equation: (5) $$H_{i,j}^n=h^n(R_i).h^n(R_j)$$

where hⁿ(R_i) and hⁿ(R_j) are the values of the nth type of physicochemical property of R_i and R_j in Eq. (1), respectively. It is noteworthy that before substituting the values of all physicochemical properties into Eq. (5), they were undergone a standard conversion as described below: (6) $$h^k(R_i)=\frac{h_0^k(R_i)-{\displaystyle {\sum }_{\alpha =1}^{20}{h}_0^k(R_{\alpha })/20}}{\sqrt{{\displaystyle {\sum }_{u=1}^{20}{\left(h_0^k(R_i)-{\displaystyle {\sum }_{\alpha =1}^{20}{h}_0^k(R_{\alpha })/20}\right)}^2}}}$$

where R_i (i = 1, 2, . . ., 20) denotes the 20-standard amino acid in the alphabetical order of their single-letter codes. $h_0^k\left(R_i\right)$ is the initial value of the kth type of physicochemical property for amino acid residue R_i. Nine kinds of physicochemical properties, namely hydrophobicity, hydrophilicity, mass, pK1, pK2, pI, rigidity, flexibility, and irreplaceability, were considered in this report.

Feature selection

Generally, not all features make an equal contribution to the prediction system. A part of features make significant contributions, while some others make less important contributions (Zhao et al., 2016). Feature selection, thus, is a critical step to reduce feature dimensionality and build a highly effective prediction model (Su et al., 2018; Tang, Chen & Lin, 2016). In this work, the fselect.py program in the LIBSVM 3.23 package was applied to evaluate each feature’s significance to the classification system (Chang & Lin, 2011). As a consequence, each feature corresponds to an F-score. The greater F-score implies the larger importance of the corresponding feature to the prediction model. We rearranged all features by F-scores in descending order. The incremental feature selection strategy was then utilized to determine the optimal feature subset (He et al., 2016b; Li et al., 2017), which can produce the maximal accuracy. Feature selection was conducted as follows: (i) investigating the accuracy of the first feature subset which included the feature with the largest F-score; (ii) examining the accuracy of the second feature subset that was generated by appending the feature with the second largest F-score; (iii) iterating the second step from the larger F-score to the smaller F-score until all candidate features were added. The best feature subset with the highest accuracy can be finally obtained.

Support vector machine

The SVM is a powerful supervised learning method, which has been widely applied in classification (He et al., 2016b; Kang et al., 2018; Li et al., 2017; Ru et al., 2014) and regression analysis. In this study, we utilized the LIBSVM 3.23 program (Chang & Lin, 2011) that could be freely available for download from http://www.csie.ntu.edu.tw/∼cjlin/libsvm/. We chose the radial basis function kernel as the kernel function. The optimal kernel width parameter γ and penalty constant C were selected by using the parameter selection tool in the LIBSVM 3.23 (Chang & Lin, 2011).

Performance evaluation

The 10-fold cross-validation was adopted to evaluate the predictive model in this study. Four commonly-used parameters, including sensitivity (Sn), specificity (Sp), accuracy (Acc), and Matthews correlation coefficient (MCC), were employed to investigate the performance of the constructed model. These measures were expressed as follows: (7) $$\text{Sn}=\frac{\text{TP}}{\text{TP}+\text{FN}}$$ (8) $\text{Sp}=\frac{\text{TN}}{\text{FP}+\text{TN}}$ (9) $$\text{Acc}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{FN}+\text{FP}+\text{TN}}$$ (10) $$\text{MCC}=\frac{\text{TP}\times \text{TN}-\text{FP}\times \text{FN}}{\sqrt{\left(\text{TP}+\text{FP}\right)\left(\text{TP}+\text{FN}\right)\left(\text{TN}+\text{FP}\right)\left(\text{TN}+\text{FN}\right)}}$$where TP and TN denote the number of true positives and negatives, respectively. FP and FN are the number of false positives and negatives, respectively. The area under the receiver operating characteristic (ROC) curve (AUC) was also calculated as a performance measure. The AUC ranges from zero to one. The AUC of one represents a perfect prediction, 0.5 a random guess.

To estimate the statistical significance of the predictive accuracy, a permutation test with 1,000 shuffles was performed by exchanging the labels of the benchmark dataset. The 10-fold cross-validation was then conducted against the label-rearranged dataset. Thus, each permutation trial corresponds to an accuracy value. The p-value was calculated by the number of permutations that the Acc produced by the permuted dataset was higher than Acc based on the un-permuted dataset divided by the overall shuffle times. p-values of <0.05 were referred to as statistically significant.

Standalone version implementation

The standalone version of PhD7Faster 2.0 was developed with open source Qt 5.7 under the GPL & LGPLv3 licenses, which uses standard C++ for developing multiple-platform applications. Both graphical user interface (GUI) and command-line versions of PhD7Faster 2.0 were implemented. We provided different versions for Windows and Linux systems with little or no modification. All versions and source code are freely available at http://i.uestc.edu.cn/sarotup3/download.html.

Parameter optimization

Two important parameters: λ and w in Eq. (3) were necessary to be optimized before building the model. To obtain the best parameters, multiple experiments were performed according to the following standard: (11) $$\{\begin{array}{l}1\le \lambda \le 6\text{ with step }\Delta =1\hfill \\ 0.05\le w\le 0.70\text{ with step }\Delta =0.05\hfill \end{array}$$

Thus, a total of 6 × 14 = 84 individual combinations were obtained. Then, we used the 10-fold cross-validation to investigate the accuracy of the model, which was built with SVM and the feature set without feature selection. λ = 3 and w = 0.15 produced the highest accuracy, which was considered as the best parameter combination.

Performance of PhD7Faster 2.0

The optimal feature subset with 644 features was determined through feature selection against 8,027 features including 8,000 tripeptide features and 27 PseAAC features. The SVM-based model was then trained with the optimal feature set. The results from the 10-fold cross-validation showed that the Acc of the predictive model was 81.84% with MCC of 0.64, Sn of 84.51%, and Sp of 79.17% when the threshold to distinguish between predicted positives and negatives (tp) was set to be 0.5. The ROC curve for model tuning is shown in Fig. 1, where the AUC is approximately 0.90. The permutation test resulted in a p-value of < 0.001. The above results indicated that PhD7Faster 2.0 achieved a promising performance.

Web and standalone versions of PhD7Faster 2.0

For the convenience of users, the SVM-based predictive model was implemented into a user-friendly web server, called PhD7Faster 2.0, which is freely available at http://i.uestc.edu.cn/sarotup3/cgi-bin/PhD7Faster.pl. The standalone GUI and command-line versions of PhD7Faster 2.0 for Windows and Linux systems were also provided. The interface, as well as the utilization of the GUI version, is remarkably similar to those of the web version (Fig. 2). A dataset with 20,000 peptides from the Ph.D.-7 phage display library was constructed (see testdataset.fasta in Supplementary Data). The standalone PhD7Faster 2.0 can complete analysis of the dataset within 60 s on a regular computer with Intel Core i3 Processor and 4GB RAM, which suggests that PhD7Faster 2.0 is highly efficient in processing massive datasets. PhD7Faster 2.0 was integrated into the SAROTUP 3.0 suite, which contains a series of computational tools to identify TUPs.

Comparison between PhDFaster 2.0 and 1.0

Parasitic sequences were identified significantly enriched in the amplified Ph.D.-7 phage display library by differential enrichment analysis of naïve and amplified Ph.D.-7 phage display libraries (Matochko et al., 2014). These parasitic peptides were grouped into the positive dataset of PhD7Faster 2.0. However, the positive dataset of PhD7Faster 1.0 was constructed based on threshold in copy numbers after one round of amplification in one replicate of sequencing data, irrespective of copy numbers in the naïve Ph.D.-7 library. Peptides with high abundances in both the naïve and amplified Ph.D.-7 libraries may also be selected as fast-growing sequences. Therefore, the positive training dataset of PhD7Faster 2.0 is more reliable than that of PhD7Faster 1.0.

PhD7Faster 2.0 was developed based on the combination of PseAAC and tripeptide composition, whereas only dipeptide composition was employed to build PhD7Faster 1.0. We also tried to use dipeptide composition to encode each peptide in the training dataset of PhD7Faster 2.0, but only 64% accuracy was obtained in the 10-fold cross-validation after feature selection. PseAAC coupled with tripeptide composition has been used in multiple protein prediction fields, such as predicting the subcellular localization of mycobacterial proteins (Zhu et al., 2015) and predicting apoptosis protein subcellular location (Liao et al., 2011). They contain more sequence-order information than dipeptide composition and hence can better reflect the feature of a peptide sequence. Thus, PhD7Faster 2.0 has 5% sensitivity, 2% accuracy and 0.04 MCC higher than PhD7Faster 1.0 (Table 1).

The standalone PhD7Faster 2.0 is empowered to identify PrTUPs within output of conventional phage display as well as large next-generation sequencing data, whereas PhD7Faster 1.0 can only work with small-scale data sets (several hundreds of peptides). This important improvement makes PhD7Faster 2.0 as an enhanced and powerful tool for scanning and reporting PrTUPs from the Ph.D.-7 phage display library. The emergence of PhD7Faster 2.0 highlights the significance of high throughput sequencing of different types of phage display libraries and developing bioinformatics tools for identifying PrTUPs from these libraries.

PhD7Faster 2.0 cannot predict the censorship in the Ph.D. libraries

It is possible that some peptides are likely to be censored from being displayed on the phage in the first place. The censorship of positively charged amino acids has been reported since these residues suppress proper insertion of pIII into the inner membrane of Escherichia coli, thus decreasing the efficiency of the assembly and extrusion of phage clones (Peters et al., 1994). Rodi, Soares & Makowski (2002) also observed that peptides of α-helix or β-sheet conformations were censored in Ph.D.-12 and Ph.D.-C7C libraries. Steiner et al. (2006) have shown that maturely folded proteins are displayed poorly via the Sec translocation pathway. However, this censorship is a completely different phenomenon from that of phage growing faster. Therefore, PhD7Faster 2.0 is not able to predict this censorship.

PhD7Faster 2.0 predict PrTUPs in the Ph.D.-7 library

The PrTUPs have significantly higher proliferation rates than normal-growing phage and are favored during the amplification steps. The proliferation advantage of some PrTUPs have been verified to be intrinsic to mutations in the 5′-untranslated region (UTR) of gene II in M13 phage (Brammer et al., 2008; Nguyen et al., 2014; Zygiel et al., 2017). Zygiel et al. (2017) have also described the likelihood that these mutations compensate for the replication defect afforded by the lacZα insert present in the M13 bacteriophage-based vector upon which the Ph.D.-7 (and Ph.D.-12) library was based. Thus, the particular peptide displayed (e.g., HAIYPRH, GKPMPPM, AKIDART) is merely a stowaway on a clone that propagates fast due to its gene II 5′-UTR mutation(s). In these clones, the peptide itself is completely arbitrary; it just happens to be the peptide displayed on a clone that picked up a mutation prior to or during library construction. As these mutations in the phage genome are unrelated to the displayed peptide, PhD7Faster 2.0 may not be able to predict this type of PrTUPs in the Ph.D.-7 library. In addition, Smith et al. indicated that the enhanced propagation rate of some PrTUPs may be due to the displayed peptide (Thomas, Golomb & Smith, 2010), and PhD7Faster 2.0 can be used to predict this type of PrTUPs in the Ph.D.-7 library. However, no direct evidence supports that displayed peptides allow the phage to propagate faster, and the biological mechanism remains to be further examined.