Simple Markov chains
Due Friday 2 Nov 2007 (11 a.m.).

This assignment is intended to deepen your understanding of first-order Markov chains and of relative entropy. I will provide you with a data set consisting of many protein sequences in a FASTA file (see Perl assignment 1). You will build stochastic models of these sequences as zero-order and first-order Markov chains, and measure the information gain of the first-order model over the zero-order model.

I have provided two FASTA files: perl2_f07_1.x-seqs and perl2_f07_2.x-seqs for training and testing the models. You will do 2-fold cross-validation (training on each of the sets and testing on the other). These can be accessed on SoE machines directly from the directory /cse/classes/bme205/Fall07/ so that you don't need to make your own copies.

The combined set was derived from Dunbrack's culledpdb list of X-ray structures with resolution 1.8 Angstroms or better and R-factor <= 0.25, reduced so that no two sequences had more than 30% identity. I used only those sequences that were in our t06 template library, which primarily means that short sequences were removed. I randomly split the list of chain ids into two partitions. I also x-ed out the lower-case letters of the sequences, so that amino acids not present in the 3D structure file are not listed. For full details on how the set was created, see the Makefile.

Note: this is a different train and test set than the ones used in previous years---make sure you get the right set.

For simplicity, we will not be modeling the length distribution of the strings---instead we will build length-dependent stochastic models. That is, our models will be probability functions that sum to one for strings of a given length, rather than for all strings.

The simplest model we discussed in class is the i.i.d. (independent, identically distributed) model, or zero-order Markov model. It assumes that each character of the string comes from the same background distribution P_0, and the probability of a string is just the product of the probabilities of the letters of the string.

Step 1: Write a program that reads a FASTA file from STDIN and outputs the counts of amino acids in the sequences to STDOUT. This output should be both human- and machine-readable, as you will be using it later. Document the output format!

Run this program on each of the sequence sets to get 0-order models.

Step 2: Write a program that reads in a table of background counts (in the same format as you output them) from a file and FASTA-formatted sequences from STDIN, and computes the encoding cost in bits for each sequence in the file. The encoding cost of a character x is -log_2 P_0(x), and the encoding cost of a sequence is the sum of the costs for the characters (note: this is equivalent to taking -log_2 of the product of the probabilities).

Note that you need to convert the counts into probabilities. The simplest approach is to normalize the counts to sum to 1. Because the FASTA training file will be large, you will not need to worry about pseudocounts---this simple maximum-likelihood estimate will be fine.

Produce a table that reports for each sequence the ID, the encoding cost for that sequence in bits, and the average encoding cost per character (encoding cost/length). Also report for the whole data set the total encoding cost, the average encoding cost per sequence, and the average encoding cost per character.

Run this program on both sets of sequences using both models. This should give you 2 files of results on training sets and 2 files of results on test sets. Name the files sensibly, so that you can tell which result is which!

Step 3: Write a program that reads a FASTA file from STDIN and builds a first-order Markov model of the sequences. Since we are not modeling the length distribution, there will be no stop state, but we will have a start state, so a first-order Markov model consists of 21 probability distributions over the 20 amino acids. For this assignment, output counts for each of the 21 conditions, rather than probabilities. Document the output format!

Run this program on both sets to get two first-order models.

Step 4: Write a program that reads a first-order Markov model from a file (in the same format as you output for Step 3) and FASTA-formatted sequences from STDIN, and computes exactly the same table as for the zero-order Markov model. Note: it is possible to have a pair of characters occur in the test set that had zero counts in the training set. Your program should not provide infinite results in this case. Document what you do to avoid infinity.

Run this program on the training and test sets and compare the average encoding cost of the first-order model with the zero-order model.

To turn in

• Combine the two training programs into one that provides both zero-order and first-order models (in one file or in two files).
• Combine the two testing programs into one that evaluates both zero-order and first-order models and make the table have the difference in encoding cost (and difference in per-character encoding cost) as new columns in the table.
• Note: it is possible to have one program that does all 4 tasks, but make sure that the count outputs are available as a file, as these may be useful for other tasks.
• Compute the relative entropy between each of the 21 probability distributions for the possible prior states (start and 20 amino acids) and the background distribution. Where are the biggest information gains coming from?
• Measure the "reversibility" of a set of sequences by trying to encode the reversed sequences with the Markov chain. How much more cheaply encoded are they in the forward direction than the reverse direction?
• Instead of just using the amino-acid alphabet, allow the user to specify some other alphabet at run time---for example "-alphabet ACGT". I have provided a couple of other sets, one using the str2 secondary-structure alphabet, and one using the near-backbone-11 burial alphabet: perl2_f07_1.str2s perl2_f07_2.str2s perl2_f07_1.near-backbone-11s perl2_f07_2.near-backbone-11s The str2 alphabet includes the following "normal" characters: ABCEGHMPQSTYZ. The character X is also possible (as a wildcard) and I (as an alias for H). The near-backbone-11 burial alphabet has letters ABCDEFGHIJK.
• To get even fancier with alphabets, allow the user to specify the alphabet by name, and look it up in the /projects/compbio/lib/alphabet/ files we have defined for our various other programs. The ExtAA alphabet gives the definition of the amino acids. Note the specification of wildcards. Another alphabet we use a lot is the EHL2 alphabet, in DSSP.alphabet---note its use of aliases to collapse several different letters in a sequence into the same equivalence class.
• Use higher-order Markov chains as models. Note that each time you increase the order of the Markov chain, the number of bins you have to count into goes up by a factor of 20, which will soon result in very small counts, necessitating a better way to estimate the probabilities than the simple maximum-likelihood estimate. One good method would be to use the probabilities of the next lower-order Markov chain as pseudocounts.

Hints

• Share the FASTA-reading routines between all your programs. You may want to put the routines into a PERL module for easy access in different programs, and to avoid having to correct the same bugs in multiple places.

  You need to think about about upper and lower case (they should be treated as equivalent) and about letters that are not the 20 standard amino acids (wildcards B,Z, and X, and spaces "." and "-").

  • X is a code for any amino acid
  • B is a code for aspartic acid or asparagine (D or N)
  • Z is a code for glutamic acid or glutamine (E or Q)
  The B and Z codes occur mainly in protein sequences that have been determined by old-fashioned chemical methods, not from translating DNA sequences.

  There are several ways to handle wildcards and non-standard amino acid letters:

  1. treat them the same as spaces, tabs, ".", and "-". That is ignore them completely on input. This will make one erroneous transition: ABF would look like AF.
  2. Arbitrarily pick one of the possible amino acids (maybe replace B by D and Z by E).
  3. Treat the character as a gap in the sequence---handle the subsequences on each side as different sequences. There are two variants of this---one which starts a new sequence after the B like any other sequence (so ABF would have a start->F transition) and one that doesn't add the start->F transition. Note that adding a start->F transition may throw off probabilities for what the first residue of a sequence really is.
  There are undoubtedly other methods also. Method 3 with no extra start transition is probably best for this program, but any of the methods is acceptable. You must document how you handle B and Z.

  Think about the FASTA-reading requirements of the first-order Markov model first---it needs to process each sequence separately, so you might want to write a routine that reads one sequence from the file, and returns it as a single string, returning "undef" if there are no more sequences to read. (Don't use the empty string to mark the end of the file, since an empty string is a legal sequence, and actually occurs fairly frequently in real data sets.)

• For representing a count or a probability distribution over 20 amino acids, PERL's hash type is convenient. You can use the amino acid letters as keys into the hash. For output, you can use commands like

      foreach $letter (sort(keys(%aa_prob)))
      {   print "$letter\t$aa_prob{$letter}\n";
      }
  

• Because real sequences may have length 0, be sure your program can handle the empty sequence properly, without choking over a divide-by-zero error. It doesn't matter what you report as the average encoding cost for a zero-length string, just that it not result in errors.
• Strings can be manipulated in perl, but when you want to get at the individual characters, it is often better to split the string into an array of characters with

      @chars = split(//, $seq);
  

• The reversal of an array is standard perl function, as is the reversal of a string:

      @rev = reverse(@chars);
      $rev_seq = reverse($seq);
  

• An observation from previous years: A lot of students managed to make their programs more complex to write and harder to use by insisting on file I/O instead of using STDIN and STDOUT. Try to avoid that this time. Also, error messages should be printed to STDERR (to avoid corrupting the main output in STDOUT). You may want to use the "die" function in perl for fatal errors.

Things learned for next year

• Remember that Markov chains are conditional probabilities P(x_i=s | x_{i+1}=t), not joint probabilities.
• It is cleaner to use a hash with a string as a key, rather than a hash of hashes with letters as keys.
• Use "next", "last", "return", "die", ... to reduce the depth of indenting needed in your code. When choosing the polarity for an if test in an if-then-else, see if one of the blocks ends with a statement that does not continue sequentially, and make it the then block. The "else" can then be eliminated, reducing the depth of the code. If both blocks fall through, then put the shorter, simple block as the then-block.
• Break code into paragraphs, separated by blank lines, and start each paragraph with a comment as a topic sentence for the paragraph.
• Don't read all the data before processing any. If you want to make a fasta-reading module, make sure its interface returns one sequence at a time.
• A sequence may span many lines. Make sure that you don't start a new sequence with each new line.
• Don't have a logarithm computation in the inner loop. Convery you counts into costs before running through the test set.

SoE home

Kevin Karplus's home page

BS, MS, and PhD programs

BME 205 home page

Karplus's lab page

UCSC Bioinformatics research

Questions about page content should be directed to