Chapter 1.5. Data mining. II. Prediction of protein structure and optimization of protein crystallizability

Various strategies for overcoming the size and time-scale limitations have been proposed, including the development of novel coarse-grained (CG) methods and implicit solvent models. The CG models reduce the complexity of each amino acid by representing them by fewer numbers of pseudo-atoms (nodes) than the actual number of atoms, which allows a significant decrease in the number of degrees of freedom (Liwo et al., 2011; Kolinski, 2004; Kmiecik et al., 2016; Wabik et al., 2013), especially since coarse graining may typically be performed at the level of one geometric point per amino acid. Despite the fact that the improved computational performance of CG models comes at the cost of loss of some structural accuracy, CG models are still powerful and successful tools for de novo structure prediction (Kmiecik et al., 2016).

There are various alternative ways to improve conformational sampling efficiency: one can use enhanced sampling techniques such as replica-exchange (Hansmann, 1997; Kouza & Hansmann, 2011) and umbrella sampling (Boczko & Brooks, 1995), or reduce the conformational space further by implementing a fragment-based approach (Zhang, 2008; Simons et al., 1997; Blaszczyk et al., 2013, 2016).

The fragment-based approach to predicting protein structures was adopted early on in the Rosetta suite (Simons et al., 1997). Rosetta creates a library of possible conformations of short fragments of the target protein sequence. Those short fragments, derived from a set of known structures, are assembled by a Monte Carlo (MC) minimization together with a scoring function rewarding native-like properties. The success of Rosetta in structure and function predictions shows the importance of the information in structural fragments for theoretical predictive studies of proteins. Newer versions using sequence correlation information predict contacts within structures with great success (Ovchinnikov et al., 2017).

It is widely assumed that a protein's three-dimensional structure is encoded in its sequence of amino acids. Although there are counterexamples, as sometimes one mutation in a sequence can lead to very different protein folds and functions (Alexander et al., 2009; Kouza & Hansmann, 2012), a sufficiently high similarity between protein sequences generally implies similarity between protein structures. This observation is the basis of many structure prediction algorithms that use sequence similarity searches to identify homologous proteins. While protein crystallography may be the method of choice for structure determination, it is too expensive and too slow to deal with the exponentially growing number of protein sequences being produced by high-throughput genome sequencing. Computational methods are thus going to be needed to solve this problem (Baker & Sali, 2001). One of the key steps in the prediction of tertiary structures of proteins is the prediction of secondary structures. Typically, in ab initio methods this is the first step. Prediction of secondary structure is a one-dimensional problem, which is much easier than three-dimensional structure prediction, and has been a subject of interest for a long time. However, de novo protein structure prediction (and secondary structure prediction) still remains a challenging problem in molecular biology.

Recent advances in technology have made storing huge data sets extremely cheap. More than 95% of protein structural data has been generated and stored in the last decade and there seems to be no end to new data generation. Data-mining techniques for extracting useful information from huge data sets is a maturing field. To analyse and process these enormous amounts of biological data, new `big data' methodologies have been developed. They are based on mapping, searching and analysing for patterns in sequences and data mining of three-dimensional structures. Mining for information in biological databases involves various forms of data analysis such as clustering, sequence homology searches, structure homology searches, examination of statistical significance and so on.

The most basic data-mining tool in biology is the basic local alignment search tool (BLAST; Altschul et al., 1990) for examining a new nucleic acid or protein sequence. BLAST compares the new sequence to all sequences in the PDB database to find sequences that are the most similar to a query sequence. This is normally the way genes and proteins are assigned functions.

Data-mining tasks can be descriptive, by uncovering patterns and relationships in the available data, or predictive, based on models derived from these data. Most popular and frequently used are automated data-mining tools that employ sophisticated algorithms based on statistics or machine learning to discover hidden patterns in biological data. Data mining can be applied to a variety of biological problems such as analysis of protein interactions, finding homologous sequences or homologous structures, multiple sequence alignment, construction of phylogenetic trees, genomic sequence analysis, gene finding, gene mapping, gene expression data analysis, drug discovery etc. All of these different problems can be studied by using various data-mining tools and techniques.

In this chapter, first we discuss the application of data mining to the problem of protein secondary structure prediction. We describe an original method (fragment data mining, FDM), which has proven to be a good predictor of the secondary structure of proteins. The FDM method mines structural segments in the Protein Data Bank (PDB) and utilizes structural information from the matching sequence fragments, with the aim of improving predictions of secondary structure. Next, we briefly describe the consensus data mining (CDM) method, which combines FDM with the classical GOR V secondary structure prediction method, which is based on information theory and Bayesian statistics, coupled with evolutionary information from multiple sequence alignments. Finally, we discuss the application of the data-mining approach to the problem of optimizing protein crystallization conditions.

Taking inspiration from the approach implemented in Rosetta, we developed the FDM method. Because a sufficiently high similarity between protein sequences implies similarity between protein structures and conserved local motifs may assume a similar shape, we searched for a local alignment method to obtain structure information for predicting the secondary structure of query sequences. In the FDM method, BLAST is applied to query all the structures from the PDB. Then, evolutionary information is introduced to the prediction of secondary structure by analysing the fragments of the alignments belonging to proteins in the PDB. In the prediction and evaluation part, each query sequence from the data set is subjected to the following procedure.

The procedure starts from the weights assignment to the matching segments obtained from BLAST. In this step various types of parameters are considered, including different substitution matrices, similarity/identity cutoffs, degree of exposure to solvent of residues, and protein classification and sizes. In the second step, we calculate normalized scores for each residue. This is followed by prediction of the secondary structure for that residue according to the normalized scores. It should be noted that in order to predict the secondary structure according to the normalized scores of residues, we applied two different strategies. The first strategy is to choose the highest-scoring structure class as the prediction, and the second one is to use machine-learning approaches to choose a classification based on training. Finally, we calculate Q₃ (the sum of the fractions of successful predictions for α-helix, β-strand and coil) and the Matthews' correlation coefficients.

Fig. 1 shows a scheme of the multistage procedure employed in the fragment database mining (FDM) method.

Figure 1 Multistage procedure for the prediction of secondary structure using the fragment database mining (FDM) method.

For training we used the CB513 benchmarked data set of 513 non-redundant domains developed by Cuff & Barton (2000, 1999). Local sequence alignments were generated by using BLAST with different substitution matrices (Henikoff & Henikoff, 1992), including BLOSUM-45, BLOSUM-62, BLOSUM-80, PAM-30 and PAM-70.

Unlike in DSSP (Kabsch & Sander, 1983), where an eight-letter structural alphabet is used, we instead used a reduced three-letter structural alphabet of secondary structures: α-helix (H), extended (β-sheet) (E) and coil (C). H, E and C in our three-letter structural code correspond to helices (H, G and I), strands and bridges (E and B) together with turn, bend and coil (T, S and C) in the DSSP code, respectively.

For weight assignments we defined identity scores and their powers (id^c, where c is a positive real number) as the weights of matching segments. Here id is the ratio of the number of exact matches of residues to the total number of residues in the matching segment. Weights were then adjusted to obtain the best match. This is illustrated in Fig. 2. At each position, the predicted secondary structure is determined by the secondary structures of the matches at that position. Each match is assigned a weight according to the similarity or identity score of the alignment from BLAST. At each position, the weights are normalized, and the normalized scores for each position within each of the secondary structure states are calculated.

Figure 2 An example showing a query sequence and its matching segments based on sequence matches (sequences not shown). The matching segments are expressed as secondary structure elements. The weights are shown for each segment.

We defined s(H, i) as the normalized score for position i to be in the state H,Here w(H, i) is the weight for one matching segment with a residue at the ith position in a helix, and w(E, i) and w(C, i) are similarly defined. For example s(H, 2) = 0.2/(0.1 + 0.2 + 0.4) = 0.29, and s(E, 5) = (0.1 + 0.2)/(0.1 + 0.2 + 0.3 + 0.4) = 0.3. The secondary structure state having the highest score is chosen as the final prediction result for a given position in the sequence. For the ith position of a query sequence, we have three normalized scores for each secondary structure state: s(H, i), s(E, i) and s(C, i). In our prediction scheme, we always choose the highest score among these three to determine the secondary structure prediction at the ith position.

We have used two different substitution matrices: PAM and BLOSUM. The PAM (percent accepted mutation) matrix was introduced by Dayhoff et al. (1978) to quantify the amount of evolutionary change in a protein sequence, based on observation of how often different amino acids replace other amino acids during evolution. BLOSUM (blocks substitution matrix) was introduced by Henikoff & Henikoff (1992) to obtain a better measure of differences between two proteins, specifically for distantly related ones. The BLOSUM matrix is derived from observations of the frequencies of substitutions in blocks of local alignments in related proteins. We used several versions of these matrices (BLOSUM-45, BLOSUM-62, BLOSUM-80, PAM-30 and PAM-70) in BLAST. For example, PAM-30 refers to 30 mutations per 100 amino acids of sequence, and BLOSUM-62 means that it was derived from sequence blocks clustered at the 62% identity level.

We tried different combinations of matrices and identity powers. The best results were obtained by using BLOSUM-45 and id³ as the weight assignment method. [See Fig. 3 in Cheng et al. (2005) for prediction accuracies computed for the CB513 data set using five different substitution matrices and several powers of identity scores.]

To test the concept of fragment assembly for secondary structure prediction, we set limits to the extent of similarity or identity for fragments to be included by using different cutoffs: 99, 90, 80, 70 and 60% of similarity or identity scores. The matches having similarity or identity scores higher than a cutoff were eliminated from the matching lists of segments that are used for calculating normalized scores. Fig. 3 shows results obtained for BLOSUM-45 and id³. We observe a constant decrease in Q₃ with a decrease in the cutoff.

Figure 3 Prediction accuracy Q3 using different identity cutoffs.

Matches with the highest identity scores (greater than the id cutoff) were filtered out, but `reasonably' high id matches were kept. We define the `reasonable' high-id matches to be those that have relatively high identity scores (greater than the id cutoff) and are not too short (>5 residues) or too long (less than 90 or 95% of the length of the query sequence).

The prediction accuracies Q₃ were compared for three cases, always using BLOSUM-45 with an id cutoff of 0.90. In case (1), all high-id matches (matches with identity scores higher than the identity cutoff) are filtered out. This case is used as a control. In case (2), sequences with identity scores greater than the id cutoff (0.90 here), lengths longer than 5 residues and lengths less than 90% of query sequence are retained. In case (3), sequences with identity scores greater than 90%, lengths longer than 5 residues and lengths below 95% of query sequence are kept.

Table 1 shows the accuracies for all three cases. We observe again that id³ gives the best accuracy.

It is interesting to include approximate tertiary information in the computations, as recent studies (Adamczak et al., 2005; Faraggi et al., 2012) have demonstrated that secondary structure prediction can be improved by taking into account solvent-exposed surface area. To do so, we use the software NACCESS (Hubbard & Thornton, 1993), which allows calculation of the degree of exposure to solvent for all available PDB sequences. The idea was to differentiate buried and exposed residues by assigning them different weights according to the following rule. If the computed relative accessibility of a residue is less than 5.0%, then it is regarded as a buried one, while the residue is considered as exposed if its computed relative accessibility is greater than 40.0%. The others are regarded as intermediate cases. Buried residues were weighted most heavily. We made the following linear changes to the weights of residues: if a residue is buried, its weight is doubled; if it is intermediate, the original weight is multiplied by 1.5; if exposed, the weight remains unchanged. However, our results obtained for weights id³, id cutoff = 0.99 and the BLOSUM-45 substitution matrix show that differentiation of buried and exposed residues does not have a significant effect. Alternatively, computations could use a fast ASA (accessible surface area) predictor (Faraggi et al., 2014, 2017) to estimate the benefits of including solvent-exposed surface area possibly to increase the accuracy of secondary structure prediction, but the present results do not make this look very promising.

In order to evaluate whether the protein size leads to any changes in the accuracy of prediction, we divided all proteins from the CB513 data set into four groups according to the sequence length: very small (n ≤ 100 residues), small (100 < n ≤ 200 residues), large (200 < n ≤ 300 residues) and very large (n > 300 residues). The very small, small, large and very large groups contained 154, 216, 84 and 58 sequences, respectively. As in the previous case, we used the BLOSUM-45 substitution matrix with weight function id³. No optimization was applied. The prediction accuracies for proteins of different size are presented in Table 2. With increasing sizes, we noted a small increase in accuracy from 0.911 for the smallest group to 0.948 for the large category.

The success of FDM largely depends on the availability of fragments similar to the target sequence. In practice, however, the availability of similar sequences can vary significantly. In order to analyse the relationship between the performance of CDM and the sequence similarity of fragments, we methodically excluded fragment alignments with sequence identities above a certain limit and called this limit the upper sequence identity limit. The upper sequence identity limit is not an additional parameter in the CDM method; these results demonstrate what expected results would be in the absence of fragments with similarities above the sequence identity limit.

The performance of all secondary structure prediction methods can be improved with multiple sequence alignments: the GOR V method tested with the full jack-knife methodology (a popular leave-one-out resampling technique) yields an accuracy of 73.5% when multiple sequence alignments (MSAs) are included; otherwise the accuracy is about 10% less.

One of the significant advantages of FDM is its applicability to various evolutionary problems, because the algorithm does not rely exclusively on the sequences with the highest sequence similarity, but assigns weights to BLAST-aligned sequences that apparently capture divergent evolutionary relationships. As a result, CDM, which incorporates FDM, can be successfully used even when there is a range of sequence similarities for the BLAST-identified sequences (Kandoi et al., 2017).

Recent efforts to improve protein sequence matching are yielding significant gains (Jia & Jernigan, 2021) by incorporating structural information into the sequence matching. This is done by utilizing correlated sequence pairs that are in contact in protein structures to increase the number of allowed substitutions. By using this approach, the sequence fragment library is significantly expanded and therefore the results for the methods discussed above should be significantly improved.

Note that the results shown here are not conventional machine-learning results in the sense of using a training set and a test set, but are just from trials to investigate how successful the variations on the approach can be, all on the same set of known experimental structures. The next secondary structure prediction method discussed below is a machine-learning method, having training and test sets of data.

The first five steps predict torsion angles (both φ and ψ). SPINE-X first generates a position-specific scoring matrix (PSSM) using PSIBLAST and seven representative physical parameters (PPs): a steric parameter (graph shape index), hydro­phobicity, volume, polarizability, isoelectric point, helix probability and sheet probability. In the first step, a neural network predicts a secondary structure SS0 employing PSSM and PP as input. In the second step, another neural network is built to predict residue solvent accessibility (RSA) with PSSM, PP and predicted SS0 as input. Then, in the third step, the predicted RSA and SS0 together with PSSM and PP are used to predict the torsion angles. In the fourth step, a new round of secondary structure prediction (SS1) is performed based on the previous predictions. Then, in the fifth round, new torsion-angle predictions are performed based on predictions from previous iterations. In the sixth (final) step, a neural network is trained to predict secondary structure using PSSM, PP and predicted values from the first five rounds.

In each step, the general architecture of the neural networks is the same, consisting of two hidden layers with 101 hidden nodes. Training and initial testing for all neural networks were performed on the SPINE data set of 2640 proteins from the PDB and on its subset of 2479 proteins with lengths less than 500 residues.

For secondary structures the accuracy is between 81 and 84% depending on the data set and choice of tests. The Pearson correlation coefficient for accessible surface area predictions is 0.75 and the mean absolute error for the φ and ψ dihedral angles are 20° and 33°, respectively. All details are available in Faraggi et al. (2012) and Faraggi & Kloczkowski (2017).

The ultimate goal of protein science is to determine the three-dimensional structures of all proteins and to determine their functions and interactions with all other proteins and ligands. X-ray crystallography (Kendrew et al., 1960) and nuclear magnetic resonance (NMR) spectroscopy (Bax & Tjandra, 1997) are the main techniques that are used to determine three-dimensional protein structures, although applications of cryo-EM are growing rapidly. NMR spectroscopy allows proteins to be studied in solutions rather than in crystals; however, this method has disadvantages in that it is limited to relatively small and medium-sized proteins. For larger proteins, the NMR spectrum becomes highly crowded with many overlapping peaks and therefore very difficult to interpret.

Ninety per cent of the experimental structures in the PDB were determined by X-ray diffraction. X-ray crystallography is based on diffraction of an X-ray beam by a crystal lattice, as discovered in 1911 by Max von Laue. The incoming beam is scattered, and the directions and intensities of such reflections depend on the types and distribution of the atoms within the crystal. In 1912 Bragg showed that X-ray scattering can be used for structure determination. The first crystal structures of myoglobin and haemoglobin were determined by John Kendrew in 1957 and Max Perutz in 1959, respectively (for which they shared the Nobel Prize in Chemistry in 1962), and since then the amount of structural data deposited in the PDB has been growing at a rapid pace. While it took Max Perutz 22 years to determine the crystal structure of haemoglobin, the structure deposition rate has been increasing from about one structure per year in the 1960s to about 1 structure per hour in the last decade. Now (as of March 2022), the number of structures solved experimentally and deposited in the Protein Data Bank (Berman et al., 2000) is ∼188 000, while the number of known protein sequences in March 2022 in UniProtKB/TrEMBL (The UniProt Consortium, 2015) is 230 million. The number of protein sequences continues to grow at a faster rate than the number of protein structures deposited in the PDB. Note that in 2014 the number of solved molecules and known protein sequences were ∼114 000 and ∼80 million, respectively (The UniProt Consortium, 2015).

There are many steps on the path to crystal structure determination. Most of the stages have been improved, optimized and automated by using robots, faster data-collection devices, better fitting procedures, remote control of the crystallization process etc. However, the main bottleneck in the crystallization process that still needs to be improved is the way in which the conditions for crystallization are chosen in order to obtain a sufficiently large high-quality crystal for X-ray diffraction. Unsuccessful attempts (the success rate is generally well below 10%) increase the cost of the crystallization process by up to 70%, as reported by the Joint Center for Structural Genomics (Jahandideh et al., 2014). A better understanding of all the factors governing the crystallization process has attracted intensive experimental and theoretical interest, as this could provide insight into the critical relevant variables, i.e., the crystallization space, for optimizing crystal size and morphology. The ability to predict the critical variables (the concentration and the nature of the protein, salt concentration, type of salt, pH, buffer, additives, temperature, precipitant type and concentration, etc.) for a given protein could significantly speed up the structure determination process. Recently, data-mining approaches have become powerful tools for characterizing the crystallization space that can help to control and improve the crystallization process.

Publicly available databases such as the Biological Macromolecule Crystallization Database (the BMCD) (Gilliland et al., 1996, 2002; Tung & Gallagher, 2009) and the PDB provide crystallization information for macromolecular structures. The PDB is not only the largest database of experimentally determined three-dimensional structures of biological macro­molecules, but also contains a wealth of information about the crystallization process. There have been attempts to broaden our understanding of the crystallization process by analysing this information (Peat et al., 2005; Kirkwood et al., 2015; Pérez-Priede & García-Granda, 2017).

The BMCD (Gilliland et al., 1996, 2002) includes the protein name, protein concentration, crystallization precipitant, pH, temperature, unit cell and resolution, which are the parameters describing how a given protein has been successfully crystallized. In version 7.0 of the BMCD, the content has been expanded to include 99 211 crystal entries and the following features have been added: macromolecule sequence, which can enable the use of more elaborate analysis of relations among protein properties, crystal-growth conditions, and the geometric and diffraction properties of the crystals (Tung & Gallagher, 2009). The BMCD is available as a server at http://bmcd.ibbr.umd.edu/ .

It should be noted that there are two main limitations to the BMCD database that restrict its application for unbiased data mining. First, it only includes positive data by collecting data from successful experiments. Negative results are not reported in this database and only rarely reported in the literature in general (Newman et al., 2012). Such negative results could be a key to deeper understanding of the crystallization process. A second limitation is that preparation methods for different entries in the database are not always fully described and can be significantly dissimilar. However, despite these limitations, the BMCD database has been used for successful crystallization outcomes (Schiefner et al., 2015; García-Fernández et al., 2012).

In order to improve crystallization success rates, there have been attempts to obtain the minimal sets of conditions that can be used to crystallize most proteins in a given data set. At the Joint Center for Structural Genomics (the JCSG), by mining a data set of 539 T. maritime proteins, the ten most effective conditions were identified to crystallize 196 proteins, while 108 best conditions led to 465 successful crystallization outcomes. The authors identified 67 conditions that were the most productive in promoting protein crystallization, and referred to these as the core screen. Together with the next 29 most effective conditions these formed the expanded core screen, which has been widely used for initial crystallization trials (Page & Stevens, 2004).

Kimber et al. (2003) mined a data set of 775 proteins with 48 conditions using different sample preparation processes and crystallization conditions. It was demonstrated that proteins from different species typically exhibit different crystallization behaviour. A minimal screen with six conditions produced 205 crystals, while an extended screen with 24 conditions led to 318 crystallized proteins out of the 338 trial proteins in the set. These results also show that many of the conditions used in sparse matrix screens are redundant, and that current screens are not sufficient to crystallize all proteins, since many proteins failed to crystallize at all (Kimber et al., 2003; Page & Stevens, 2004).

Segelke (2001) and Rupp (2003) suggested a random sampling of crystallization space and DeLucas et al. (2003) applied incomplete factorial screens to streamline the crystallization procedure. For some specific target proteins, those strategies are more appropriate and efficient to identify crystallization conditions. Oldfield (2001) applied data mining of protein fragments for molecular replacement models that are used for solving the phase problem in X-ray crystallography.

Babnigg & Joachimiak (2010) analysed protein properties for more than 1300 proteins that are well expressed but insoluble, and for ∼720 unique proteins for which structures had been solved by X-ray diffraction. They showed that a protein's isoelectric point and grand average hydro­pathy (GRAVY) correlate with its propensity to crystallize. Additional physiochemical properties of amino acids from the AAindex database (Kawashima et al., 2008) were considered and used for data-mining purposes. The set of attributes that were most strongly correlated with protein crystallization propensity was identified and has been successfully used by incorporation into a support vector machine (SVM) classifier. Using the proposed SVM method, conditions for insoluble proteins and for proteins with solved structures, deposited in the PDB, were predicted at 56% and 75% accuracies, respectively.