key: cord-0002886-jbmj21c5 authors: Du, Qi-Shi; Wang, Shu-Qing; Xie, Neng-Zhong; Wang, Qing-Yan; Huang, Ri-Bo; Chou, Kuo-Chen title: 2L-PCA: a two-level principal component analyzer for quantitative drug design and its applications date: 2017-08-01 journal: Oncotarget DOI: 10.18632/oncotarget.19757 sha: f580401828d99aa0874fc10a731aa34a8cfb69b3 doc_id: 2886 cord_uid: jbmj21c5 A two-level principal component predictor (2L-PCA) was proposed based on the principal component analysis (PCA) approach. It can be used to quantitatively analyze various compounds and peptides about their functions or potentials to become useful drugs. One level is for dealing with the physicochemical properties of drug molecules, while the other level is for dealing with their structural fragments. The predictor has the self-learning and feedback features to automatically improve its accuracy. It is anticipated that 2L-PCA will become a very useful tool for timely providing various useful clues during the process of drug development. With the fast developments of computer-aided drug design (CADD) [1] [2] [3] [4] [5] , currently a number of drug design approaches are developed, and several computer software packs [6, 7] are available that can speed up the discovery of new chemical and biological drugs in more efficient and economical procedure. However, so far we still have no perfect theories, ideal technologies, and faultless software tools that can guarantee complete success of the designed drugs due to the complicity of the interactions between medicinal drugs and their biological targets [8, 9] . The factors and parameters that may affect the bioactivities of drugs are not only from the structure of drug itself, but also from its biological target, coenzymes, and interaction environment [10, 11] . Principal component analysis (PCA) [12] [13] [14] is a useful tool that has been widely used in chemistry, biology, environment, and many fields of social science. The PCA approach has also been used in drug design for many years. Traditionally PCA is a single level and one direction prediction and analysis technique, described as the following equation where { } , x i k � are the physicochemical parameters of the i-th molecule, � � { } a k are the coefficients of molecular parameters, and w i is bioactivity of the i-th molecule [15, 16] . The bioactivity w i could be logarithm of IC 50,i (pIC 50,i =-logIC 50,i ), or binding free energy ΔG°i between drug and receptor. After the coefficients � � { } a k of parameters are solved from the linear equations Eq.1 in a training set of drug candidates, the parameter coefficients � � { } a k can be used to predict the bioactivities of the designed or newly synthesized drug compounds, where K is the total number of molecular parameters. Currently hundreds even thousands of molecular parameters are available for drug design [17, 18] . However for certain drug-receptor interaction system, these parameters are not equally important; actually too many parameters may cause the over correlation problem [19, 20] . In PCA technique only the principle components are selected to describe the bioactivities of drug molecules, and to predict the bioactivities of drug candidates. In the present study, an improved principal component analysis method, the so-called two-level principal component analysis (2L-PCA), is proposed to deal with the extreme complexity and huge amount of parameters in drug design and discovery. In the 2L-PCA predictor, the 1 st level is to deal with the physicochemical properties of drug molecules, and the 2 nd level is to deal with the fragments of molecular structures. The proposed two-level model can not only significantly enhance the prediction power, but also yield more useful information for in-depth analysis. According to Chou's 5-step rule [21] that has been widely used by many investigators (see, e.g., [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] ), to develop a really useful statistical predictor, one should consider the following five procedures: (1) benchmark dataset; (2) sample representation; (3) operation algorithm; (4) cross validation; (5) web-server. Below, let us describe how to deal with them one-by-one. However, to comply with the Journal's rubric style, they are not exactly following the aforementioned order. As an example to show the advantage of 2L-PCA, we applied it for predicting the binding affinity of epitopepeptides with class I MHC molecules HLA-A*0201 [38, 39] . HLA-A*0201 is one of the most frequent class I alleles found in many different species and populations, which plays a critical role for antigen presentation in both viral antigens [40] and tumor antigens from a variety of cancers [41] [42] [43] [44] , and is expressed in approximately 50% of Caucasians population [45] . The epitope-peptides consist of nine amino acids [38, 39] . In the 2L-PCA study for the epitope-peptides, the nine side chains of the nine amino acids are the nine fragments. Eight physicochemical properties are used as the descriptors of the 20 natural amino acids. Four of them are the HMLP parameters [15, 16] , describing the lipophilic character, hydrophilic character, surface area with lipophilic potential, and surface area with hydrophilic potential, respectively. The fifth property is the volume of amino acid side chains. The remaining three properties are the secondary structural potency indices of amino acids: the α-potency, β-potency, and coil-potency [46] . Listed in Table 1 are the eight physicochemical parameters of 20 amino acids used in this study. In this study the HMLP parameters were used to describe the lipophilicity and hydrophilicity of molecular fragments. In peptides the HMLP parameters of the 20 natural amino acid side chains are available from literatures. However, the HMLP parameters of common chemical molecular fragments have to be derived using complicated calculations. In such cases other hydrophobic parameters can be used, e.g., the atom-based hydrophobic parameters in [47] . To reduce computational time, the cross validation in this study was performed via the independent dataset test [48] , as described as follows. The sequences and experimental binding affinities of the 90 peptides were used as the training dataset to train the model, while those of the 40 peptides taken from [49] as the independent dataset to test the model. Actually, such 40 peptides had also been compiled in a series of publications [41, 42, [50] [51] [52] [53] [54] [55] . The logarithms (pIC 50 ) of IC 50 were used as the bioactivity, because they are related to the changes in the free binding energy [55, 56] . Listed in Table 2 are the sequences and the experimental pIC 50 of the peptides used in the training set. The binding strength of the 90 training peptides and 40 testing peptides covers the low, intermediate, and high affinity. The following two criteria were applied in the choice of the testing peptides: (1) the range of binding affinities in the testing dataset should not exceed the range of affinities in the training set; (2) the amino acid at each position in the testing dataset should also be present at that position in the training set of peptides. These two conditions make the 130 peptides to be the ideal benchmark dataset for 2L-PCA method. The iterative 2L-PCA technique described in Method section is used for the binding affinity study of peptides based on the sequences and experimental data listed in Table 2 of fragment parameters were assigned to 1, implying that all fragment parameters are equally important. Shown in Figure 1 are the curves of correlation coefficients R vs iterations , where the curve R a is for the iteration of coefficients � { } a k , and the curve R b is for the iterations of Table 3 . In the iterative solution precedure the correlation coefficien increases from the first value R A (1) =0.4167 to the converged value R A (98) =0.8871, and the prediction residue decreases from the first value Q A (1) =0.7223 to the converged value Q A (98) =0.0387. The predicted pIC 50 of the 40 queried peptides in the testing set are given in Table 4 , which were predicted It is expected that, with more experimental data available, the predictive power of 2L-PCA will be further improved. Actually, 30 prediction servers for human MHC-I peptide molecules were evaluated in a review article [57] . Among the 30 existing servers, 16 were ranked as the first class that provided the most accurate prediction results for MHC-I peptide molecules with the correlation coefficients ranging from r = 0.55 to r = 0.87. It has been shown in this study that the prediction correlation coefficient yielded by our 2L-PCA method is r = 0.868, being ranked around the very top of the first class. 2L-PCA neither needs knowing the exact comformations of the peptides nor needs aligning the peptides according to a template. The two steps are necessary but quite difficult for CoMFA [58, 59] and CoMSIA [60, 61] owing to that there are numerous possible conformations for peptides and that the experimental crystal structure for serving as a template is often not available. 2L-PCA method provides an alternate way for design of the chemical drugs and peptide drugs. The eigenvalues and contributions of physicochemical properties and amino acid positions in peptides are summarized in Table 5 and shown in Figure 3 . In Table 5 the eigenvalues are normalized. The eigenvalue portion of the first three property eigenvectors is almost 100%, and the eigenvalue portion of the first eigenvector alone is larger than 99%. Most contributions are made by the three properties: side chain volume (Vol), lipophilic surface area (S L ), and hydrophilic surface area (S H ), as The molecular structure is divided into 4 fragments according to the substitutes being investigated. The fragments F 1 , F 2 and F 3 are three substituent groups, and the fragment F 4 is the remaining part of the molecular parent. (B) In short peptides each side chain of amino acid residue is a fragment. shown in Figure 3b and Table 5 . The contributions of other 5 properties seem very small. The eigenvalue of the first peptide position eigenvector is larger than 98%. In Table 5 the contributions of the nine amino acid positions are different. However the differences are not big, implying that all positions are almost equally important. The detailed computation results are given in Supplementary Information 1. We are often facing two kinds of challenges in theoretical prediction for drug design: one is over-correlation problem, and the other is lack of information and explanation for the predicted results. The over-correlation problem is caused by large amount of parameters used in the prediction model, which may yield quite good correlation results in self-consistency test [62, 63] , but very poor predicted results in independent dataset test owing to the high dimensional disaster [19] or "curse of dimensionality" problem. To solve this problem, the pseudo amino acid composition (PseAAC) was introduced [64] . Ever since then, the concept of PseAAC or the general PseAAC [21] has been widely used in drug development and biomedicine [65, 66] and nearly all the areas of computational proteomics (see, e.g., [67] as well as a long list of references cited in [68, 69] ). Actually, the physicochemical properties used here can be regarded as some optimal pseudo components [70] . It is through such a PseAAC approach to remove the trivial parameters (or reduce the feature vector's dimension) and grasp the key ones. Besides, the traditional prediction methods fail to provide a good explanation for the predicted results; i.e., how do the physicochemical properties and the structural changes affect the bioactivities? In contrast to that, the proposed "2L-PCA" method can provide more information about the impact of the physicochemical properties and molecular fragments to the bioactivities of drug candidates. In practical drug design and development, usually the basic structure of drug candidates keep constant, only small modifications are made on several fragments. The structure parameters of the entire molecules cannot clearly describe the detailed characters of the small changes at individual fragments or substitutes. In the 2L-PCA model the molecular structures are separated into several fragments, and are described by a set of fragment parameters. An example of molecular structure and its fragments is shown in Figure 4A . The idea of molecular fragments also can be applied to the peptide drugs, in which each side chain of an amino acid is a fragment, as shown in Figure 4B . In the 2L-PCA prediction model the bioactivity w i of molecule i is the summation of contributions Δg i,l from all molecular fragments; i.e., where Δg i,l is the contribution of fragment l to the bioactivity w i of molecule i, b l is the prediction coefficient of fragment l, and L is the total number of molecular fragments. The contribution Δg i,l of fragment l is the summation of the contributions from all physicochemical properties of fragment l, namely where x i,l,k is the physicochemical property k of fragment l in molecule i, a k is the prediction coefficient of physicochemical property k, and K is the total number of physicochemical properties. Inserting the Eq.4 into Eq.3 we get the general equation of 2L-PCA prediction model as given by where N is the total number of molecular samples. Eq.5 can be expressed in vector and matrix form as given below where X N,L,K is the three dimensional (3D) data matrix of molecular parameters, W N is the bioactivity column vector of molecular samples, B L is the coefficient vector of fragments, and A K is the coefficient vector of physicochemical properties. The general three-dimensional 2L-PCA equation of Eq.6 can be reduced to two 2D equations with the following algebra operations, where H N,L is the 2D data matrix of molecular fragments. Substituting H N,L into Eq.6, we obtain the following fragment 2D equation Likewise, the property 2D equation can also be expressed as and F A W where F N,K is the 2D data matrix of physicochemical properties. The fragment 2D equation Eq.8 and the property 2D equation Eq.10 can be solved using the standard algebra method. Both sides of the fragment 2D equation of Eq.8 are multiplied with the transposed matrix H t N,L from left, it follows that and H W S Thus, we get the following symmetrically square matrix equation of fragments Since the fragment square matrix equation of Eq.13 is multiplied by its inverse matrix U -1 L,L , the prediction coefficients B L for the fragments are obtained, as given below where the inverse matrix U -1 L,L can be obtained by solving the eigen equation [71] where Ψ L,L is the eigenvectors and α L is the eigenvalues of fragment square matrix U L,L [72, 73] . Similarly, left-multiplying both sides of property 2D equation of Eq.10 with F t N,K , we have and From Eqs.17-18, we get the following square matrix equation of properties Multiplying Equation Eq.19 with the inverse matrix V -1 K,K , will give the solution of property prediction coefficients A K ; i.e. Thus, the inverse matrix V -1 K,K is obtained by solving the eigen equation of property square matrix V K,K : and where Φ K,K is the eigen-vectors and β k is the eigen-values of the property square matrix V K,K . In the training dataset for drug candidates the two prediction coefficients set A K and B L in the 2L-PCA general equation Eq.6 are solved in an iterative procedure [74, 75] . Firstly the initial fragment coefficients B (0) L are assigned to 1 {b i =1, i=1,2…,L}, implying all fragments are equally important. The initial B (0) L are used in the property 2D equations Eq (9) and (10) L is used to find the A (2) K . Above iterative procedure is repeated for n times, until to reaching a threshold value ε; i.e., Illustrated in Figure 5 is the iterative solution procedure for the 2L-PCA predictor. The property eigenvectors � � { } ϕ k are orthogonal and normalized; i.e., and www.impactjournals.com/oncotarget ϕ ϕ k k j k where the term φ 2 j,k is the component of the j-th property in the k-th eigen-vector φ k . The first K′ property eigenvectors are the principal components whose eigen-values are larger than a threshold (e.g., ε=90% or 95%); i.e., The total contribution γ j of the j-th property to the bioactivity of molecular samples in training set is defined as the following summation, The property eigen-vectors � � { } ϕ k span an orthogonal multiple space, in which a drug molecule P i is a vector, and its projection J i,k on the k-th property-eigenvector φ k is calculated by where f i is the i-th row vector of the property matrix F NK of Eq.9. In the projection J i,k of molecular sample P i on the k-th property-eigenvector j k the component of the r-th property is α k φ 2 r,k , therefore the total contribution of r-th property to the sample P i is the summation of components from all principal property eigenvectors, namely Similarly, the fragment eigenvectors ψ l span an L-dimensional orthogonal space. The first L′ fragment eigenvectors are the principal components. The total contribution factor λ j of the j-th fragment to the bioactivity of peptide set is given by In the same way the projection I i,l of sample P i on the l-th fragment-eigenvector ψ l can be calculated by where h i is the i-th row vector of the fragment matrix H NL of Eq.7. In the projection I i,l of molecule P i on the l-th fragment-eigenvector φ l the component of the r-th fragment is α l ψ 2 r,l , therefore the total contribution of r-th fragment to the sample P i is the summation of components from all principal fragment eigenvectors; i.e., As pointed out in [76] , user-friendly and publicly accessible web-servers represent the future direction for developing practically more useful predictors or any computational tools. Actually, user-friendly web-servers as given in a series of recent publications [23-25, 30, 32, 34-36, 69, 70, 77-91] will significantly enhance the impacts by attracting the broad experimental scientists [66, 92] . We will do our best to establish a web-server for 2L-PCA as soon as possible. Once it has been done, an announcement will be made thorough a publication or our webpage. The 2L-PCA predictor proposed in this paper is a very useful tool for drug design. Its advantages can be summarized as follows. (1) With 2L-PCA, the molecular structures of drug candidates can be separated into several fragments described by physicochemical parameters of the molecular fragments, thus the small modifications on individual fragments can be clearly shown. The authors are very much indebted to the eight anonymous reviewers, whose constructive comments are very helpful for strengthening the presentation of this paper. This work was supported by grants from the National Science Foundation of China (NSFC http://www.nsfc.gov.cn/) under the contract numbers 31360207, 31370716, and 31400079, Guangxi Science and Technology Development Projects (no. 14123001- We thank the National Supper Computing Center (NSCC www.nscc-tj.gov.cn) in Tianjin for the valuable help in the modeling calculations using TH-1A super computer. REFERENCES 1 The many roles of computation in drug discovery Computer-based de novo design of drug-like molecules The role of structure-based ligand design and molecular modelling in drug discovery Structural bioinformatics and its impact to biomedical science Assessing scoring functions for protein-ligand interactions QSAR models for predicting octanol/water and organic carbon/water partition coefficients of polychlorinated biphenyls Docking and scoring in virtual screening for drug discovery: methods and applications Application of the three-dimensional structures of protein target molecules in structure-based drug design Computational methods for biomolecular docking Automated docking with grid-based energy evaluation Amino acid principal component analysis (AAPCA) and its applications in protein structural class prediction Warmuth MK, Kuzmin D. Randomized online PCA algorithms with regret bounds that are logarithmic in the dimension Heuristic molecular lipophilicity potential (HMLP): a 2D-QSAR study to LADH of molecular family pyrazole and derivatives Heuristic molecular lipophilicity potential (HMLP): lipophilicity and hydrophilicity of amino acid side chains AAindex: amino acid index database, progress report AAindex: amino acid index database Predicting membrane protein types by the LLDA algorithm iCTX-type: a sequence-based predictor for identifying the types of conotoxins in targeting ion channels Some remarks on protein attribute prediction and pseudo amino acid composition iSNO-AAPair: incorporating amino acid pairwise coupling into PseAAC for predicting cysteine S-nitrosylation sites in proteins iRSpot-PseDNC: identify recombination spots with pseudo dinucleotide composition iPro54-PseKNC: a sequencebased predictor for identifying sigma-54 promoters in prokaryote with pseudo k-tuple nucleotide composition iEnhancer-2L: a two-layer predictor for identifying enhancers and their strength by pseudo k-tuple nucleotide composition iSuc-PseOpt: identifying lysine succinylation sites in proteins by incorporating sequencecoupling effects into pseudo components and optimizing imbalanced training dataset pRNAm-PC: predicting N-methyladenosine sites in RNA sequences via physical-chemical properties pSuc-Lys: predict lysine succinylation sites in proteins with PseAAC and ensemble random forest approach Predicting antimicrobial peptides with improved accuracy by incorporating the compositional, physico-chemical and structural features into Chou's general iATC-mISF: a multi-label classifier for predicting the classes of anatomical therapeutic chemicals Unb-DPC: identify mycobacterial membrane protein types by incorporating un-biased dipeptide composition into Chou's general PseAAC iRSpot-EL: identify recombination spots with an ensemble learning approach OOgenesis_Pred: a sequence-based method for predicting oogenesis proteins by six different modes of Chou's pseudo amino acid composition iATC-mHyb: a hybrid multilabel classifier for predicting the classification of anatomical therapeutic chemicals iRNAm5C-PseDNC: identifying RNA 5-methylcytosine sites by incorporating physicalchemical properties into pseudo dinucleotide composition Prediction of the aquatic toxicity of aromatic compounds to tetrahymena pyriformis through support vector regression 2L-piRNA: a two-layer ensemble classifier for identifying piwi-interacting RNAs and their function Frequency of HLA-B alleles in a Caucasoid population determined by a two-stage PCR-SSOP typing strategy HLA-DPB1 DNA polymorphism in the Swiss population: linkage disequilibrium with other HLA loci and population genetic affinities Influenza virus-specific cytotoxic T lymphocytes recognize HLAmolecules. Blocking by monoclonal anti-HLA antibodies Tumorspecific lysis of human renal cell carcinomas by tumorinfiltrating lymphocytes. I. HLA-A2-restricted recognition of autologous and allogeneic tumor lines Induction of tumor-reactive CTL from peripheral blood and tumor-infiltrating lymphocytes of melanoma patients by in vitro stimulation with an immunodominant peptide of the human melanoma antigen MART-1 Identification of new HER2/neu-derived peptide epitopes that can elicit specific CTL against autologous and allogeneic carcinomas and melanomas Identification of a shared HLA-A*0201-restricted T-cell epitope from the melanoma antigen tyrosinase-related protein 2 (TRP2) Breast and ovarian cancer-specific cytotoxic T lymphocytes recognize the same HER2/ neu-derived peptide Conformational parameters for amino acids in helical, beta-sheet, and random coil regions calculated from proteins Mapping hydrophobicity on the protein molecular surface at atomlevel resolution Prediction of protein structural classes Toward the quantitative prediction of T-cell epitopes: coMFA and coMSIA studies of peptides with affinity for the class I MHC molecule HLA-A*0201 Binding of a peptide antigen to multiple HLA alleles allows definition of an A2-like supertype Identification of subdominant CTL epitopes of the GP100 melanoma-associated tumor antigen by primary in vitro immunization with peptidepulsed dendritic cells Comparison of cytotoxic T lymphocyte responses induced by peptide or DNA immunization: implications on immunogenicity and immunodominance Recognition of multiple epitopes in the www.impactjournals.com/oncotarget human melanoma antigen gp100 by tumor-infiltrating T lymphocytes associated with in vivo tumor regression Peptide binding to the most frequent HLA-A class I alleles measured by quantitative molecular binding assays The relationship between class I binding affinity and immunogenicity of potential cytotoxic T cell epitopes Prominent role of secondary anchor residues in peptide binding to HLA-A2.1 molecules Evaluation of MHC class I peptide binding prediction servers: applications for vaccine research Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins Comparative Molecular Field Analysis (CoMFA) and Comparative Molecular Similarity Indices Analysis (CoMSIA) studies on α(1A)-adrenergic receptor antagonists based on pharmacophore molecular alignment Comparative molecular similarity index analysis (CoMSIA) to study hydrogen-bonding properties and to score combinatorial libraries Molecular similarity indices in a comparative analysis (CoMSIA) of drug molecules to correlate and predict their biological activity Recent progresses in protein subcellular location prediction Cell-PLoc 2.0: an improved package of webservers for predicting subcellular localization of proteins in various organisms Prediction of protein cellular attributes using pseudo amino acid composition Molecular science for drug development and biomedicine An unprecedented revolution in medicinal chemistry driven by the progress of biological science Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology Pseudo nucleotide composition or PseKNC: an effective formulation for analyzing genomic sequences Pse-in-One 2.0: an improved package of web servers for generating various modes of pseudo components of DNA, RNA, and protein sequences Pse-in-One: a web server for generating various modes of pseudo components of DNA, RNA, and protein sequences An eigenvalue-eigenvector approach to predicting protein folding types Eigenvalue computation in the 20th century On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations Multiple field three dimensional quantitative structure-activity relationship (MF-3D-QSAR) Fragment-based quantitative structure-activity relationship (FB-QSAR) for fragment-based drug design Recent advances in developing web-servers for predicting protein attributes iCar-PseCp: identify carbonylation sites in proteins by Monto Carlo sampling and incorporating sequence coupled effects into general PseAAC iHyd-PseCp: identify hydroxyproline and hydroxylysine in proteins by incorporating sequence-coupled effects into general PseAAC iOri-Human: identify human origin of replication by incorporating dinucleotide physicochemical properties into pseudo nucleotide composition iRNA-AI: identifying the adenosine to inosine editing sites in RNA sequences Pse-Analysis: a python package for DNA/RNA and protein/peptide sequence analysis based on pseudo components and kernel methods iACP: a sequence-based tool for identifying anticancer peptides iPhos-PseEn: identifying phosphorylation sites in proteins by fusing different pseudo components into an ensemble classifier iROS-gPseKNC: predicting replication origin sites in DNA by incorporating dinucleotide position-specific propensity into general pseudo nucleotide composition iPGK-PseAAC: identify lysine phosphoglycerylation sites in proteins by incorporating four different tiers of amino acid pairwise coupling information into the general PseAAC pSumo-CD: predicting sumoylation sites in proteins with covariance discriminant algorithm by incorporating sequence-coupled effects into general PseAAC identifying DNase I hypersensitive sites by fusing three different modes of pseudo nucleotide composition into an en-semble learning framework iPTM-mLys: identifying multiple lysine PTM sites and their different types iPreny-PseAAC: identify C-terminal cysteine prenylation sites in proteins by incorporating two tiers of sequence couplings into PseAAC iRNA-2methyl: identify RNA 2'-O-methylation sites by incorporating sequence-coupled effects into general PseKNC and ensemble classifier pLoc-mPlant: predict subcellular localization of multi-location plant proteins via incorporating the optimal GO information into general PseAAC Impacts of bioinformatics to medicinal chemistry The authors declare no conflicting interest.