The fragfinder program has not had extensive testing—we selected an algorithm and HMM based on intuition and minuscule tests on one target sequence. The undertaker program can be used to evaluate the fragments as fragments (looking at the histogram of RMSD, or RMSD normalized by length, for example), and we can come up with various figures of merit for fragment libraries, allowing us to compare different fragment finding methods.

Sol Katzman wrote routines for undertaker to filter fragfinder output and produce Rosetta fragment libraries. Some of his work provides more measures for fragment-library quality. It also opens up the possibility of comparing fragfinder-generated libraries with the libraries that Rosetta normally uses.

Evaluating quality measures

Some of these measures are known to have problems, particularly for structures that are not very similar to the real structure. For example, RMSD and RMSD_CA are both very sensitive to the worst-done part of the structure, and interpretation of the result is very dependent on the length of the chain. GDT (and smooth_GDT) evaluates only how much of the chain is correct, giving no penalty for grossly wrong parts. All of these superposition-based measures tend to over-reward structures that have been compacted too much—they have no penalty for unreasonably dense proteins. In some cases, the point model (all atoms at (0,0,0)) scores reasonably well.

The clens function was devised to try to handle the problems of the other measures, be rewarding for correct contacts and penalizing for incorrect ones. In preliminary testing using models generated by undertaker, clens and GDT were very highly correlated with each other, indicating that they may be nearly interchangeable in their ability to recognize decent models. We'd like to do more extensive testing of "clens", using all the models available from previous runs of CASP. In particular, we'd like to see if we can find any outliers, where GDT and clens evaluate models quite differently. We could then examine such outliers, and see which metric does a better job of capturing the quality of the model.

It is also possible to modify the clens definition in undertaker, to parameterize the definition of contact or change the scoring function for missing and extra contacts.

Rotamer optimization

A possible approach is to read a pdb file, randomize the rotamers, optimize with undertaker with only the rotamer-changing operators turned on, reoptimize with scwrl (callable from inside undertaker), then reoptimize with undertaker, reporting the three costs and all-atom rmsd for each of the conformations. One could also gather all-atom rmsd for a huge number of sidechain conformations and try to do linear regression to set weights for the cost function, but this is likely to have serious problems due to correlation between the components of the cost function.

Tweaking undertaker's cost function

I did some preliminary work on this, and (as expected) linear regression works terribly when the cost functions are highly correlated. We'll probably need to do some dimension reduction first (perhaps with principal components analysis?) before linear regression stands a chance of working.

A minor programming project that may need to be done before we can try tweaking the cost functions: we need to clean up cost functions to handle incomplete conformations. When reading in true, but incomplete, conformations (the usual case for PDB structures) the cost functions often report the conformations as terrible, since the missing atoms are all at location (0,0,0). The cost functions need to be modified to be aware of the "contains_atom()" member of the Conformation class. Some cost function components have already been fixed, but not all.

Protein design by back propagation

More random thoughts on using back-propagation. There are lots of variants that could be tried in the sampling of sequences from the back-propagation. For example, one could resample the sequence points where the back-propagated derivative for the current residue is negative rather than resampling all points. Instead of sampling from the profile one could take the most probable residue of the profile. After generating a few hundred sequences in this way from random starting points, one could make a profile from the sequences, and score each sequence with that profile. Further tweaking could take just the best 100 or so sequences and the profile generated from them. Each could be used with that profile as a starting point for the back-propagation and resample process. A tertiary scoring system (such as Rosetta's) could be used to pick out the best sequences from the list, and use the resulting sequences and profile to do more sampling.

Another approach to using neural nets would be to build a neural net whose input was many structural properties and whose output was a sequence. Such a neural net could be trained in the same manner as current nets for prediction (with inputs and outputs swapped), and would provide a probability distribution for amino acids at each position. This would require a few changes to the neural net code (mainly providing an input layer that accepts several alphabets, instead of just one alphabet), and would more directly provide probability distributions to sample sequences from. This method is less likely to provide a variety of very different solutions to a protein design problem, though, so I have not pursued it.

One protein that it may be interesting to apply the back-propagation plus Rosetta approach to is the "miniAGRP" peptide, a small cysteine knot that Glen Millhauser's group is interested in. They would like to remove the disulfide bridges that currently hold the structure together, without greatly disrupting the structure. Certain of the residues are believed to be functionally important, and should not be touched in the redesign. Millhauser's group may be interested in synthesizing a design, if a better one can be found than the design that Craig Lowe did (using just Rosetta) last year. For more information, see miniAGRP paying particular attention to the README file and the design/ subdirectory.

In addition to that particular design, it would be good to experiment with various networks and local structure alphabets, so see which combinations result in the best designs. Since we do not have the time nor the facilities for fabricating and testing large numbers of proteins, we will have to rely mainly on computational methods. One simple test is to look at the percent sequence recovery: given just the backbone for a real protein, try to design amino acids sequences that would fold to that conformation. The closer the designs are to the real protein, the more convincing the design method. Such a study should be done on hundreds of different backbones from radically different structures.

Projects in this section generally require some programming, database, or web-design skills.

SAM-T05 Web site

The SAM-T05 web site was created last spring and summer, but it is still in beta release form, and we'd like to make it public before CASP starts. Some things that need doing include the following:

• Fix the muscle profile-profile alignment script, which keeps crashing, probably because of bugs in sam_t04_conf.pm
• Add target05 multiple alignment (currently does target2k and target04). This requires training new neural networks for t05 alignments and building a somewhat more complete t05 template library.
• Add the new hbond alphabets to the secondary structure predictions.
• Change the weighting for different multi-track HMMs, based on fold-recognition tests done by Grant (and possibly others), as we are currently putting too much weight on less accurate predictors.
• Reduce the number of different target-model scoring runs, to spped things up. Try to merge the calibration and scoring runs, if SAM now supports that.
• minor updates in citations and error messages:
  • In RunProgWeb.pm, change error messages about Target99 to SAM_T05, if relevant. (Or eliminate RunProgWeb.pm) Also in CheckSeqs.pm and ScoreLibT99Web.pm
  • In IndexPage.pm (maybe elsewhere) eliminate references to HMM-apps.
  • Update citations in annotate_target_scores and annotate_template_scores (possibly elsewhere) to include the CASP6 paper and to explain when to cite each paper.
• Debug gromacs options (running on desktops, but not farm-cluster machines).
• Separate SAM_T05 and SAM_T05test directories so that dev2live can be used effectively. Look for any hardwired references to SAM_T05 directory.
• Change constraint generation scripts to give higher weights to helix and strand constraints (provide external multiplier).
• Automatically generate sheet constraints from the better alignments using an undertaker script.
• Make sure that the frozen version of programs is used, not the /projects/compbio/bin version. This is occasionally a problem with nested calls (such as prettyalign from a2m2html).
• The LIBSIZE macro should be different for t2k, t04, and t05 libraries, but currently is not.
• Make sure that the target name cleanup leaves only alphanumerics, minus sign, and underbar in the name, as the target name gets used for generating file names. [I believe this has been done now.]
• The recursive make for creating alignments is far too slow. We need to have a more efficient, but still easily modified, way to generate the 30 or so alignments per template that we need.
• The single-track-alignment and muscle-profile-profile alignments fail when something is in the t2k template library, but not the t04 template library. The errors are currently ignored, but the checks should be refined when a better alignment script is created. The problem will only increase when we add the t05 template library, since it will have less redundancy than the older libraries.
• Richard suggested reordering the .html file by using includes, then generating the included .html-part files in a different order than they occur in the master file.
• Maybe provide a link for a jmol viewer for the pdb files (assuming jmol can handle gzipped pdb files).
• The text for the pictures should be different for the pictures from undertaker-align and the pictures from the try1-opt2 model (this should be fixed in add_jpeg_views_html). There should also be an explanation for the multiple models in undertaker-align.
• For CASP, we will need to create files in the official "TS" format, which just requires a minor wrapper around the PDB files. We'll probably want to submit the try1-opt2 model as our first model, then the first 4 models of undertaker-align as our next 5. The necessary scripts are available in the casp6 directory, but their use has not been fully automated as they relied on manual editing of superimpose-best.under to select which 5 models to submit. This should be fully automated. There may still be some remnants of the older web service (which submitted alignments, not 3d models, to CASP), and these need to be cleaned out of SAM_T05.
• Some metaservers would like to have secondary structure (DSSP_EHL2 in CASP format) and the final models e-mailed, so these should be options on the form.

It may be worthwhile to do some other small genomes, such as the archeal genomes that Todd Lowe is building DNA chips for, or the bacterial genomes that Fitnat Yldiz and Karen Ottemann are interested in.

Index the yeast predictions

This could be done by adding a combining layer to the outputs of a normal neural net. For each output of the neural net we associate two parameters (center and spread). The predicted mean would be the neural-net weighted sum of the center values, and the predicted standard deviation would be the standard deviation of a mixture distribution, with mixture coefficients from the neural net.

The center and spread parameters could initially be set by partitioning the outputs in the training data into roughly equal weight discrete classes (as is now done) and computing a mean and standard deviation for each class. They could be tweaked a bit by a few passes of EM to get the output described as a mixture. After a few epochs of training the neural net, the center and spread parameters could be trained by back propagation also.

We would need to compare the encoding cost of the actual burial value using this continuous method (which requires some care in the discretizing, since our underlying data is usually integer-valued) versus the encoding cost of the discrete alphabet method, which assumes a flat distribution within each of the discrete bins.

Recoding neural net inputs

The resulting vectors can be clustered (with the square of the dot product as the appropriate similarity measure) to determine the most common directions. Alternatively, Principal Components Analysis (PCA) could be done to reduce the dimensionality. We can attempt to identify the properties these directions represent by looking at the correlation with the indices in the AAIndex database.

We can also compute the dot products of the profile with the centers of the most populated clusters and use them as inputs to the neural net, in addition to or in place of the profile input. (Similarly, the factors that come from PCA can be used for recoding.)

The above residue-at-a-time analysis does not get at the neighborhood patterns that the neural nets learn. People have tried to determine some of the patterns by looking at Fourier transforms and wavelet transforms of numeric properties (mainly hydrophobicity). We could look at each hidden unit on the first neural layer and see if its weight matrix can be roughly matched by taking a single direction vector and a neighbor weight function. That is, decompose
sum_{a,-2<=j<2} w_{a,j} Prob(amino acid a in position i+j)
into
sum_{-2<=j<=2} f(j) sum_a dir(a) Prob(amino acid a in position i+j)
or w_{a,j} = f(j) dir(a).
If there are hidden units that can be approximately described this way, it may be worth providing input features that correspond to this neighborhood property also.

The neighborhood analysis may be a bit limited in our current networks, since we have gone to using rather short windows (only 3 wide) in the first layer, but the residue-at-a-time recoding analysis may be more revealing than in networks with wider input windows.

Another possible residue-at-a-time recoding scheme would be to replace the 1-hot encoding of 20 amino acids for the guide sequence with a binary encoding of amino acids. There are some simple greedy algorithms and random search optimization algorithms that could be applied to find 5 or 6 different partitions of the amino acids so that each of the weight vectors could be well approximated by adding just 5 or 6 weights to the appropriate subsets of the amino acids.

Entropy profiles

To avoid artifacts, we should use our own testing methods and our own computation of entropy. Note that entropy is not a property of the protein structure (unlike most of our alphabets), but a property of a multiple alignment. This means that we do not have all our usual tools for evaluating an alphabet (in particular, conservation and predictability are not meaningful concepts), so we are limited to alignment tests and fold-recognition tests.

Even if entropy profiles alone are not very useful, they could make an interesting extra track in profile-profile HMMs.

Converting RDB files to HMMs

Although there is not much code needed for adding this functionality to the SAM package, a lot of code would need to be read in order to understand the internal data structures and parameter passing mechanisms of SAM.

Combining fold recognition results

One possible project is to incorporate the "product-of-pvalues" method into our automatic predictions. This requires designing a way to update the calibration parameters for the combining method automatically as new templates are added to the library, and converting the fold results back into selection of good templates. (This needn't be a weekly update, as the SCOP classification is only updated once or twice a year.)

Another possible project is to devise a different combining method. The "product-of-pvalues" method treats all templates in a fold class as equally informative and does not use information from templates in competing fold classes. One can envision combining methods such as logistic regression that could use the extra information to get better predictions. One danger is that multi-domain proteins correctly have multiple correct folds, so simple competition between fold classes is not quite the right model.

Handling multi-domain proteins

We have a script available for breaking up long proteins and doing separate predictions on each part (split-into-domains). The script can be used either manually to pull out a predicted subdomain for separate prediction or automatically to predict overlapping windows of some fixed length.

In CASP6, it was obvious that we had managed to put together some multi-domain predictions without explicit domain prediction. It is not clear to me that explicit domain prediction is as valuable as I thought in 2002. It is clear though, that a subdomain with many homologs can result in masking of other domains in the sequence, due to sequence-weighting artifacts (among other problems).

Projects involving multiple-domain proteins include

• Better annotation of hits in the web server (reporting only those SCOP domains that overlap significantly with the target, rather than all SCOP domains of template).
• Automatic parsing into subdomains, with separate predictions for (overlapping) subdomains.
• Automatic combining of secondary and tertiary predictions from subdomains, to guide tertiary prediction of whole chain. This is not as difficult as it sounds, as undertaker can take both constraints and alignments from subdomains without modification. There may be a little more difficulty in piecing together tertiary predictions, as separately predicted subdomains tend to be incompatible without some modification.

• One of the problems is that simple i.i.d. models of "random" sequences work ok for amino-acid sequences, but not for secondary structure sequences, which often have long runs of the same letter. We have a sophisticated sequence generator for amino-acid sequences released as open-source code: gen_sequence, but no good random sequence generator for secondary structure sequences. There are several ways to try to generate random secondary structure sequences: i.i.d. models, Dirichlet mixture for composition with i.i.d. from the composition (used in gen_sequence), first-order Markov chains, higher-order Markov chains, variable-order Markov chains, mixtures of Markov chains, hidden Markov models, Markov chains of segments with length models for each segment, and so forth. Implementing one or more of these methods and testing how well it can encode secondary structure sequences (over various alphabets) in a train/test experiment would be very useful. I'm particularly intrigued by the idea of generating sequences of segments, with a separate length model for each segment type, and have done a little preliminary work on the idea.
• Even if we get a good random sequence generator for the various local structure alphabets, we still have some trouble with the failure of the "reversibility assumption" for some alphabets that otherwise look quite useful, such as de Brevern's protein blocks. For multi-track HMMs using these alphabets, we need a different way of calibrating the E-values that does not use reverse-sequence null models. One possibility would be to fit a Gumbel extreme-value distribution to scores on randomly generated sequences (using a good generator for the local structure alphabet). Other parameterized distributions could also be explored, as there is some evidence that the Gumbel distribution does not have a fat enough tail.

Model drift often results in the HMM being better at recognizing some subfamily different from the one containing the target protein. For example, in one test, the SH3 domain for 2abl was given as a seed for the SAM-T02 web site, and the self-hit was about 80 down from the top (partly because 2abl is not in the template library). Blast finds several sequences containing exact matches to the 55-residue seed: 1oplA, 1opkA, 2abl, 1ju5C, 1abq, 1abo[AB], 1awo, and 1bbz[ACEG]. Four of these were in the template library (1aboB, 1awo, 1bbz[AE]), but the best scoring match using the template library was 1ng2A, which has the SH3 domain twice. If we just use the w0.5 model, we can see that it has drifted to center on 1udlA, 1gfc, and 1gcqA, rather than the sequence we started with.

Various fixes are possible.

• One is to look at the simple pairwise scores (ignoring the multiple alignment) when choosing close homologs—for this example, the best alignment scores came from the alignments that had identical sequences, which would have worked fine in this case. This reduces to using rather crude tools though, so might not work as well when sequences are somewhat diverged.
• Another possibility is to look for pdb files in the t2k.a2m or t04.a2m alignments (there are 78 of them in this example), and make a phylogenetic tree based on the multiple alignment. The closest templates in the tree are the most likely ones to be worth using. This would work well here, since there are many PDB files in the alignment, and the right ones would cluster no matter what tree-building method was used.
• A slightly more powerful method would be to align all the top hits from PDB to the w0.5 model, then build a tree from that alignment. This may get a few more distant sequences than the t2k alignment building (for example, there are 95 SH3 domains found, not just 78 in the 2awl example). The method could be applied to fold-recognition targets as well as the targets with really close homologs in the templates.
• We could also look at the ordering of the predictions with different weights on the second (or third) track. It seems that the predicted local structure tracks become more valuable for more distant templates. If we find close templates, we may want to use parameters that put more weight on the AA track.

Some possible projects include

• Studying measures of correlated mutation that have been proposed, but which George has not had time to implement. These include fairly simple ones (like mutual information on reduced alphabets) and more complicated ones (like mutual information on "blurred" contingency tables).
• Investigating other definitions of "contact". We are currently using CB-CB<8Å, but are interested in looking at a few other definitions—one involving sidechain proxies that should get at sidechain-sidechain interactions better, one looking for beta-sheet neighbors.
• Cleaner scripts for creating the predictions, visualizing them, and incorporating them into undertaker as constraints.
• Improved neural nets for making predictions (this may be too close to what George is currently working on.)

There was a tool almost usable in 2004: ProtoShop or ProteinShop. Information is available at http://proteinshop.lbl.gov/ http://graphics.cs.ucdavis.edu/~okreylos/ResDev/ProtoShop

ProteinShop seems to have some fairly nice manual-move capabilities, but was not usable last time I tried it (in 2004). I understand that David Bernick did some work installing the program and debugging it, and that the publicly available open-source version may now be usable.

The class project would be to install (and possibly debug) ProteinShop, then try to interface it to undertaker, so that one could look at partially optimized structures generated by undertaker, manipulate them with ProteinShop, and hand them back to undertaker for further optimization. The two programs could probably interact by passing PDB files of the conformation back and forth.

A more advanced GUI would not only allow manual manipulation of the conformation, but would allow applying specific undertaker conformation-change operators in specific regions. Visualization of the "ternary edges" that hold together structures that lack a contiguous backbone would also be useful, though these are not currently represented in the PDB output by undertaker.

One operator that would be nice to have is the ability to insert a particular alignment, or a subdomain from a particular PDB file. There *was* an operator in undertaker for this (ApplyFragmentsToConform), but this command is now obsolete, since it used an old "slicer" format for specific fragments, which is no longer used by undertaker. If we added a NameToPtr hash table to the AlignmentLibrary, we could refer to alignments by name and create a new command "ApplyAlignment".

Save and restore AlignedFragments library

The fragfinder output is as an A2M file, which is fairly compact and easily produced by the SAM software, but which does not include 3D information. Undertaker can read the a2m file and the necessary PDB template files and generate fragments by sidechain replacement, but does not currently have a way of storing these fragments in a file. Existing undertaker file formats could easily be adapted to save the fragments, speeding up multiple undertaker runs and making it easier to run undertaker on the kilocluster (by eliminating the need for access to the PDB library).

Although we can output fragments in Rosetta format, the Rosetta format does not contain enough information to reconstruct undertaker fragments (which include sidechain information, as well as backbone information).

Undertaker cost functions from local structure predictions

What we currently do for alpha is to convert the 11-dimensional neural network prediction vector into a nearly continuous function function of the alpha angle (360 1-degree bins). This involves combining the neural network prediction with the histograms for individual amino acids. The current algorithm starts with the histograms for alpha conditioned on the residue and conditioned on the next residue, and takes their sum to get an initial continuous curve, each of the 11 sections is scaled to get the appropriate probability according to the neural net prediction, and the resulting curve is smoothed to remove the discontinuities at the edges of the sections. The scaling and smoothing is repeated several times.

We could do a similar thing for burial, but the burial curves are simpler than the alpha curve, so it might be better to use a simpler representation, such as just a mean and standard deviation parameter for a family of curves such as the gamma distribution.

In addition to the various burial alphabets, we have some backbone alphabets. These are discrete alphabets and cost functions could be implemented just by having a table lookup of -log(Prob(letter)) for each position—simpler that what we do for alpha.

Currently the "Bystroff" alphabet is understood by undertaker, and it would not be hard to add the de Brevern protein-block alphabet. The new hbond alphabets that we have developed locally are also understood by undertaker, though there is no scoring scheme for them implemented yet. Some of the more conventional alphabets (dssp, stride, and our dssp extension str2) would be more difficult to implement, as we do not want to duplicate the innards of DSSP or STRIDE.

Other ways of using predicted local structure in undertaker

Predictions of the new n_sep and o_sep alphabets could be used to create specific hydrogen bond constraints, to supplement the helix and strand constraints from other local structure alphabets.

The use of Hbond constraints (either from helix and sheet constraints or directly provided) could be improved by changing the TargetA and TargetB functions to provide better estimates of where the donor or acceptor should be. Currently the Hbonds use the same somewhat random target location that other constraints use.

Besides using constraints, local structure information could also be used to affect where fragments get inserted, with higher probabilities for loop regions and for regions where the current secondary structure does not match the predicted secondary structure. With a little more programming effort, we could even make the probabilities of different fragments depend on how well they match predicted structure. For example, we could

• Look for fragment that matches peaks in predicted alpha values.
• Label fragments with DSSP (or other) labels and look for matches to strongly predicted regions.
• When breaks are inserted (say by jiggle or opt subtree), make the more likely to be inserted in regions which have lower probability of being in a helix or strand. Possibly deliberately insert a break in each such region to act as a rapid hinge.

One operator that might be nice to have is an ABA crossover operator. Like the current (AB) crossover operator, it would start from two parent conformations. But instead of having one crossover point, it would have two, using two peptide splices to join a chunk from conformation A, a chunk from conformation B, and the rest from conformation A. This could be done either as two peptide-splice operations (based on how the current crossover operator works), or by extracting a chunk of conformation B as an incomplete AlignedFragments conformation, and applying it to conformation A. This would do a better job of keeping the two parts of A in their current relationship to each other.

Another useful operator (or set of operators) would be ones that jiggle the residues without breaking the backbone. One fairly sophisticated way to do this would be to tweak one torsion angle, then use inverse kinematics to compute torsion angles needed to keep a later peptide plane in position. Somewhat simpler and cruder would be to find pairs of bonds that are almost parallel and to make "crankshaft" moves that tweak the two torsion angles in opposite directions. This is done in a simple way by the TweakPsiPhiSubtree and TweakPsiPhiSegment operators. A variant of this would be to do a "crankshaft" move between residues A and B, add theta to psi(A-1), -theta to phi(A), -theta to psi(B-1), theta to phi(B).

There is plenty of room for thinking up new conformation-change operators!

Improved handling of clashes and bond angles in undertaker

One big disadvantage of the current scheme is that there is no attractive van der Waals term in our cost function. Atoms that don't come into contact score the same as ones that do contact without clashing. While a traditional Lennard-Jones potential gives too high a penalty to clashes for our purposes (making the landscape too rough), the attractive well of the potential is important in determining whether a protein is properly packed. The Lennard-Jones potential is easily parameterized, having only one parameter for each pair of atoms. We need to come up with a similarly simple cost function that is easily calibrated from available data, with a good fit to the observed distance distributions, but with a softer clash penalty.

It is quite likely that bond lengths and H-bonds will not fit the same shape length distributions as van der Waals contacts, so we may want to provide completely separate mechanisms for bonded atoms (we already have a separate H-bond cost function). Atoms that are 2 bonds apart are also not governed by van der Waals contacts; their distances are basically measures of bond angles. The current clash detection scheme has different thresholds for bonded atoms and for atoms 2 or 3 bonds apart (because of the different thresholds for same residue, adjacent residue, and arbitrary chain separation), but does not have different shapes to the cost function.

• Home

Quick Links

• Home
• About Us
• Advising
• Departments
• Research
• News & Events
• Courses
• People
• Administration
• Contact Us

Page Not Found

Home

The web page you requested could not be found.

Our web site is constantly being updated. If you clicked a link or a bookmark to get here, the link or bookmark probably needs to be updated.

Search

You can use our search engine to find the page you are looking for.

Thanks for your patience and understanding!

© 2009 UC Santa Cruz • All Rights Reserved
1156 High St, Santa Cruz, CA 95064 • 831-459-2158 •

UCSC Home • BSOE Home • Web Mail • Forums • SSH • Support • Wireless Registration • Employment • Basic Menus • Log In