Approximately half the grade for this course will be based on the successful completion of a substantial project. Students are expected to turn in

Fri Jan 11: half-page proposal

Students should select a project within the first two weeks of the class, and make a very short proposal outlining what they want to do.

Fri Jan 25: detailed research plan

By week 4, students should have researched previous work on the topic of their project, done some preliminary assessment of the difficulty of the project, determined how much they can get done in a quarter, made a list of milestones, scheduled deadlines for various parts of the project, and perhaps started preliminary coding. The research plan should be roughly the equivalent of the introduction and prior work sections of a research paper, with an outline (with deadlines) for the rest of the paper.

Weekly drafts of final report, each Friday at 4pm

Students are required to turn in a weekly report starting Jan 18. Perhaps the best way to keep the amount of writing to a minimum is for the weekly report to be drafts of the final project report, with change bars in the margin indicating any changes from the previous week's report. The research plan due Jan 25 is one of these progress reports.

Final report, Tuesday, 18 March 2002, 8am.

The final report is primarily what will be evaluated, but the timeliness and usefulness of the weekly reports will also be considered. It is expected that as many as half the projects will be suitable for conference or journal publication (perhaps with an extra quarter's work).

Poster Session, Tuesday 19 March 2002, 8 a.m.--11 a.m.

In addition to the written report, we will have a poster session during the scheduled exam time for the course.

Each student is expected to have a 30-minute weekly meeting with the instructor. These meetings will be used to discuss the project, aid in debugging, and help with homework.

Sufficiently large projects may be tackled by a group of students, but past experience has favored individual projects. Projects can be supervised by more senior graduate students or faculty other than the instructor---these are often the most successful projects. External supervision is in addition to, not in place of the weekly meetings with the instructor.

This is a list of potential projects for students in Bioinformatics I. There are, of course, many other potential projects, and each student will be required to submit a detailed proposal for the project that he or she intends to do.

The best projects are usually ones that students are most motivated to do, so don't consider yourself constrained by this list---if you have a project idea you would like to pursue, propose it!

Project list--Winter 2002

Fold recognition

We need to create a new fold-recognition server to replace the SAM-T99 server. There are many improvements that we have made in the past couple of years that have not been put into a server yet, and many ideas that have not yet been implemented.

In order to make the server easy to extend, we are planning to use a Makefile (using gnu's version of make) to control the computation. This makes it very easy to add new features to web site, separating the flow of control from the web interface. Rachel Karchin will be creating the web Makefile and web scripts, based on one of the Makefiles in /projects/compbio/experiments/protein-predict/ and on the SAM-T99 web site.

Improved homolog search

SAM-T2K needs to be updated to SAM-T02. For example, the last iteration of the search procedure should use a calibrated model. This means adding a calibration step, and adjusting the thresholds for the earlier steps so that the errors in the uncalibrated models don't make the early steps include too many sequences. (One could calibrate at each iteration, but this seems unnecessarily expensive.)

Improved search/multiple alignment by parameter adjustment

Fold-recognition tests should be done to tune some of the parameters---particularly the threshold of the final iteration and the transition regularizers needed. The current fold-recognition test will take about 40-80 CPU days to build all the models and do the test, but there are some details still to work out in doing the scripting. The model building is the slow part---testing a set of models takes only about 36 CPU hours.

A more modest version of parameter tweaking would not change the underlying multiple alignments, but just the weighting of the amino-acid track and the predicted secondary structure track.

The fold-recognition tests I'm currently using are in ~karplus/pcef/fssp-test-01, and are mainly summarized in the Makefile, though there are descriptions of it in some of my recent talks.

Improved search/multiple alignment using psiblast

The prefilter we have for SAM-T2k uses wu-blastp, and it would probably be a good idea to provide a prefilter for NCBI blastp as well. With NCBI blast, we could also speed up the extraction of the FASTA sequence files, which is woefully inefficient in the current implementation. (We would also get rid of the sequence counting step, since blast reports how many sequences it searches.)

It might be interesting to change the pre-filter method to generate psi-blast profiles from the HMMs, and use psi-blast as a prefilter at each iteration. This should allow us to detect slightly more remote homologs than our current pre-filtering strategy, without enormous increase in cost.

One warning: people think of psi-blast as very fast, but only the blast part is fast. When psiblast has to realign many sequences to a profile, it gets quite slow also. It may be that we only need to use the BLAST part of psi-blast, if we generate psiblast profiles from our HMMs. (Our HMM aligning is no faster than psi-blast's but it does seem to produce better alignments.)

Improved multiple alignment using T-coffee

The best multiple alignment program currently is T-coffee (at least based on the rather poor balibase alignment benchmark). We may get some advantage from using T-coffee to produce multiple alignments of the homologous sequences instead of buildmodel. T-coffee has a problem with large numbers of sequences (it is O(n^3)), so we may have to do a quick-and-dirty alignment with an HMM, thin with uniqseq, select a small enough subset of sequences, align the subset with T-coffee, convert back to the alignment format we use, then extend the alignment to all sequences using an HMM.

Improve secondary structure prediction

Our secondary structure predictor is one of the world's best, but David Jones's PsiPred is doing slightly better in the EVA evaluations. The difference is not significant according to the EVA criteria, but it has been fairly consistent---it may be due to their training on DSSP and our training on STRIDE.

I have several ideas for improving the secondary structure predictor, some of which are trivial and some of which are difficult. Some of the difficult ones are already being implemented by grad students (in their spare time), but a collaboration to get the methods finished and tested would useful even on these. Here are a few of the project ideas:

• It would be good to have an evaluation program for secondary structure predictors that is independent of the predictor itself. We have a standard RDB file format for secondary-structure predictions which all our current predictors use, so it would be fairly easy to write an evaluator that reports Q3 (3-state accuracy), the SOV (segment overlap values), and encoding cost for both DSSP and STRIDE assignments. Most of the desired evaluation reports are available from predict-2nd, so one way to create the evaluation program would be to add new commands to predict-2nd that read in external predictions.
• Analyze the failures of the secondary structure predictors. Are there some common features of the sequences, the multiple alignments, or the structures for hard-to-predict sequences? Do the different methods on the web fail on the same sequences or different sequences? Are errors in prediction random or are errors correlated for nearby residues (along the sequence or in 3-space)? How common are serious errors (mis-labeled segments in core structures, broken or joined helices, ...) compared to minor errors at the edges of segments or in regions of the protein not conserved across homologs?

  Our current methods are trained to minimize encoding cost or maximize the Q3 measure (percent of residues with correct classification). The EVA evaluation system uses the segment-based SOV score as its primary criterion. We need to figure out what sorts of errors contribute most to low SOV scores (missed strands, missed helices, broken helices, broken strands, merged helices, merged strands, ...). We can then increase the training weight for particularly important examples (ones that contribute more to the SOV score being bad), trading off small losses of Q3 for (we hope) big gains in SOV. The weighting could be based on the true structure, the predicted structure, or errors in the prediction.

• Test the effect on the neural-network predictor of changes to the multiple alignment algorithm. For example, is it better to include more remotely related sequences (which may contaminate some multiple alignments with noise) or to use a stricter cutoff on similarity? Do the SAM-T2K alignments give better secondary structure prediction than the SAM-T99 ones? In one set of pairwise alignment tests, we found that thinning the multiple alignment to remove highly redundant sequences made the sequence weighting work better. The uniqueseq program can be used to thin multiple alignments before passing them to the neural network. (Currently the best network uses SAM-T2K alignments thinned to 90% identity, but it is not clear that this is anywhere near the optimum.)
• We have a train/cross-train/test set available, so that we can test a variety of different neural net architectures to find which ones generalize the best. We have done no systematic exploration of the space, but have arbitrarily chosen a few architectures with about the same number of parameters. The original idea was to keep the network small, so that it would avoid overtraining, but if we use the cross-training set appropriately, we can avoid overtraining even on large networks. (Note: PsiPred uses a much larger network than we do, so it is probably worth investigating some bigger networks.)

  I have done some testing of different architectures with one part of this data (training on part 1, crosstraining on part 2, final test on part 3). Right now it looks like increasing the number of parameters will cause more overtraining, but we can distribute the parameters more efficiently between the different layers. Perhaps, with a better distribution of the parameters, we could increase the number somewhat.

• Our current SAM-T99 web server combines 3 secondary structure predictors (all working from the same multiple alignment) by averaging the 3-dimensional (EHL) probability vectors of the three methods. Unfortunately, the three methods are not all calibrated the same, so this averaging may not produce the best prediction. Find a better way to combine the methods.

  Simple idea---use a linear combination of the probability vectors and optimize the weights. More complicated idea---create a calibration curve for each method to rescale the probabilities before averaging (or doing linear combination).

  Another idea: take a linear combination of the logs of the probability vectors, and convert back to probabilities by exponentiating and rescaling.

  Another idea: provide a tool for re-calibrating the output of a secondary-structure predictor. That is, produce a monotonic function from the predictor output to the observed probability that the prediction is correct. There may even be a separate function for each of the predictor outputs. The predictor outputs can then be rescaled by this function, and renormalized to sum to 1. Combining these re-calibrated outputs may be more robust than combining the uncalibrated ones. Note: for the neural net, this recalibration does not seem to be necessary---when the output of the network is p, then the fraction of time it is correct is very close to p. The other techniques we are using do not have such a good internal calibration. [I believe that Spencer Tu has done this for his method, but Mark Diekhans method has not been recalibrated.]

• Other predictors often use a jury of multiple neural nets to get a prediction. It might be worth experimenting with the current neural network, training several neural nets (perhaps with different internal structures) on the current data, and seeing if the average of their probability vector outputs is any better than a single network. It would also be interesting to compare with a network containing as many parameters as the combined sizes of the juried networks. Is jurying the best way to use more parameters?

  The jurying could be built into the predict-2nd program by adding a new MultipleNet class that can have combine several networks that have identical input and output interfaces.

  I have built at least a dozen networks for predicting STRIDE EBGHTL alphabet, trained from scratch on the same training data. I have not yet determined whether combining multiple networks would improve the predictions.

  We could investigate several different jurying methods:

  • Averaging the probability vectors.
  • Linear combination of probability vectors.
  • Linear combination of ln P(E | pos)/P(E), exponentiated and rescaled to a probability vector.
  • Neural net with probability vectors from nets as inputs.
  • Neural net with ln P(E | pos)/P(E) vectors as inputs.
• I had some other ideas about ways to combine methods that don't just look at the outputs. For example, both Mark's method and Spencer's method could start from the neural network output and use it in doing the forward-backward computation that produces their output. I suspect, though I can't guarantee, that using one or both of their methods to bootstrap the the neural net output would produce a better result than simply averaging independently created outputs. Mark has done some experiments with using his method to improve the output of the neural net predictor. Initial results are quite promising, but we need to do a proper train and test, and (if it works) we need to install it as part of the standard Makefile.

  We could also create a different HMM-based method that is only useful for this sort of post-processing. It would be trained on secondary structure strings alone using labeled states (no amino acid information and states emit exactly one of E, H, or L). The probability vectors from the neural net could be aligned to the 2ry HMM using forward-backward, and new probability vectors read off by posterior decoding (in each position use the sum of the probabilities of all states with a particular label). This should help eliminate the one-long helices and other ridiculous results that sometimes come up from residue-at-a-time prediction. If the HMM is well deigned, it may even be able to look for common patterns, such as the helix-loop-strand-loop repetition of alpha-beta proteins.

• The neural network predictor currently used has only about 3500 parameters, but we have about 500,000 training examples (2500 protein chains averaging about 200 long). This data/parameter ratio is large enough that we did not have to worry about overtraining the network, and did not hold out any of the data for crosstraining. A possible project is to build larger neural nets using a clean cross-training method, and see whether a larger net could actually get noticeably better performance. (I suspect that the larger neural net used by Psi-pred may be responsible for some of their slight advantage.) The most obvious place to increase the size of the network is the window size for the first layer. [Note: doing proper train and test on a somewhat smaller, less-redundant data set indicates that we may be running into overtraining already with our fairly small network.]

  We could also train on a larger (but more redundant) set of alignments---such as all the SAM-T2K template alignments we've built (currently about 5575 chains).

• Our current methods are trained to predict the STRIDE classification of secondary structure, while EVA evaluates on the DSSP classification. The two differ in about 4.6% of the columns on the EHL reduced alphabet (mainly having to do with the exact ends of helices). They disagree on 16% of the columns when using an EBGHTL alphabet. It might be good to retrain the methods to match both DSSP and STRIDE (perhaps aiming for a probability vector of 0.7 STRIDE, 0.3 DSSP where they disagree).

  To handle both DSSP and STRIDE training data, we need to modify the input routines, since DSSP and STRIDE assign secondary structure for slightly different subsets of the residues. We should probably accept a secondary structure string that has either the same number of characters as the number of columns in the multiple alignment (the current requirement) OR has the same number of characters as the guide sequence. Then the secondary structure string can be aligned to the input multiple alignment by pairing it up with the guide sequence.

  Another approach is to build separate predictors for the STRIDE and DSSP alphabets and report both for the web site. We could then use a 3-track HMM for fold recognition. This is the direction that Rachel Karchin and I are currently pursuing.

• It might better to train the network to (0.98, 0.1, 0.1) (or to (0.69, 0.29, 0.2) if STRIDE and DSSP disagree), rather than (1.0, 0, 0) probability vectors. I hope that this would free up some of the weight-space that is currently used for trying to get very accurate predictions on easy sequences to do slightly better on the hard ones. It would also prevent getting predictions of 0 probability as output, which could mess up the 2-track HMMs.
• In the multi-track HMMs, we use predicted secondary structure to improve alignments and (possibly) fold recognition. For this to work well, we really want to predict the possible secondary structures of homologous and analogous proteins in PDB---a perfect prediction of the true secondary structure may not be of much use if that structure is not identical in the homologs. It would be interesting to put together a training set where the correct secondary structure probabilities are taken from structural alignments (FSSP or SCOP domains), rather than just the target sequence.

  There are some difficulties in deciding how to handle gaps in the multiple alignment. The simplest is to use FSSP representative sequences as the targets, and to use only the secondary structures that are aligned at a particular residue. Another approach would be to have the neural net predict a gap character as well as the normal {E, H, L}, though it is not clear how that gap character would be used in the multi-track HMM.

Predicting Local Structure

Several people have observed that the current secondary structure classification is not a very complete description of the local structure of proteins. There are definitely common conformations of short segments of protein backbones that are not well-expressed in the popular 3-state classification of residues. Some of these are highly predictable (see the ISITES library, for example).

The goal of this project is to come up with other structural features and methods for predicting them. We are looking for a one-dimensional representation here (some feature for each residue) that can be used with HMMs. This project can vary enormously in scope, from quick checks of a single proposed feature set to Ph.D. dissertations. In fact, Rachel Karchin has proposed a substantial investigation of local structure features for her Ph.D. thesis (see /projects/compbio/papers/karchin/thesis-proposal.ps).

Here are some possible criteria for choosing features:

• Features should be well conserved between corresponding residues of homologous proteins. If the features are from a small alphabet, then an information theoretic test can be used---what is the cost of encoding the features of protein A given a structural alignment to protein B? A simple substitution matrix approach can be used to get the probabilities.
• Features should be well-defined, even for somewhat sloppy structures. One problem with the DSSP and STRIDE definitions of secondary structures is that small errors in atom positions can change the classification, particularly for DSSP which relies on finding the H-bonds.
• Features should be at least somewhat predictable from multiple alignments.

Here are some ideas for possibly useful structural features to predict:

• Divide up the phi/psi space of the Ramachandran plot into different regions and predict the region for each residue. Note that some regions correspond to helical and some to strand configurations, but this classification is not equivalent to the standard 3-state secondary structure prediction, since we are not looking for sheets and helices but only the conformation of individual residues. Conjecture: this feature is more predictable and/or more conserved between homologs than the standard EHL feature. We have started some work on this, using Bystroff's clustering of the phi/psi angles (in the HMMSTR paper).
• Predict the numerical phi and psi angles directly, giving both a mean value and standard deviation. (Charlie Strauss indicated that this would be very useful for fragment packing methods like Rosetta.) One possibility for using neural nets is to train the output of an n-output neural net to cosine(delta phi)+cosine(delta psi), for n canonical phi-psi pairs. The predicted phi-psi could be determined from the n outputs by solving a linear system to get the direction vectors. The length of the vector may be useful for determining the standard deviations of the predictions.
• Predict the C-alpha(i-1), C-alpha(i), C-alpha(i+1) dihedral angle.
• Superimpose the backbones based on the peptide bond atoms (C(i-1), N(i), C-alpha(i)) and cluster on the vectors consisting of the locations of O(i-1) and O(i). This is yet another encoding of the phi-psi space, but it also includes the cis/trans configuration of the peptide bond. I think that the location of O(i) may cluster more neatly than phi/psi angles do. The alphabet that results from this clustering may be another interesting one to try. The clustering should probably be done as a Gaussian mixture.
• Predict the sidechain rotamer---at least the first torsion angle. The rotamer hypothesis provides an obvious 3-letter alphabet to map try. It would be interesting to see whether one can predict the rotamer from the sequence---it is suprisingly difficult to predict the rotamer even when the backbone is known. Results should be compared with Dunbrack's SCWRL program, which uses a known backbone.
• Instead of giving the probability of a structural feature directly, give a conditional probability based on some other structural property. For example, predict accessibility conditioned on secondary structure, or phi-psi angles conditioned on a neighboring phi-psi pair. Ideally, we would like to predict a feature that is not found early during fragment packing conditioned on a feature that is generally found fairly early. (For example, fragment packing usually determines what secondary structures it thinks are reasonable fairly early, but the packing of the structures together takes a long time.) The conditional probabilities could be used to guide fragment replacement or scoring in a fragment packer.
• Divide the beta-strand prediction into 9 classes: (left-side: edge 0, parallel strand + , anti-parallel strand -) x (right-side: edge 0, parallel strand +, anti-parallel strand -). This would give us an 11-state prediction instead of a 3-state prediction, and the different beta states may be difficult to predict from multiple alignment information. If some information about the tertiary status of the strand residues could be predicted, it could be enormously useful in scoring for fragment packers. A substantial part of this project would be producing the training set, as I do not know of an existing program that determines the desired classification. (For that matter, can we distinguish "right" and "left" sides, or should we reduce the beta states to 6: 00, 0+, 0-, ++, +-, --?)

Parameter tweaking for improved use of local structure prediction

In CASP4, it looked like a lot of our improvement between the SAM-T99 automatic server and the SAM-T2K hand results came from the new 2-track HMMs that use both amino acid and secondary structure. We need to test different weightings of the two tracks for fold recognition and for alignment. These tests can be done more quickly than other parameter adjustments, since the underlying HMMs do not need to be rebuilt.

Multi-track HMMS for template models

The 2-track HMMs can currently only be used with target models, not template models, since we need an accurate secondary structure for the sequences and probability vectors for the model. If we could create reasonable probability vectors for the templates (perhaps by using context-sensitive substitution matrix on the real secondary structure string), we could use 2-track HMMs with template models to score the predicted secondary structure of the target.

The idea here is to look at what the probability is for the secondary structure predictor predicting a H, E, or L as the most probable letter in a specific context (such as mid-strand, beginning of strand, end of strand, ...). The appropriate probabilities can then be used in building the secondary structure HMM. The creation of the substitution matrix should probably be automated so that (1) the definitions of the contexts is easily changed and (2) new secondary structure predictors can have substitution matrices calculated for them quickly. The substitution matrices would actually be quite useful in designing and debugging secondary structure predictors, so should be a stand-alone tool.

Another way to create 2-track HMMs for template models is to use probability vectors in the sequence and labels in the states (instead of vectors in the states and labels in the sequence). In this way we could use the predictions of secondary structure for the target sequence with the known secondary structure sequence for both the target and the template models. This method would require some (fairly modest) changes to the inner loop of SAM's hmmscore program, and a number of I/O changes.

Calibrating HMMs

We have discovered that some assumptions made about the way statistical significance (the E-values) were computed in SAM are incorrect, and we have created a -calibrate option to hmmscore to calibrate the HMMs properly. Calibration consists of two parts: (1) generating and scoring sequences and (2) fitting the distribution of their scores. The quality of the calibration depends on the sequences chosen and on the quality of the fitting.

I have created a generator for amino-acid sequences that I believe creates better "random" protein-like sequences than standard methods. This generator is now built into hmmscore.

We could still use a generator for paired sequences (amino acid and secondary structure) for calibrating multi-track HMMs. There are several approaches for generating such pairs, and I don't know which will be best:

1. Generate AA sequence as currently, and generate a 2ry sequence using the same approach of choosing a composition distribution (another dirichlet mixture) , and doing independent draws from it. This is what hmmscore does, but the secondary structure strings generated are terrrible decoys. The problem is that secondary structure letters usually come in runs---they are not even approximately independent from one position to the next. As a result the calibration from this process is worse than doing nothing.
2. Generate AA sequence as currently, then use a neural net or other 2ry structure predictor to produce a secondary structure string. This would provide good correlation between the AA string and the 2ry string, but I would expect the 2ry strings to be unusually high in "loop" letters, since the random AA sequence is not likely to have the characteristics of structured proteins.
3. Generate AA sequence as currently, but use a "generalized" HMM to generate the secondary structure string. Each state generates a segment, rather than a single residue. The states could be pure labeled states (just E, just H, or just L) and a length distribution would be kept for each state (either explicitly or as a parameterized distribution). I think that a small HMM of this type could be created that would match the existing data quite well. It might be hard to get overall length right, but probably easier than with residue-at-a-time HMMs. There would be no connection to the AA generator, except through the sequence length---it would probably be easier to generate the secondary structure string from the generalized HMM then generate an AA sequence to match the length. I have segment length histograms and a "segment composition" Dirichlet mixture, so this would be fairly straight-forward to implement for the EBGHTL alphabet.

   We would probably need a set of tools for dealing with these segment-based HMMs, so that we could train them on real 2ry strings and build a generator based on them. There may also be some applications of such HMMs to secondary-structure prediction, post-processing the output of residue-at-a-time predictors.

   A simpler version of this approach would generate the segment labels as almost independent draws (prohibiting duplicating the same letter twice in a row), then use the segment length distributions to generate the whole strings from the segment labels.

   This approach should work-well for 2-track HMMs (since the amino acid track and the secondary-structure track have only about 0.1 bits of mutual information when considered a column at a time), but not for multi-track HMMs with multiple variants of the secondary structure alphabet (since these will be highly correlated).

4. Yet another approach is to generate the AA string and the secondary structure string simultaneously from an HMM, such as the HMM used by Spencer Tu to predict 2ry structure. It should be easy using Spencer's code to write a sequence generator that generates a random path through the HMM with random emissions from each state. The two outputs could be sent to a pair of FASTA files, to match the input expectations of the programs that handle multi-track data. I suspect that this generator would produce the most realistic decoys, but the sequence lengths might come out too short. It should be easy to check that hypothesis---just generate a 1000 or so sequences and make a histogram of their lengths.

In addition to generating the sequences (ideally entirely with C code, so that the sequence generator can be easily incorporated into SAM), we have to fit the parameters of the E-value computation to the observed distribution of scores. I'll cover in class the derivation of the E-value computation we use for the reverse-sequence null, and show how to fit lambda very easily. We can generalize the fitting, introducing an extra variable as an exponent on the scores. Fitting this two-parameter model requires some iteration, but is still fast. The two-parameter model we are fitting has no theoretical justification. it would be nice to come up with a heavy-tailed distribution that fit the data well and had some justification.

Our current best results are from moment-matching to fit the distribution of database sequences, but it would be nice to use a maximum-likelihood (or maximum posterior probability) approach instead. This would probably require a different family of heavy-tailed distributions.

Improve template HMMs by using structural alignments as seeds

An idea that we tried an rejected a couple years ago was to build SAM-Txx models starting with a structural alignment as a seed. The ideas was that structural alignments are better at pulling out the important residue conservation patterns than sequence alignments. This did not work as well as we expected. I think that the problem was that we were using fold-level multiple alignments, so the sequence signal was too diluted to be useful. We might be able to use superfamily-level structural multiple alignments more effectively. We could select the superfamilies either by using SCOP or by using some combination of the FSSP Z-score and %identity.

We have found that HMMs built from structural alignments do a better job of aligning remote homologs (not included in the structural alignment) than HMMs built from SAM-T99 alignments.

Certainly any method that selects between alignments should include alignments build from structurally aligned templates.

Labeling template hits by SCOP family and FSSP representative

It is useful when examining template matches to cluster similar templates together. There are four standard databases that are useful for doing this: SCOP, FSSP, Dali Domain Dictionary, and CATH. It would be useful to annotate each hit in the template library with the appropriate identifiers from the databases.

Note: this annotation is straightforward, but not completely trivial. The databases generally contain DOMAINS, not whole chains. Therefore, to avoid incorrect labeling, one has to examine the alignment between the target and the template and determine which domains of the template are matched by the target. The SCOP database may be the best to use for this purpose, though we will have some templates in the library that are not yet in SCOP.

Improve template selection by dividing into domains

Multi-domain targets often cause problems for fold-recognition techniques. One problem is that a particular domain matches fairly strongly and the other domains are either ignored (no sufficiently strong scores seen to compete with the good domain) or are mistakenly matched to domains that are sometimes paired with the strongly matching domain.

To avoid this problem, we could try to split the target up into separate domains, scoring each one separately. The hard part comes in figuring out where to put domain boundaries. Some possible methods include

• Using a domain template library to score the target sequence and split off parts that score well. The scoring should not be local but domain-global (matching all the template to part of the target).
• Using a non-structural HMM library (such as Pfam) to recognize domains.
• Chopping up the target based on the alignments obtained by the whole-chain models (template and target).
• Chopping rather arbitrarily based on the distribution of domain lengths.

One could try several different (possibly even overlapping) domain divisions, then select the one that gets the best matches to domain templates.

Improve template selection by combining information from multiple templates

Timothy Bailey and William Grundy have created a method for combining correlated p-values to get improved fold recognition. They applied their method to improve a method based on pairwise alignment (a rather weak initial method). It would be interesting to see if their method would improve fold recognition when applied to a more powerful method (such as SAM-T2K or SAM-T02). I suspect that it would improve the performance, since the extra information that several models for the same superfamily score well has been used in our human-intervention methods, but not so far in the automatic ones. A crude test of this was done last year in 243, and it seemed to help, but the test needs to be redone with better controls, and the full template library calibrated. The software for maintaining the template library as new templates are added also needs to be developed.

One problem with the product-of-p-values method is that all the templates in a fold family are treated equivalently. In practice there may be several templates from one subfamily and only one from another. Ideally, a scheme should be developed that can weight each template appropriately.

Improve template selection with functional information

In doing structure prediction for CASP, we often look at the functional description of the protein to choose between possible templates. It would be good to build some automatic methods for this. One promising preliminary approach is the SAWTED tool, which is part of the 3D-PSSM fold-recognition server.

The algorithm is quite simple: keywords and text are extracted from the best-matching Swissprot entries of the target and templates, and dot-products of their frequencies of words (excluding non-discriminating words) are used as a similarity measure. Actually, one usually normalizes to a "cosine" similarity measure:
a.b / (sqrt(a.a) sqrt(b.b))

The project would be either to re-implement and test the SAWTED algorithm here or to develop a similarly simple algorithm and test it. A possible variation: use the SAWTED score as a kernel function for an SVM multi-class classifier. The fold prediction method should be evaluated by itself and as a complement to an existing method (such as SAM-T99). (Possible Ph.D. project: developing a more sophisticated way to handle functional information.)

Some preliminary results were obtained last year and this method looks promising. A possible starting point for a scoring function is to take the E-value from the 2-track t2k-1-.3-ebghtl fold recognition results and divide by the "cosine" similarity measure raised to some power.

We need to look into more sophisticated representations than just a "bag of words". Possible improvements include stemming and multi-word phrases. The version last year only extracted information from the single best-matching Swissprot entry, when we probably should look at the top several Swissprot entries.

Improved template selection and alignment with HMM-HMM alignment

We have had moderate success with HMM-based structure prediction, but we could probably do better with it. One thing that we have not experimented much with is to use HMM-HMM alignment. It might be good to try doing local alignment with matching state i of one HMM to state j of the other costing -log sum_a P(a_i)P(a_j) (the co-emission probability). Figuring out what to charge for insertions and deletions is more difficult, as the different HMMs have different transition costs (we could simply ignore the costs from one of the HMMs, for example). Local alignments could be tested either for alignment quality (using Melissa Cline's alignment tests) or for fold-recognition capability (using the test set that we have used for testing SAM-T2K).

Another approach to aligning two sequences is to take the component-wise product of the posterior decoding matrices for sequence A against HMM B and sequence B against HMM A, and doing posterior decoding of the resulting matrix.

Improved target/template alignment

We have not yet determined the best way to to align target with template. There is a pairwise alignment test set documented in ~cline/for_kk/aligntest.html and ~cline/for_kk/comparison.html. These tests have not been applied to multi-track HMMs yet, but we are VERY interested in knowing whether the multi-track HMMs lead to better alignments (and not just better fold recognition).

Selecting good target/template alignments

So far, the only automatic ways we have of selecting alignments are to look at the E-values reported by the HMMs that created the alignments or to select a general method (such as aligning the target to an HMM based on a structural alignment of templates) based on some alignment tests. In the hand-selection of alignments we examine the 3d structure that the alignment predicts for the target.

In particular, we look for features like the compactness of the resulting model, whether internal strands of beta sheets are missing (though surrounding strands are present), whether active-site residues are preserved, whether identical residues cluster in 3-space, whether they form stripes across beta sheets, and so forth. Some of these checks can be automated (such as the compactness test). It would be good to automate some of the checks, and test them to see if they improve the fold-recognition capabilities of the SAM-T2K method. The first checks I would attempt to automate are radius of gyration and contact order, both of which are easy to compute from a structure. We'd have to get histograms of these measures for real structures and for "good" predictions, to see what to expect.

Other 3D scoring techniques (such as threading score functions or fragment-packing score functions) could also be used to filter the top hits from the HMM scoring.

One would have to be careful with any scoring function, since the sizes of the alignments can be quite different and that can affect the scores more than the quality of the prediction does.

Viewing predictions

For viewing our predictions, we have been using rasmol as the graphics engine, and a locally produced program (saen) that shows the alignment and allows selecting residues from the alignment as arguments for rasmol commands. It also allows editing the alignment, highlighting the identical residues. This program really needs to be rewritten to be more easily ported between systems (it currently relies on the AceDB package for doing graphics, but should probalby use something like OpenGL instead). It is working on the old Alpha systems, but has not been ported to the Linux boxes yet---it would be very useful to have it rewritten and improved.

Undertaker

If several people pick projects from this list, we will probably have to start using cvs to manage the undertaker source code, or we'll never get all the pieces integrated.

• Add more consistency checks to the code, to make debugging easier. For example, the ConformationPool class should probably check that a new Conformation being added is compatible with the existing Conformations (either the same Master, or one with the same atoms in the same order). This sort of code cleanup, though important, is probably too low level for a class project.
• The PDB reader included in undertaker (in the ChainsResiduesAndAtoms class) is very crude. It doesn't even understand the notion of multiple models in a PDB file. It would be good to improve it to understand (and select from) multiple models, to include temperature information in the ChainsResiduesAndAtoms class, to include the unit cell information, and perhaps to understand MODRES records also.

  Along the same lines, undertaker currently cannot do much with one of the most common modified residues in PDB files, pyroglutamic acid (PCA) which is usually modified from glutamic acid (GLU). It might be interesting to examine what would would be required to handle it as well. (Of course, we rarely have the information that a particular residue has been modified until after the structure is solved, so I'm not sure when undertaker will be able to use the information in a real prediction.)

• Two years ago, Eric Savage wrote a very useful program (slicer) that converts an alignment of a target to a PDB template into a ChainsResiduesAndAtoms file. Unfortunately, he used a somewhat different c++ programming style than the rest of undertaker, and I have not attempted to integrate his code into the main body of undertaker. It would be nice to have everything "under one roof", so that alignments can be read in by undertaker directly.

  This integration is a relatively small project, but other improvements could be made at the same time, such as cleaning up the way that temporary files are created and using the improved PDB reading suggested above. A bigger change would be to provide side-chain replacement without having to call SCWRL. It turns out that one can do almost as well for an initial guess just by using the most probable rotamer from Dunbrack's backbone-dependent library (except for the residues that are identical, where you just copy the residue from the template). It might be good to use Dunbrack's library directly. In any case, the version of SCWRL we are using is a bit out of date, and should probably be replaced by a newer one from Dunbrack.

• In summer 2000 Rachel Karchin wrote a program to look for short fragments using SAM posterior decoding. It would be good to integrate this program with slicer, so that the fragments are available ready to splice in place. The program might also need to be a little more tightly coupled to the SAM HMM code, so that only the fragments are output, rather than the whole posterior decoding matrix. I'm not sure whether Rachel's program uses 2-track HMMs, but I suspect that we could get much better fragment libraries if we did use them.

  Also, right now undertaker does not really distinguish between the short fragments found by Rachel's program and the big alignments slicer is really intended for. It might be good to keep the two sorts of chunks in separate data structures in undertaker (probably still using the AlignedFragments class, but with two different variables in the Globals class, and with different operations in the genetic algorithm that uses the alignments or fragments).

• Some of the "generic" fragments that are taken from a large training set of PDB files have bad bond lengths or angles. It would be good to add some more checks to ChainsResiduesAndAtoms::bad_residues() for residues with bond lengths more than 2 standard deviations from normal or for residues with high temperatures (see the PDB reading project above). It may also be desirable to have an external file format for marking bad residues, so that programs with more sophisticated notions of badness could be used to filter the residues. For example, Christian Barrett has been using a program that fills in the hydrogens and looks for clashes with the hydrogens and not just the heavy atoms. This eliminates some incorrect rotamers, and could be a good prefilter for the generic fragment library.

  If the residues are marked, then the FragmentLibrary will not include those residues in generic fragments. The "rotamer" library would also exclude the bad residues.

• The "burial" function is currently computed on the kestrel processor, but uses an API (application programming interface) that has not been rewritten for the Linux-based drivers. Only one board is installed with the old Windows driver, making something of a bottleneck. There are two solutions:
  1. Rewrite the API for the new Linux drivers (which should be fairly easy, since the drivers were written with the API in mind). Some improvements could be made in the API, and undertaker modified to use them. There should be a pretty big gain in speed from using the new drivers, as the old one was quite inefficient. [I think there is a new Linux API, so this may only involve slight rewriting of the code for the difference between the APIs.]
  2. Write a serial version of the burial function using a better data structure. The computation that kestrel is currently doing is to count the number of atoms within a specified radius of each of a large number of points. On the kestrel, this is done most efficiently as an all-against-all computation, since the time is mainly determined by how fast we can feed the atom locations through the array.

     In software, we might be better off binning the atom locations and only checking atoms in those bins that are sufficiently close to the center of the sphere we are interested in. I serialized the clash computation (Clash.cc) in this way and got much better performance with the serial version than the parallel one (of course, the parallel algorithm there was a much slower one than the burial computation).

     A serial version of the burial function is probably desirable even if it is slower than the parallel one, as it is easier to distribute software than hardware. We could run a serial version on the new cluster, but we don't have 1000 kestrels.

• Undertaker does not currently use secondary structure prediction, though we have one of the world's best secondary structure predictors (and ideas for improving it---see above). There are several three ways I can see introducing this information:
  1. Use 2-track HMMs for selecting probable fragments (see the section on Rachel's fragment program above).
  2. Compute the secondary structure of the conformations and use the predicted probabilities as part of the score function. This might be better if we used a different definition of local structure that could be determined more quickly and locally (see Predicting Local Structure above).
  3. Keep track of the secondary structure for the fragments (both generic and the specific ones found for a target) and score conformations based on the probabilities of these. This is less expensive than recomputing the secondary structure, but may introduce errors, since non-local conformation changes can change the "secondary" structure (say by moving a strand away from the sheet it is part of).
  4. Add a new way of generating fragments that uses a "canonical" backbone for predicted helices or strands. I suspect that this might work ok for helices (which are fairly superimposable), but not so well for strands (which tend to twist in different ways).
• The 3d filtering proposed for the fold-recognition problem, using radius of gyration and contact order (and possibly other easily measured properties) could also be used for filtering predicted backbone conformations for undertaker.
• Undertaker could, in theory, handle quaternary structure (that is, multiple chains) as well as it does tertiary structure. The main task would be to generate conformations for multimers, and make sure that the scoring functions handled the conformations correctly. The simplest case is probably homo-multimers in which each molecule can be obtained by repeated rotation and translation. This covers a very large class of homo-dimers, homo-trimers, ... . Only one molecule's conformation needs to be explicitly stored and modified. The only other information needed is the transformation to get from one molecule to the next. For an n-mer, the transformation must be an nth root of one, which puts strong constraints on the quaternion that specifies the rotation.

  The project here is to add the necessary data structures and operations to handle quaternary structure, and to come up with some heuristics for generating plausible transformations. One simple heuristic would be to determine a transformation from a plausible multimeric template. Another might be to conjecture that certain exposed hydrophobics need to touch.

• It might be good to recompute the burial distributions used for scoring taking into account all the atoms of the crystal structure in the training set. Currently, only the atoms of the single chain in the training set are counted, which undoubtedly provides too many low-burial hydrophobics. Not only would it be good to include the atoms of the other chains in the asymetric unit cell, but it might be good to include the closer atoms from neighboring unit cells (the "crystal contacts") (see the PDB reader project above).
• The current scoring function in undertaker has a "side-chain cost" component. The side-chain cost is not modified by most of the operators in the genetic algorithm, but is recomputed anyway. It might be better to cache the costs (both in the conformations and in the fragments) and only recompute when necessary---the residues immediately adjacent to the peptide bond where a fragment is spliced in, or when the rotamer is changed.
• The "rotamers" in our "rotamer" library are not really rotamers, but the residue conformations that come closest to the centers of clusters found by a simple clustering algorithm (they are called "archetypes" in the undertaker source code). It would probably be a good idea to examine this library closely and compare it with real rotamer libraries. Are we missing any rotamers? Have we assigned too high a probability to unusual conformations? To aid in this analysis, it would be a good idea to be able to output the archetypes in PDB format, check the bond lengths and angles, convert to the torsion-angle coordinates that rotamer libraries are usually expressed in, and so forth.

  We could also read in Dunbrack's Backbone-Dependent Rotamer Library and evaluate our score function on each rotamer, and compare our assessment of probability with his. Where there are big discrepancies, we'll probably have to improve the way we handle archetypes, either by using a different representation or by clustering differently. The quality of the rotamer library could depend also on how well "bad" residues are rejected (see the project on rejecting bad residues above).

  Undertaker has the ability to call SCWRL to set the rotamers. I have done some informal comparisons of SCWRL and our "rotamer" selection method, and found that SCWRL does somewhat better with the correct backbone. I have not examined which does better when the backbone is a little bit off (a more realistic test). We will need to install the latest version of SCWRL for this to be a fair comaprison.

• The ReduceBreak operation in the genetic algorithm attempts to improve the score by reducing the size of a gap in the backbone. Right now it has a crude heuristic that tries to estimate how much a fragment will help in closing the gap, but it does not compute the new gap correctly if the fragment actually crosses the gap and moves it to a different place in the chain. Improving this heuristic might be valuable.
• One way that homology modelers have used for closing gaps is to have a library of loop fragments. The library is indexed by the length of the loop, the takeoff angles, and (perhaps) by the chain positions of prolines and glycines in the loop. It might be useful to add a loop library to undertaker. The library, like the generic fragment library, is probably best created on the fly from a set of training chains, rather than having yet another file format.

  Using a loop library would probably require having something that does sidechain replacement on the backbone fragment (see the project above about integrating slicer into undertaker).

  Note: Baker's group did very well at loop modeling in CASP4 by using fragment packing without a loop library. They did however, have the ability to freeze parts of the backbone and use the fragment packing only on the loop regions. This should not be too hard to add to undertaker.

• The current data structure for conformations (actually the AlignedFragments class) can introduce new chain breaks, and we have operations in the genetic algorithm for reducing the size of chain breaks, but we have no operator currently for recognizing that a chain break has been essentially closed and removing it from the data structure.
• The scoring function used for undertaker has bonuses for bring cysteines together (on the theory that they might be part of a disulphide bond). One weird effect is that it sometimes brings three cysteines together, rather than pairing them. This could be fixed by using a maximum matching algorithm on the cysteines, and only giving bonuses for the pairs that match.

  The max-matching algorithm may have some undesirable side-effects. For example, CYS residues often cluster because they are coordinating a metal ion, rather than because they are in disulphide bridges. This might be fixed by adding another bonus value for cys clusters that are more widely spaced than disulphide bridges (and excluding any known or predicted bridges from being in such clusters).

• One of the slowest parts of undertaker right now is in the generation of new conformations once an alignment with chain breaks has been inserted. In order to maintain the connected tree of fragments, we make a minimum spanning tree using the distances between atoms to get the distances between segments. The recomputation of so many atom-atom distances is unnecessary and slow. There are several approaches to speeding this up, including using only c-alpha distances, binning the atoms and only computing distances for atoms between segments that may be closer than the minimum distance found so far, finding a different way to maintain the connectivity of the tree, ... .
• The score function in undertaker currently has two "burial" terms: one for "dry spots" and one for "wet spots". I have convinced myself that I need to generalize this to any number of terms, and initially try three terms: small wet spots, small dry spots, and large dry spots. The current data structures assume exactly two spots in several places, so generalizing to n spots will touch a lot of the code, but be pretty straight-forward.
• The current score function has several terms and we will undoubtedly add others. I have no idea how best to combine these terms. At the moment, I'm using a linear combination of terms, and I play with the weights on each different target, to see what happens. It might be good to have a more principled way to set the weights. For example, we could look at some badness function of a conformation (say, sum of square Angstrom deviations plus number of clashes) and do a least-squares fit of the score function to the badness for a large decoy set.

  It would be particularly good if we could do the fitting across several different structures, to avoid target-specific biases. To do this well, we would need a badness function that has similar scaling in different targets. Since the score function is a sum of terms (with terms for each residue), it will tend to grow linearly with the chain length. That means we probably want the badness to grow linearly with chain length also.

  For the least-squares fitting, we probably want to include a QR decomposition routine into undertaker. In the past we used lapack++, clapack, and blas, but we might want to find a smaller, more easily incorporated way to do the QR decomposition. For a very simple proof-of-concept, a preliminary version could be done using an external least-squares fitter, but this might make some things more difficult. For example, we might want to fit the weights of the score terms, then generate a bunch more decoys by trying to optimize with the "improved" score function. Certain terms may not have mattered much in the original set of decoys, because the phenomena they were supposed to detect did not occur, or were masked by other problems.

• Another least-squares fitting project would be to try to improve the burial function by weighting different atom types differently. Instead of just counting atoms, we'd add the weights of the atoms in the spheres. One way to try to set the weights is to compute the mode or mean of the burial distribution for each residue type, then do a least-squares fit of the weighted burials to the "typical" burials for the residues. If the weights come out different from all ones, then we have found a set of atom weights that give more peaked distributions of the burial, and may provide more information.
• One of the genetic algorithm operations in OptCommands.cc (ReduceBreak) tries to select a fragment that makes some progress toward its goal (reducing one gap). It might be good to introduce other operators that are somewhat goal-directed, or that quickly reject obviously bad things. For example, the InsertFragment operator could check for bad clashes and reject fragments that would introduce too many. There is a delicate balance here between speeding things up by not having to do full evaluation of conformations, and cutting off major parts of the conformation space (possibly making the true conformation unreachable).
• There are other ways of creating conformations besides the genetic algorithm that I have been playing with in undertaker. The very successful Rosetta program uses simulated annealing with just fragment insertion as the operation. Other decoy-creation programs use a backtracking search to add one fragment at a time to the end of a chain, creating only clash-free conformations. The traditional physics-based methods use molecular dynamics or energy minimization of atoms in 3-space to create conformations.

  It might be interesting to create one of these other conformation generators, and see how it performs on the same score functions as are used in undertaker. Doing so may require that scoring and conformation generation be a bit more decoupled than they currently are in undertaker, but the code is fairly close to decoupling the generation from the scoring.

Protein projects that are NOT part of our CASP effort

Creation of protein web pages

I would like to create a script for quickly building a web page of useful information and links for a specific protein. You can get an idea of what sort of web page I'm thinking of by looking at the mockup for DNA-repair proteins. When the script is completed, we should offer to make up web pages for each of the proteins being studied by biologists or chemists here at UCSC.

A variant of this project for the non-computer people: Pick a protein or protein family whose structure has not been solved yet, and use existing tools (here and on the web) to do the best job you can of predicting the structure, homologs, evolutionary relationships, active site residues, function, ... of the protein. Prepare the result either as a proposal for specific wet-lab work to verify computational hypotheses or as comprehensive web page for a lab that is studying the protein. This project will require a lot of library work and use of web tools (for example, the metaserver at bioinfo.pl), but not much original programming. This is an end-user's view of bioinformatics, rather than a tool-developer's view.

Last year, Lydia Gregoret expressed an interest in creating a web site for "all things cold shock protein related". She said

There might be links to

• labs studying cold shock protein function
• labs studying cold shock protein folding
• labs studying cold shock protein structure
• evolution of the cold shock fold
• automated PubMed search for cold shock articles
• automated PDB search for cold shock protein structure
• nice images of cold shock proteins

among other things. An example site on lectins can be found at http://www.cermav.cnrs.fr/databank/lectine/ This is not exactly what I had in mind, but a start.

The student would learn about the relationship between protein sequence, structure, and function, what kinds of information is out there and also gain experience in scientific web page design.

Lydia is no longer at UCSC, but other researchers in biology or chemistry may want web pages for their favorite protein families---ask around!

Support Vector Machines for subfamily classification

The hidden Markov models (HMMs) used by SAM have been our main tool for protein family recognition. They work pretty well at the superfamily level, but are often not discriminating enough to distinguish subfamilies.

Rachel Karchin has used support vector machines (SVMs) to build subfamily recognizers for the subfamilies of the GPCRs (G-protein-coupled receptors). She started with the methods and code developed by Tommi Jaakkola and Mark Diekhans, but has since rewritten the SVM programs in Java for easier modification and portability.

Projects here include building SVMs for classifying other families of proteins (using existing tools), trying to improve the SVMs by reducing the large feature vectors, experimenting with different ways of generating feature vectors, trying other ways to use the feature vectors besides SVMs, ... .

Classifying proteins by SVMs is an attractive problem, because we usually have a fairly small set of positive training examples. The big problem is that we often have an extremely large number of negative examples. We might be able to make faster, more accurate SVMs by doing pre-screening with the HMM used to create the feature vectors. If we pick a threshold for the HMM score that is loose enough to include all the positive training examples and about 100 false positives, and use the SVM only for classifying the sequences that pass this HMM test, then the SVM only needs to learn the difficult cases that really are near the boundary. This should result in smaller classifiers (fewer support vectors) that do as good a job or better and are much faster to train [training time is generally quadratic in the number of training examples].

One particularly useful project would be to build a library of classifiers (or a multi-family classifier) for all the SCOP superfamilies. This could be used as a fold recognizer, followed by HMM-based aligners to produce the alignments. If we could get this automated, I think it would make an excellent new structure prediction server. In version 50 of SCOP there are 934 superfamilies, ranging in size from one structure (160 such superfamilies) to 1822 (for superfamily 2.1.1--immunoglobulins). Other heavily populated superfamilies include lysozyme-like (4.2.1, 648 members), globin-like (1.1.1, 612 members), trypsin-line serine proteases (2.44.1, 510 members), NAD(P)-binding Rossmann-fold domains (3.2.1, 454 members), Thiamin diphosphate-binding fold (THDP-binding) (3.31.1, 421 members). [Note: a SCOP classifier would be of use in the fold-recognition experiments for CASP.]

New kernels for Support Vector Machines

The "Fisher score" vector introduced by Tommi Jaakola seems to be a pretty good feature vector when used with a properly scaled radial-basis kernel. It is a rather large vector, though, and many of the features are probably just adding noise. It would be interesting to try some experiments in selecting features out of the vector or weighting each feature differently. Also, the radial-basis kernel only looks at the magnitude of differences between features, not at the sign of the difference. It might be interesting to explore other kernels that could use the sign information as well.

This project would probably consist of a detailed exploration of the feature vectors for one, difficult classification problem (perhaps one of the GPCR subfamilies that Rachel Karchin has worked on). We would want to look for ways to measure the value of a particular feature, then either weight features or select a subset of them, train an SVM using the modified feature set, and compare its performance with an SVM trained on the full feature vector.

Here are some preliminary thoughts on ways to look for informative features:

• Compute the variances for differences in the feature value between positive examples, between negative examples, and between positives and negatives. For a good feature, we would expect a low variance between positives, and a high variance between a positive and a negative. A great feature would also have a low variance between negatives. Perhaps we could weight features by var(pos to neg)/var(pos to pos), or var(pos to neg)/sqrt(var(pos to pos)var(neg to neg)).
• Look at the distribution of the signs of the feature differences between positive and negative examples. One sort of good feature would have mostly positive or mostly negative signs. Look at the entropy of the distribution of the signs P(>=0) log P(>=0) + P(<0) log P(<0). A low entropy indicates a valuable feature. Note that this analysis is more appropriate for a kernal that uses the signs of the feature difference, unlike the radial basis, which only uses the magnitude of the difference.
• Look at the distribution of the feature differences. Try multiplying each feature by the mean of the pos-to-neg differences for that feature (and perhaps divide by the standard deviation of the pos-to-pos differences).

Another variant on the radial-basis kernel that I want to explore (perhaps with Rachel Karchin) is replacing the "sigma^2" scaling factor with a separate scaling factor for each vector in the training set: K(i,j) = exp(-0.5 ||fi-fj||^2 / (sigma_i sigma_j)). The scaling factor could be set to the distance from the vector to the closest vector with a different label. For unlabeled (test) data, the scaling factor could be set to the maximum over the possible classes of the minimum distance from the test vector to the training vectors in the class. Having different scaling for different training vectors may help generalize the classification better, since the kernel won't drop to zero so quickly away from the data, except in the neighborhood of the boundary.

TeXshade

There is a package, TeXshade, for displaying multiple alignments in LaTeX documents. This package does not currently accept the A2M format that we use for multiple alignments. Adding that capability to TeXshade would be quite useful. We might also want to add to our existing web servers a TeXshade output option for alignments.

Compare different protein family diversity measures

There are several ways to measure the diversity of a family of proteins. For example, we can use the total sequence weight needed to get 0.5 bits/column saved using the recode3.20comp regularizer on a multiple alignment of the family. Another method is to use the ESIZE algorithm described in this paper, which relies on fitting the distribution of scores from a number of random queries.

One could investigate the relationship between effective family size on the one hand and percent sequence identity, pairwise sequence p-values, and other measures of sequence similarity. Characterize existing databases (Pfam, SCOP, etc.) in terms of the sizes (using different measures) of the families represented. The ESIZE algorithm is available as part of the FPS software, available here.

Predicting protein-protein and protein-nucleic-acid interaction sites

It may be useful, once a structure has been determined, to predict what locations on the protein might interact with other proteins. One approach is to predict for each surface residue whether it will interact with other proteins. This could be based on the identity of the residue, its chain neighbors, its neighbors in 3-space, or local structural features of protein. Huan-Xiang Zhou claims a 70% success rate on interface/non-interface predictions (not published? at UCSC Bioinformatics Symposium Dec 2000). The biggest problem with this project is getting training data. Zhou had to use crystal contacts to get enough training data, and it is not clear to me that crystal contacts are representative of biologically interesting protein-protein contacts.

Postdoc Yael Mandel-Gutfreund is particularly interested in predicting protein-nucleic-acid interations and may be willing to help supervise projects in this area.

Open source projects

There are many open-source projects in bioinformatics: bioperl, biojava, biopython, bioxml, biocorba, ensembl, GAME XML, DAS, E-CELL, FAKtory, Telegraph, ... , but UCSC has not been part of these. You can find a partial list under the BOSC 2000 conference site. If you want to contribute something useful to one of these projects--pick the project, identify the need you plan to fill, and write a proposal for the class.

There are some specific open-source projects that would be interesting to me:

Open source library for Dirichlet mixture functions

We have done a lot of pioneering work with the use of mixtures of Dirichlet distributions, but the ideas have been slow to catch on, because many bioinformatics people are scared of the mathematics. An open-source library of useful routines would go a long way toward easing the incorporation of Dirichlet mixtures in programs. We have implemented Dirichlet mixtures in C (for SAM) and C++ (for the "ultimate" library), but neither of these implementations is really modular enough to easily incorporated in other people's code.

I have started an open-source C library (with a random generator for Dirichlet mixtures), and would like to see it expanded into a more complete Dirichlet mixture library, while maintaining a very modular, minimalist approach to the code.

There have been requests for source code for our sequence-weighting techniques (which rely on Dirichlet mixture computations). It would be very useful to extract these sequence weighting techniques from predict-2nd, estimate-dist, or the SAM software suite, and incorporate them into a open-source project. The biggest difference between our different implementations is the different implementations of the multiple-alignment data type. If you started with an open-source project that had an adequate representation, it could be quite straight-forward to port the code.

Regularizer library for BioJava

We have a moderately large library of useful C++ classes in /projects/compbio/ultimate. It might be useful to redesign some of these and release them as part of BioJava. I think that the Alphabet and Regularizer types are particularly useful, though there may need to be some redesign of Alphabet to be compatible with existing BioJava design.

Add to Telegraph open-source dynamic programming project

There is an effort to make an open-source general-purpose dynamic programming package (called Telegraph). Ewan Birney, the author of Dynamite and Genewise, has described Telegraph as "Dynamite done right". Unfortunately, Telegraph is currently far from completed, and exists mainly as PERL code, which is too slow to be of any real use. Creating PERL-callable C modules for Telegraph could make it into an important and useful tool. This is an excellent way to improve your understanding of dynamic programming algorithms.

Yet another structure-structure aligner

It seems that every group wants its own structure-structure aligner---no one is ever happy with the existing ones. The "standard" method now seems to be to extract the secondary structure elements and align them with double dynamic programming, then refine the alignment to minimize RMSD or some function of the difference in the distance maps. Some methods are optimized for making very rapid similarity decisions, others are optimized for high quality alignments. Some use rigid superposition, others try to take into account small amounts of flexibility (say by using distance maps), though almost all structure-structure aligners have trouble handling hinged regions (whether domains or loops). Some methods prefer to align large numbers of residues (even if the correspondence is loose), others align only the most similar parts of the structure (even if this leaves many gaps).

Dali,Vast, CE, and Kenobi are worth looking at for inspiration.

Few of the structure-structure aligners use sequence information to help disambiguate the structural alignment (one exception is Jung and Lee's program Sheba).
Use of residue pairs in protein sequence-sequence and sequence-structure alignments
Protein Science (2000), 9: 1576-1588.
Jongsun Jung And Byungkook Lee
It might be interesting to incorporate information from a posterior decoding matrix (which is generally better than a sequence-sequence alignment) into a Dali-like structure-structure comparison.

For an example of where ignoring sequence can result in poor structural alignments, consider DALI's alignment of 1c1gA, 1c1gB, 1c1gC, and 1c1gD. These are all identical sequences, but 1c1gA and 1c1gB are aligned with only 11% residue identity in FSSP, because the coiled-coils are slightly kinked in the solved structure.

Multiple alignment of protein structures

Most of the multiple alignments of protein structures have been created as "master-slave" alignments. It might be interesting to create multiple alignments of superfamilies or folds by using T-coffee to combine pairwise structural alignments. T-coffee has the ability to accept any number of different alignments of the same pair, so different structural aligners could be used to build the library of pairwise alignments. One problem with this approach is that T-coffee does not seem to have a notion of unaligned residues, which is clearly necessary when aligning structures. It may be necessary to re-implement T-coffee with a better underlying model of alignments (allowing unaligned regions) in order to use the approach to get good multiple structure alignments. Another approach is to use a genetic algorithm to adjust the multiple alignment, as is done for pairwise alignments by Kenobi.

Note: the Dali Domain Dictionary uses T-coffee with DALI structural alignments to create multiple alignments, so there may not be much new that needs to be done here.

Protein legos

The success of the Rosetta fragment packer has led to a lot of interest in finding small, reusable pieces of proteins. Some people have concentrated on pieces with continuous backbones (which I call "fragments"), others have looked for conserved chunks (such as helix-helix packing and helix-sheet packing). Zhipeng Wang (Boston University) has called these repeated 3D chunks "protein legos", and has attempted to find them automatically using the Kenobi structure-structure aligner. The idea is to chop a pair of protein structures up into pieces and form alignments of the pieces, then expand to alignments of parts of the protein. The resulting protein chunks can then be put in a library and all-against-all alignment done to cluster them. Highly populated clusters are common "legos", and their frequency can (perhaps) be used in scoring functions for fragment packers or to guide fragment packers toward popular ways to pack helices and sheets.

The project here could be to build a small "legos" library, or to try to find a use for Zhipeng Weng's library within undertaker, Rosetta, or some other fragment packer.

Build a threader

We have worked with 1D (HMM) methods for protein structure prediction, and have started working on fragment packing, but have never built a traditional threader, which tries to align sequences to fixed 3D frameworks. It would be useful to have a threading program for evaluating score functions and to round out our suite of protein-structure prediction tools. This project would be a good way to learn bout things like the frozen approximation and double dynamic programming (see, for example, David Jones's Threader programs).

The scoring functions for fragment packers and threaders have different constraints, as fragment packers have reason to believe that the local structures of the backbone are reasonable, but need to determine of the overall conformation is protein-like, while threaders know that all conformations are protein-like, but that the residues may be in highly unlikely conformations (locally or in terms of tertiary contacts).

Because the objectives of threading and fragment-packing score functions are different, it may be useful to evaluate score functions in both contexts (before throwing them out as useless).

Hubbard plots

At CASP3 in 1998, Tim Hubbard introduced a new way of looking at the quality of a prediction of protein structure, now called "Hubbard Plots". The plots for CASP3 were very useful in determining whether a particular prediction was good relative to others submitted. Similar plots, using a slightly different method, were created for CASP4, but the results are not public yet.

I would like for us to be able to generate such plots efficiently, using the c++ library of routines that we have developed for handling protein structures, transformations, and so forth.

Note: this project was done last year, but the results were not installed. It may be a very small project to clean up the code and install it, but there are certainly many more features that could be added to the program.

Web interface for SAM-T99 library

We have a large library (over 5700 entries) of multiple alignments and HMMs built using SAM-T99 and SAM-T2K. This library is currently not available to the public, though it is used by the SAM-T99 web page for searches. You can browse the library (with no useful assistance) currently at http://www.soe.ucsc.edu/~karplus/sam-alignment-library/

The project is to provide a web interface that allows retrieving individual multiple alignments and HMMs, by giving PDB sequence ids, keywords, or protein sequences. When we do not have exactly the requested sequence, similar sequences should be offered (perhaps using FSSP's TABLE2 to suggest sequences). The search by sequence should probably use some version of blastp, though an option for doing the more expensive search by scoring the sequence against all models should be provided.

It would be good to offer the alignments in the same formats that the SAM-T99 web page offers (a2m, prettyalign, html, sequence logo) and even a2m_with_dots. In the HTML format, we should show what sequences are structurally aligned by FSSP and give DALI Z-scores, We should also link any PDB sequences to their entries in PDB, FSSP, DSSP, SCOP, ... . We should give DSSP and STRIDE secondary structure in the alignment part of the HTML file.

The interface should be easy to use, with easy maintenance as new multiple alignments are added. There should probably be a queue of old alignments to rebuild, to take advantage of updates to the sequence databases. When an alignment more than x months old is retrieved by a user, it should be queued for rebuilding. Instead of maintaining a consistent FSSP database, we might want to update TABLE2 and TABLE3 weekly, but only update the FSSP files themselves when an alignment is rebuilt.

It would also be a good idea to provide some way to download the entire set of alignments using FTP. This would probably require a script that uses tar to bundle all the existing alignments into a single file and gzips it. Whether this script is called on demand or is run as a cron job on a regular basis depends on how frequent ftp access for the whole set are expected to be. There should probably also be FTP access to individual .a2m files.

We will soon have other such large libraries (such as predictions for the whole human genome). It would be good for the mechanism to be easily modified to incorporate different libraries.

Non-protein projects

Tutorials

If you have very strong writing skills and would rather write English than programs, you can write a few tutorials as your class project. The BME 100 (undergraduate bioinformatics class) still lacks an adequate textbook.

Here are some possible topics (with rough size in lectures for each), though naturally not all these topics could be covered in one course:

1. 1 Review of fundamental dogma of biology, DNA bases, RNA bases, Amino acids
2. 2 DNA sequence databases on the web: genomic, unfinished genomic, EST, clustered EST. Specific databases include Genbank, NCBI's NR, dbEST, and (Unigene, STACK, EGAD, or other clustered EST database) Error rates in different DNA databases due to sequencing error
3. 2 Protein sequence and structure databases on the web: Swiss-Prot, TrEMBL, NCBI's NR, PDB, FSSP, SCOP, CATH Visualization tools: RASMOL or CHIME
4. 3 Pairwise alignment: substitution matrices and dynamic programming---both Needleman-Wunsch and Smith-Waterman algorithms
5. 3 Probabilistic models: Bayesian statistics, E-value, null models This will cover what E-values mean, but not the theory for computing them.
6. 2 Fast search: BLAST algorithm and interpreting BLAST results. A cursory overview of the BLAST algorithm (not enough to re-implement it), with more thorough coverage of the uses of different options to BLAST. Visualization tool: NCBI's pictorial BLAST results.
7. 4 Multiple sequence alignment: clustal, HMMs, regularizers Information and relative entropy Visualization tools: belvu (or equivalent) and sequence logos
8. 4 Protein structure prediction: HMMs, threading, 2ry structure prediction, physics-like ab initio prediction
9. 3 DNA chips: how they work, design of oligos, extracting data from images
10. 3 mRNA expression analysis: clustering of profiles
11. 4 Genefinding: HMM structure, signals and content sensors, using EST and homology data, using cross-genome similarity
12. 4 RNA structure prediction and alignment: energy minimization (Zucker), correlated mutation, stochastic context-free grammars
13. 6 Phylogeny: parsimony, max likelihood, evolutionary models
14. 2 DNA assembly
15. 2 X-ray crystallography---how are models created
16. 3 NMR spectroscopy---how are models created
17. 2 2D gel proteomics
18. 2 Mass spec proteomics
19. 3 Finding binding sites in DNA regulatory regions

RNA alignment

At UCSC, we have used stochastic context-free grammars (SCFGs)to do multiple alignments of large RNA molecules (particularly the ribosomal RNA molecules). Unfortunately, the full dynamic programming on SCFGs is very slow. Michael Brown came up with one heuristic solution, using constraints derived from an HMM alignment to bound the dynamic programming. This approach worked well, but required a lot of work to implement.

A somewhat simpler approach has been proposed based on the observations of Michael Tanner for similar problems in decoding complex codes. Instead of doing optimum decoding, you are often better off using a better code (that is, a better model of your data) and using a non-optimal decoding scheme. The iterative scheme that Tanner has used is quite simple to implement and can be easily translated into the sort of information we get from an RNA secondary structure model (like an SCFG).

Below this line are LAST year's potential projects. Some of them may be still available, but check with Prof. Karplus.

Drosophila indexing

Bill Sullivan (Department of Molecular, Cell and Developmental Biology) suggests that it would be useful "to have a simple read of the classes of genes by cytological location--for example, what are the locations of all the serine/threonine kinases, or motor proteins, ..."

I don't know whether the Drosophila genome is fully annotated for cytological location (if it is, this would be simple query construction). If it is not (the more likely case), then location prediction would be needed (of course, location information would need to be marked to indicate the degree of certainty). There has been a fair amount of literature on predicting the cytological location of proteins, but there is still a lot that can be done, but in determining what features to consider and in determining what classification algorithms to use.

ATPases

Lincoln Taiz has suggested "something I think needs doing is to carry out an analysis of the various minor subunits of the F- and V-ATPases. (There are about 10 subunits altogether). Crystal structures are still lacking for the holoenzymes. Although functionally analogous, many of the subunits seem to be only distantly related by sequence analysis. Presumably their structures are conserved. Can students detect these structural features?"

This looks like a good application for the protein sequence analysis that Bioinformatics I tries to teach. There are at least 54 PDB chains identified as ATPase, and I'm not sure which of them are of interest. I suspect that bovine mitochondrial f1-atpase (FSSP reps: 1skyB [alpha subunit], 1skyE [beta subunit]) is one of the ones of interest. These are both 3-domain proteins, with SCOP domains 1.69.1.1.3, 2.46.1.1.3, and 3.31.1.10.5. The PDB file 1bmf has a different quaternary structure for the proteins.

Other FSSP representatives (and SCOP domains) of chains that have "atpase" in their names are 1a91 (6.2.1.1.7), 1aw0 (4.48.16.1.2), 1b8xA (1.47.1.1.19, 3.42.1.5.19), 1ba1 (3.50.1.1.1), 1bmfG (1.20.1.1.1), and 1cz4A (2.49.2.3.3, 4.27.1.1.3). I don't know if any of these are relevant, though 1bmfG clearly interacts with the alpha and beta subunits (indeed, it is part of gamma-f-atpase).

See also some earlier work on v-atpase in /projects/compbio/experiments/protein-predict/bowman/v-atpase Note particularly the importance of using reverse-sequence null modles because of the long ampipathic helices in these proteins.

Ion channel (nanopore)

Last year some people worked successfully on data analysis from the ion-channel nanopore project. There may be opportunities for more work there. Contact Mark Akeson and Stephen Winters for more info. The nanopore itself is PDB file 7AHL, and is well worth looking at. The barrel that goes through the membrane is a very pretty one, with charged groups at both ends and no charges (inside or outside) for the length of the barrel.

Evolutionary (phylogenetic) trees

It is often handy to examine a set of sequences and group them by similarity. Evolutionary trees are a popular way of expressing the relationships among sequences. We do not use evolutionary trees much in our tools---not because they aren't useful, but because it is hard to make good ones, and bad ones can be worse than useless.

I have one program (phytree) which is my reimplementation of a somewhat crude phylogenetic tree method developed by Kimmen Sjolander. Despite its rather crude trees, it has been valuable in helping select templates for protein structure prediction, when there are several templates from the same superfamily available.

This program could stand quite a bit of work (for example, eliminating the crufty readseq input routines and using zlib to read compressed a2m alignments). It might also be pedagogically valuable to try implementing some of the standard phylogenetic tree algorithms (such as neighbor joining).

A more ambitious student could try to implement some of the more modern ideas on phylogenetic trees. The Symposium on Competing Technologies for Phylogenetics (SCOPH) (April 19-21, 2001 Universite de Montreal, Montreal, Quebec) lists the following:
Olivier Gascuel, Distance based methods
Tandy Warnow, Fast converging methods
David Bryant, Quartet and combinatorial methods
Andreas Dress, Networks and split decomposition
Kevin Nixon, Parsimony
David Sankoff, Gene order parsimony
David Swofford, Maximum Likelihood
John Huelsenbeck, Bayesian methods
Mike Steel, General stochastic models; Hadamard conjugation
Tom Hagedorn, Phylogenetic invariants
Joe Felsenstein, Validation methods

The meeting is fairly cheap ($50 for students), so if you can get to Montreal just before RECOMB 2001 and are interested in phylogenetics, it looks like the place to be.

SoE home

Kevin Karplus's home page

Biomolecular Engineering Department

Karplus's lab page

UCSC Bioinformatics research

Questions about page content should be directed to