CS498CXZ Assignment #4: Pairwise Sequence Alignment
(due Oct. 11, 2005, Tuesday, 12:30pm)

1. [60 points] Consider the sequences v=CGATAAC and w=ACGTTAC. Assume that we score an alignment by awarding a match with"+1" and penalizing a mismatch and indel both with "-1".
   
   
   • [25/60 points] A. Global alignment : Fill out the following dynamic programming table for a global alignment between v and w. Draw arrows in the cells to store the backtrack information. Draw multiple arrows to indicate ambiguity in backtracking when you have multiple ways to obtain a maximum score from an alignment in the prefixes of the two sequences. What is the score of the optimal global alignment? Which global alignment(s) does this score correspond to?
     
     

     	-	A	C	G	T	T	A	C
     -
     C
     G
     A
     T
     A
     A
     C

     
     
     
   • [20/60 points] B. Local alignment : Fill out the following dynamic programming table for a local alignment between v and w. Draw arrows in the cells to store the backtrack information. Draw multiple arrows to indicate ambiguity in backtracking when you have multiple ways to obtain a maximum score from an alignment in the prefixes of the two sequences. What is the score of the optimal local alignment? Which local alignment(s) does this score correspond to?
     
     

     	-	A	C	G	T	T	A	C
     -
     C
     G
     A
     T
     A
     A
     C

     
     
     
   • [15/60 points] Affine gap penalty: Suppose we use an affine gap penalty where it costs -20 to open a gap, and -1 to extend it. Scores of matches and mismatches are unchanged. What is the optimal global alignment in this case? What score does it achieve?
   
   
   
2. [30 points] The BLOSUM scoring matrix assigns a probability score for each position in an alignment that is based on the frequency with which that substitution is known to occur among consensus blocks within related proteins. BLOSUM62 is among the best of the available matrices for detecting weak protein similarities, and is the default matrix used in BLAST 2.0, one of the most popular local sequence alignment tools. The value in a BLOSUM matrix corresponding to amino acids X and Y is computed as

                p(X,Y)
   val(X,Y)=log ----------
                p(X) p(Y)
   

   where p(X,Y) is estimated as

             number of times X aligned with Y
   p(X,Y)= -----------------------------------
             total number of aligned pairs
   

   and p(X) and p(Y) can be estimated as

          number of times seeing X in any position of a sequence
   p(X) = -------------------------------------------------------
          total number of positions in all the sequences
   
          number of times seeing Y in any position of a sequence
   p(Y) = -------------------------------------------------------
          total number of positions in all the sequences
   

   A brief note about how to compute a BLOSUM matrix is available here.
   • [10/30 points] A. Suppose we only consider 3 amino acids ( A, B, and C) and have observed the following 3 sequences in a conserved block (BLOCK1):

         BLOCK1 
     
     Seq1 = A  B  C  A
     Seq2 = A  B  B  B
     Seq3 = A  B  B  C
     

     We will have 3 possible pairwise alignments between these sequences and a total of 3x4=12 alignment pairs. Compute the BLOSUM matrix for these 3 amino acids.
   • [10/30 points] B. Now suppose we have a fourth sequence Seq4, which happens to be a duplicate of Seq2 and thus have the following block (BLOCK2):

         BLOCK2
     
     Seq1 = A  B  C  A
     Seq2 = A  B  B  B
     Seq3 = A  B  B  C
     Seq4 = A  B  B  B
     

     Compute the BLOSUM matrix for BLOCK2. Show all the computations.
   • [10/30 points] C. Compare the two matrices you have obtained. How are they different? Which amino acid pairs have higher values for the matrix of BLOCK2? How is this related to the addition of Seq4?
   
   
   
3. [10 points] BLAST is one of the most popular tools for sequence alignment. Biologists often use it to find the proteins/genes that are similar to one they are working on. The purpose of this exercise is to obtain some experience with using BLAST and to understand the meanings of the reported results through searching for protein sequences that match the sickle cell hemoglobin protein sequence. The sickle cell hemoglobin is a mulfunctioning protein due to a change in ONE nucleotide in the DNA sequence that leads to a change in ONE amino acid that changes how the hemoglobin protein folds. This change in the structure of the hemoglobin protein leads to a change in the shape of the red blood cell to a sickle shape, causing diseases as explained in this website .
   
   
   • First, please familiarize yourself with BLAST by reading this tutorial page.
   • Obtain a sickle cell protein sequence: Visit the National Center for Biotechnology Information (NCBI) site. Select "protein" database and enter "liganded sickle cell hemoglobin" in the search box. Choose the 3rd record titled "1NEJD" in the result list and click on "1NEJD". This will bring you to a PDB record. Under the "Display" menu, select "FASTA". (The default is GenPept.) You should now see a protein sequence displayed in captalized letters, each representing one amino acid.
   • Use BLAST: Open a separate window to visit NCBI website. Select BLAST from the horizontal menu bar on the top. You will see there are many variants of BLAST for different purposes. Click on "Protein-protein BLAST (blastp)" (blastp is for protein alignment while blastn is for nucleitide alignment). Cut and paste the sequence of the sickle cell hemoglobin protein to the big search box next to the keyword "search" in the blastp search interface. Click "BLAST!". When the results appear scroll down a little and Click the "Format!" button. You will see the results of local alignment returned by blastp.
   • Understand the results: Scroll down the result page to see all the alignments. You should be able to see that the alignment generally becomes worse and worse as you scroll down to see more alignments. blastp reports 5 figures for each alignment: The first two are "bits" and "Expect", which are significance values. A larger "bits" value, or equivalently a smaller "Expect" value indicates a statistically more significant match. Specifically, the "Expect" value is meant to be the probability of obtaining such an alignment score or a better score from a randomly generated sequence database. The next three are "Identities", which is the percentage of exact matches, "Positives", which is the percentage of matches with a positive value according to BLOSUM62, and "Gaps", which is the percentage of gaps.
   • Report the results: Write down the five reported figures for the most significant alignment, which is with itself. Write down the five reported figures for the least significant alignment returned. Obtain and write down the accession numbers for these two results by clicking on the corresponding result entry. Now find its alignment with a normal hemoglobin, whose NCBI internal ID is "gi:71727229" and whose accession ID is "AAZ39779". Write down the five reported figures for this alignment and describe how the sickle hemoglobin sequence and the normal hemoglobin sequence differ.
   • Refer to the BLAST guide page at NCBI if you have problems with running BLAST or understanding the alignment results.

What to turn in

Please turn in a hardcopy of your written answers at the class