Project 4.5: Protein Alignment

Preparation

The file chap4.py contains an updated version of the Needleman-Wunsch algorithm, including a number of new procedures that will help you explore Needleman-Wunsch pattern matching.

a. Make your own copy.

b. Look at the code for buildNWTables and explain to yourself (and your partners) how it has changed.

c. Scan through the file to see what other procedures seem to be new.

d. Make copies of pam250.txt and blosum62.txt, too.

1. Hydrophobicity

St. Clair and Visick suggest that one possible improvement to Needleman-Wunsch is to incorporate hydrophobicity.

We currently use four possible actions in building the Needleman-Wunsch table: Deletion, Insertion, Match, and Mismatch. We're now going to add a fifth, Similarity ('Sim').

a. Extend the default action values table and the buildTable procedure to give a score for similarity. (For example, we might use a score of 2 for a match and a score of 1 for similarity.)

b. Build a dictionary that you can use to look up whether or not each amino acid is hydrophobic. St. Clair and Visick provide a table on p. 103.

c. Make a new version of buildNWTables, called buildNWTablesH, that uses the similarity value when it encounters a pair of hydrophobic amino acids or a pair of hydrophilic amino acids. You should focus on the following code as the location for your change.

            if (seq1[col-1] == seq2[row-1]):
                valMat = values[row-1][col-1] + actionValues['Mat']
                sym = 'Mat'
            else:
                valMat = values[row-1][col-1] + actionValues['Mis']
                sym = 'Mis'

d. Update showAlignment to print a symbol (e.g., a colon or a tilde) in the summary of matches row for similar amino acids.

e. Test the two versions of buildNWTables (the old one and the new hydrophobic-sensitive one) on the sequences EREHSISIVLE and QNHKTLGFICN. (Why might it be better to make only one alignment, rather than all alignments, for this pair of sequences?)

2. Aligning Real Data

Try aligning CFTR and the E. Coli protein HisP (accession number P07109) using the basic Needleman-Wunsch algorithms. Follow the instructions on the bottom of p. 106 (last paragraph before Putting Your Skills Into Practice) for getting the sequences.

You will note that Python is unhappy about the depth of recursion necessary to to do oneAlignment. Hence, we need to rewrite it iteratively (without using recursion). The procedure oneAlignmentIterative is intended to do just that. Unfortunately, it has a bug. Fix that bug.

Once you've fixed the bug in oneAlignmentIterative, explore the differences between the alignment algorithms.

3. Incorporating Substitution Matrices

As you may have noted, you can use the readMatrix procedure to read in a substitution matrix given in a simple form. You can use lookupMatrix to look up a substitution in the matrix.

We provide two sample substitution matrices, pam250.txt and blosum62.txt.

a. Verify that you can use these procedures to read in the matrices.

b. Make a new version of buildNWTables, called buildNWTablesM, that takes a substitution matrix as the fourth parameter and uses that substitution matrix to compute the value when it has a pair of non-matching amino acids. (You should use an action of 'Sim' if the value is positive and a value of 'Mis' if the value is negative.)

c. Test your new procedure. For example,

>>> pam250 = readMatrix('pam250.txt')
>>> seq1 = 'EREHSISIVLE'
>>> seq2 = 'QNHKTLGFICN'
>>> (values,directions) = buildNWTablesM(seq1,seq2,actionValuesDefault,pam250)
>>> alignment = oneAlignment(directions, len(seq2), len(seq1))
>>> showAlignment(alignment, seq1, seq2)

d. Write a procedure that encapsulates these steps.

e. Test this procedure on other pieces of data, such as the data from 4.5.2.

f. Try the same alignments using the BLOSUM matrix. What differences, if any, do you see?