Write a program that trains a profile hidden Markov model (HMM) M from a multiple protein alignment, and then finds the proteins in your very own personal prokaryote (from Homework #0) that have the highest scoring Viterbi score given M.

You will use the multiple alignment of 17 proteins from Homework #2 results to train the profile HMM M. (See Figure 5.2 of Durbin et al. for the topology of a profile HMM.) Every column of this alignment that has 8 or fewer gap characters will correspond to a Match state in the HMM, and those sequences with a gap in this column will take the path through the corresponding Delete state. Every column that has 9 or more gap characters will not correspond to any Match state, and those sequences with an amino acid in this column will take the path through the corresponding Insert state.

The Viterbi score you will be using to score each of your personal prokaryote's proteins S is the log likelihood ratio of the Viterbi probability of S and the background probability of S. See Equations (5.1) of Durbin et al. for the recurrences to use. You will need to modify this set of recurrences slightly: you need to add a recurrence with V_End(i) on the left hand side (hint: the End state has no emissions, so this recurrence looks very much like the recurrence for the Delete states) and you need to modify the recurrences for the states I₀, M₁, and D₁ (hint: these each have only two inputs, one of which is the Begin state). The basis for the recurrence is V_Begin(0) = 0, V_k(0) = -∞ if k is any Match, Insert, or End state, and V_Begin(i) = -∞ for i > 0. All your logarithms should be computed base 2.

As in Homework #2, you can find a list of all your prokaryote's protein sequences in the .faa file(s) on the prokaryote's NCBI Microbial FTP page. Compute the background frequency q_b of each of the 20 amino acids b from these files. The word "frequency" here means the fraction of amino acids that are b in the file(s); that is, the 20 values of q_b should sum to 1. (Throughout this homework, treat any letter in the .faa files that does not represent one of the 20 amino acids as though it were an A.) You will need these frequencies for the recurrences and for the adjustments described in the next paragraph.

Because there are only 17 characters in each column of the alignment, when training the HMM some emission probabilities would be set to 0, but this is actually an artifact of the small sample size of 17. Instead of using the actual counts E_k(b) of the number of times b was emitted in state k, use the adjusted counts E_k(b) + q_b (in both numerator and denominator) when computing the emission probability e_k(b). This way, e_k(b) will be a small positive number even if E_k(b) = 0. Similarly, use A_jk + Δ_jk rather than A_jk (in both numerator and denominator) when computing the transition probability a_jk, where Δ_jk = 0.9 if j and k are both the same type of state (that is, both Match states, or both Insert states, or both Delete states) and Δ_jk = 0.05 otherwise. If j is the Begin state, use Δ_jk = 0.1 if k is the Match state and Δ_jk = 0.45 otherwise. If j is the last Delete state, use Δ_jk = 0.9 if k is the End state and Δ_jk = 0.1 otherwise. If j is the last Match or last Insert state, use Δ_jk = 0.1 if k is the End state and Δ_jk = 0.9 otherwise.

Finally, for each of the two proteins S⁽¹⁾ and S⁽²⁾ that have the two best log likelihood ratio scores, determine its most probable path through the HMM M by doing the traceback of the Viterbi algorithm. From these two paths, align S⁽¹⁾ and S⁽²⁾ to the 17 other proteins in the input alignment.

Here is an outline of your entire program, to help you get started:

Training:

1. Read the multiple sequence alignment A into an array.
2. Mark the columns of A (the columns with 8 or fewer gap characters) that correspond to match states of M.
3. For each sequence R in A, follow the path of R through M determined by step 2, incrementing appropriate counts of emissions and transitions.
4. Compute background frequencies from the .faa file(s).
5. Compute all emission and transition probabilities.

1. For each sequence S in the .faa file(s), use the recurrences of Equations (5.1) from Durbin et al. to compute score(S) = V_End(|S|).
2. Identify S⁽¹⁾ and S⁽²⁾ that have the two greatest values of score(S).
3. Trace the Viterbi maximizations backwards to find the most probable path for S⁽¹⁾ and the best alignment of S⁽¹⁾ with A, and similarly for S⁽²⁾. (Read the first extra credit problem below for more important details.)

Extra credit:

• Let P be the maximum probability path that one of your best proteins S takes through the HMM, obtained from the Viterbi traceback. Suppose that P goes through some Insert state multiple times and that a few other sequences in the original multiple alignment go through the same Insert state multiple times. Then the traceback gives you no information about how to align this insertion in S with the insertions in the other sequences. Use a separate multiple sequence alignment algorithm just on the corresponding insertions to find good gap placements. You may use some multiple alignment program available on the internet for this purpose if you like. (For the generic HW4 solution it is fine to make any arbitrary decision about gap placement within corresponding inserted sequences.)
• Use the Forward algorithm instead of the Viterbi algorithm and compare the results. In the LLR-transformed Forward algorithm you are faced with the logarithm of a summation, which is not simple to handle.
• Your ideas here.

Turn-in Instructions

You will actually run your program as described above on two separate genomes: once on your personal prokaryote and once on the "community prokaryote" Brucella melitensis biovar Abortus (Be sure to use the background frequencies from the appropriate prokaryote to train and test the HMM.) For each of these two prokaryotes, you will report the two proteins that have the two best log likelihood ratio scores. Turn in the following items:

1. The protein identifier and name for the best-scoring protein S⁽¹⁾, as given in the .faa file, e.g.,

   >gi|15826866|ref|NP_301129.1| chromosomal replication initiation protein [Mycobacterium leprae TN]
   

2. A very brief description of the function of S⁽¹⁾. A good place to start is the PROSITE database, entering the protein name in PROSITE's Search box.
3. The log likelihood ratio score for S⁽¹⁾.
4. The most probable (Viterbi) path for S⁽¹⁾, given as the sequence of state names on that path. If the path goes through an Insert state I_j x times, then I_j should be repeated x times in your path.
5. The protein identifier and name for the second best-scoring protein S⁽²⁾, as given in the .faa file.
6. A very brief description of the function of S⁽²⁾.
7. The log likelihood ratio score for S⁽²⁾.
8. The most probable (Viterbi) path for S⁽²⁾, given as the sequence of state names on that path. If the path goes through an Insert state I_j x times, then I_j should be repeated x times in your path.
9. The protein identifier and name for the best-matching protein you turned in for Homework #2, in the same format as item 1 above.
10. Your explanation for why S⁽¹⁾ is the same protein or a different protein than you turned in for Homework #2.
11. The multiple alignment of the original 17 proteins plus S⁽¹⁾ and S⁽²⁾ as determined by the Viterbi algorithm.

Your turn-in will consist of the following files, named as shown.

1. match.txt:

   [your name]

   Brucella melitensis
   [protein ID and protein name for S⁽¹⁾]
   [brief description of protein S⁽¹⁾]
   [log likelihood ratio score of S⁽¹⁾]
   [Viterbi path for S⁽¹⁾]
   [protein ID and protein name for S⁽²⁾]
   [brief description of protein S⁽²⁾]
   [log likelihood ratio score of S⁽²⁾]
   [Viterbi path for S⁽²⁾]
   [protein ID and protein name turned in for Homework #2]
   [comparison to Homework #2 and explanation]
   [original 17-protein multiple alignment aligned with S⁽¹⁾ and S⁽²⁾]

   [name of your personal prokaryote]
   [protein ID and protein name for S⁽¹⁾]
   [brief description of protein S⁽¹⁾]
   [log likelihood ratio score of S⁽¹⁾]
   [Viterbi path for S⁽¹⁾]
   [protein ID and protein name for S⁽²⁾]
   [brief description of protein S⁽²⁾]
   [log likelihood ratio score of S⁽²⁾]
   [Viterbi path for S⁽²⁾]
   [protein ID and protein name turned in for Homework #2]
   [comparison to Homework #2 and explanation]
   [original 17-protein multiple alignment aligned with S⁽¹⁾ and S⁽²⁾]

2. Source files for your program, with appropriate filenames.
3. README: a short text file explaining how to compile and run your program.
4. any files describing extra credit work, whose filenames should begin with the prefix "extra".

Submit these files to the homework drop box at https://catalyst.uw.edu/collectit/dropbox/tompa/13675.