Due: Fri. Nov. 22

SAM HMM system, PHD program, detecting remote homology

Exercises provided by Christian Barrett, with additions by D. Haussler

This assignment is intended to show you some of the basic features and capabilities of our HMM software, which is called SAM for Sequence Alignment and Modeling software. It will also give you a chance to use PHD, which is Burkhard Rost's neural net that predicts secondary structure of protein sequences. These will be demonstrated in the context of the Archae whole genome analysis project that was just started this quarter, using the M. jannaschii genome.

One aspect of this project focuses on functional annotation of the TIGR-predicted proteins from the M. jannaschii genome. This method begins with a sequence that the TIGR group suspects is a protein, based on analysis of the genome with the GENMARK genefinding program. We take this protein and build an HMM from it. Normally this is really not enough information from which to estimate all of the parameters of an HMM, but it's all we have at first. Hopefully the Bayesian method of estimating the parameters of the HMM using Dirichlet mixture priors will help, as mentioned in the lecture.

One thing we can do with this HMM is to search for related proteins within the M. jannaschii genome. It is often the case that a genome contains families of closely related genes. Thus, given one M. jannaschii protein, we will often find other similar proteins among the remaining M. jannaschi proteins. If we find these, then we can build an HMM for the entire family of related genes. Thus the next step is to search the M. jannaschii genome using the initial model built above form a single protein sequence. If we pull out any sequences that score well with the model, a new HMM is built using these sequences, along with the original sequence.

This has been done with putative protein sequence MJ0553 from M. jannaschii. Make a directory for yourself and copy all of the files from /projects/compbio/class/cmp243/hw5 into it.

One of these files is "MJ0553_2_long.dist". This is the file that contains the log-odds scores for the potentially interesting sequences in the genome. (In class we've been using using base 2 logarithms, as right thinking people ought. SAM uses base e logarithms, for "historical reasons"). (Those comments by C. Barrett. I take no responsibility! D.H.)

If you look in this "distance file", you will see some NLL-NULL scores, and above them some comments about significance. (Ignore the NLL-NULL/L scores.) First let's look at how the NLL-NULL scores are computed. For example, the 4th line gives the score for the protein sequence MJ1405. which is -123.356. Using "ln" to denote natural logarithm, this score is

   -ln (P(MJ1405|HMM)/P(MJ1405|null)) + ln (number of residues in MJ1405),

Normally, the quantity P(MJ1405|HMM) would be calculated by the forward algorithm that we studied in class. However, in this case, the actual algorithm that was used was a fast version of the Viterbi algorithm that does a Smith-Waterman-type local alignment of the sequence to the HMM. What is meant here is that the dynamic programming algorithm is set up to try to match part of the sequence with part of the HMM, rather than always matching all of the sequence to all of the HMM. The output of this method is the probability of the sequence, jointly with the best local alignment of this type. So this is really at best an approximation to the quantity P(MJ1405|HMM) as defined in class, but it has other advantages, besides speed of calculation, which we will not go into.

To understand the adjustment term "+ ln (number of residues in MJ1405)", we need to review how significance is interpreted in terms of log posterior odds. By Bayes rule, the log posterior odds are given by

   ln (P(HMM|MJ1405)/P(null|MJ1405)) =
   ln (P(MJ1405|HMM)/P(MJ1405|null)) + ln P(HMM)/P(null)

     - ln (number of residues in MJ1405 * number of residues in the database)
     = - ln (number of residues in MJ1405) - ln (number of residues in the database)

The first term is the adjustment. This is subtracted by SAM from all the negative log likelihood scores, to get the NLL-NULL score. The second term is a constant that does not depend on what sequence is being analyzed. It is mentioned in the discussion of significance, before the scores are given in the SAM MJ0553_2_long.dist file, rather than being added to all the scores. In this search there are 63467 residues in the database, and this term is the constant -ln(63467)= -11.1. Putting this together, we see that the negative log posterior odds are

   - ln (P(HMM|MJ1405)/P(null|MJ1405)) =
   -ln (P(MJ1405|HMM)/P(MJ1405|null)) + ln (number of residues in MJ1405) + 11.1 =
   NLL-NULL score + 11.1

For the match to be significant, this value should be negative. Specifically, if you want the posterior odds to be greater than 100 to 1, then the -ln posterior odds should be less than -ln(100) = -4.6. This calculation is shown in the first line of the SAM output listing of confidence cutoffs. In subsequent lines you see confidence cutoffs computed using different assumptions. We won't go into these calcs here.

Q1: Using the above calculations, how many sequences in the database have NLL-NULL scores that indicate that the posterior odds are at least 100 to 1 for the HMM model vs the null model? Repeat the above calculation assuming that the database has 100,000,000 residues in it. What is the cutoff for significance? How many of these sequences are still deemed significant at odds of 100 to 1?

Before leaving this topic, we must note that this approximation ignores the important fact that, unlike BLAST matching, we allow indels in HMM matches, which makes for a lot more possible "places" for the match to occur. Because of this, we do not have a rigorous mathematical theory of the significance of HMM scores; nobody has been able to solve the difficult mathematical problems involved in this. However, this log odds interpretation at least gives some heuristic guidance.

Since we don't want to answer the above questions for you, let us assume now that we use these scores as the simple heuristics that they are, and just look for a break in the score distribution. Upon examination of the "MJ0553_2_long.dist" file, we see immediately that there are three other sequences in M. jannaschi's protein database that score markedly better than the rest. Given this clean break in scores, we will take these three new sequences (MJ0809, MJ1331, and MJ1405) and reestimate the HMM for MJ0533 with them. (It is interesting to note that none of these sequences have any annotation in the M. jannaschi protein database. So bear in mind that if we can find a homology match to a protein with known function, or even better, known structure, then we will be the first people *ever* to know what these proteins do and look like.)

At this point we have as good a model as we can get, given the information available. The next step is to search SWISSPROT. Our aim here is to find sequences in this database that score well with our model. SWISSPROT is probably the most carefully constructed protein database, in that the sequences it contains are carefully annotated with descriptions, journal article links, other database links, and important residues (such as active sites, critical to structure or function, etc). Thus, if we can find a hit in SWISSPROT using the model built from our sequence of interest in the M. jannaschi genome, then we may be able to infer quite a bit about it. This was tantalizingly noted above.

So search SWISSPROT we do. The file "MJ0553_2_sprot33.dist" lists the top- scoring sequences from the 52,205 that are in the current release of SWISSPROT. Hmmm. Kinda on the low side of our significance range.

Q2: Repeat the same calculations as in Q1, but using the larger database size of 18531282 residues for SWISSPROT, and see how many sequences in this database search have even posterior odds (1 to 1) or greater in favor of the HMM vs the null model.

Since these scores are not convincing, we want to perform a visual inspection of the multiple alignment between these SWISSPROT sequences and the four homologous sequences from the M. jannaschi genome that we have from the search described above. For this you will use SAM's WWW interface at

     http://www.cse.ucsc.edu/research/compbio/sam_form.html

Before going on to the next step, lets see what these SWISSPROT hits are. Using your web browser go to the SWISSPROT search page at URL

      http://expasy.hcuge.ch/cgi-bin/sprot-search-de

For each PDB entry (which contains the atomic coordinates for a protein), there is the raw amino acid sequence for the protein in the Brookhaven database. In case this particular sequence isn't in SWISSPROT yet, we'll get this raw sequence from Brookhaven. It is in the file "1aps.seq".

So, go to the SAM web page SAM web page and fill in accordingly:

  Buildmodel - yes, data back - no
  Send back results of align2model only
  You have not provided an initial model
  Fill in your email address
  For training sequences, cut and paste in the sequences from the files
     1aps.seq, topSP-cleaned.seqs and MJ0553_2_long.sel
  For buildmodel parameters:
    "minmodlen 98 maxmodlen 104 Nmodels 4"
  Prettyalign - yes
  For prettyalign parameters, "-m 6"

In the next part of this assignment, we will try the PHD program. What we would like to do is run this program on the putative protein MJ0553, and see if the secondary structure it predicts shows any similarity to the known secondary structure for 1aps. If it does, then this is further evidence that MJ0553 may be homologous to 1aps. If it does not, then they might not really be homologous, or PHD may just be predicting the wrong secondary structure for MJ0553.

Go to the WWW resources page for cs243 and click on the PHD web server. Notice that there are a number of prediction programs. Click "help" for some introduction. After scanning this, go back to the main PHD page and click "submit" to submit the MJO553 sequence. Name your prediction "MJO553". Paste in the MJ0553 sequence and click "RUN PREDICTION". It may take several hours to get an email response. Just as a quality check on the PHD program, repeat the above steps with the 1aps sequence. In this case, we know the secondary structure, and we can check how well PHD performs.

Now while you are waiting for PHD predictions to be mailed, let's look at the HSSP file for 1aps, and see if we can figure out what its secondary structure is. At the top of the file are several lines that describe how the file was built, and what the notations in the file mean (lines beginning "NOTATION"). What is interesting about the HSSP file is that it also includes known homologs of the protein 1aps, and shows how these can be aligned to 1aps. In this case there are 14 homologs listed in the block after the "NOTATION" section.

Q3: How does this list of homologs differ from those in topSP-cleaned.seqs?

After this list of homologs, we get to the meat of the HSSP file (starting at "##ALIGNMENTS ..."). On the left side are the 14 homologs of 1aps aligned vertically. Reading from the right side to the left, there are several columns of information about 1aps. They are

1. SeqNo: reside number in the 1aps sequence (from Brookhaven database). Goes from 1 to 98.
2. PDB residue numbers: In this case they are the same as SeqNo, but other times PBD uses weird an flakey numbers (like -1 1 2 3 3A 3B 4 5 11 12, etc.)
3. AA: the amino acids of 1aps (vertically).
4. STRUCTURE: This is what we want. The first letter of this field is the secondary structure. As I mentioned in class, "H" represents alpha helix, and "E" represents beta strand. Other letters you see are blank, "T" and "S". These can be considered as belonging to the third class "other" for the purpose of 3-state secondary structure prediction. The remaining entries in this field give structural details we won't go into.
5. and the rest: BP1 and BP2 refer to pairing among the beta strands to form bets sheets, ACC is a measure of how exposed or buried the residue is in the structure, NOCC is the number of homologs whose alignment includes this position, and VAR is a measure of the sequence variability in this position among the homologs, i.e. how conserved this position is.

Finally, at the bottom of the HSSP file is a raw profile for this protein family built from the multiple alignment, giving amino acid probabilities as percentages. No fancy statistical estimation is done here, so these are just the raw frequencies.

Take the file 1aps.seq containing the 1asp sequence in the usual horizontal form, put it all on one line with no spaces, and create a line below it showing the known secondary structure for this sequence, using the letters "H", "E", and "L" (for "other"). Now when you get the PHD results for 1aps, read the file to figure out what is going on, then add another line to this file indicating the PHD secondary structure prediction. Turn a copy of this in with your homework.

Q4: How well does PHD do? How many residues have correct secondary structure prediction? Be aware that it is very difficult to predict the residues at the ends of beta strands and alpha helices, indeed there is some disagreement on how to classify these occasionally, even given the 3d coordinates. Look also at scores, on a scale from 0-9, that PHD gives for its predictions.

Now let us look at the PHD prediction for MJ0553.

Q5: In the critical region in the SAM alignment of MJ0553 and 1aps, along with their homologs, (say from alignment positions 8 to 98), how often does the PHD predicted secondary for MJ0553 agree with the known secondary structure of 1aps?

The last part of this assignment is for you to decide if MJ0553 (and some/all of its M. jannaschi homologs) are homologous to this family of SWISSPROT sequences. This is Q6. Justify your decision.

Here are some thoughts and hints (from Christian):

1. On first inspection, this alignment looks good because the structure sequence (1aps) aligns without any deletions or inserts to the rest of the sequences.
2. Get a colored pen or two and a straightedge. Color the completely conserved columns. Now color the columns that are plausible in that all of the amino acids seen in that column could be substituted for each other. That is, mixtures of similar amino acids like I,L,V and A,S,T or F,Y,W or D,E,N or K,R, etc. It would be best to review the Tooze chapter and Gregoret lecture notes on how the amino acids are grouped and categorized (i.e. hydrophobicity, hydrophilicity, polar. charged, size, etc.) Some other pointers:
   1. the residues that are the most important functionally tend to be certain polar ones: D, N, S, Q, E,
   2. certain structural ones are also sometimes crucial to maintain conserved patterns : G, C, W, P
   3. charged residues should be matched with similar polar residues eg, E -> Q or D->N, unless you can substitute a similar charged residue
   4. sometimes a positive residue like K is substituted (in another family having the same fold) by a negative residue (D or E), and vice versa this can be simply due to interacting with water molecules, or some kind of substrate binding
3. Look in the SWISSPROT entries for annotation indicating structurally or functionally conserved/important residues. Are the M. jannaschi sequences in agreement, or at least plausible in their substituted amino acids?
4. One of the best ways to determine if this multiple alignment makes sense is to check out the alignment of these sequences to the 3D structure. But that is probably taking this assignment a little to far.

Back to the CMP 243 Class Page.