Dynamic Programming Exercises. BME 205 Homework #8

Fall 2014
Due Wed 3 Dec 2014
(Last Update: 14:15 PST 7 December 2014 )

There is a nice tutorial on dynamic programming at http://www.sbc.su.se/~per/molbioinfo2001/seqali-dyn.html with pointers to an example at http://www.sbc.su.se/~per/molbioinfo2001/dynprog/dynamic.html. It uses the more expensive alignment algorithm that allows arbitrary gap costs for different length gaps, though only an affine gap cost is used. If you want a textbook presentation, one of the best I know is in Chapter 2 of Durbin, Eddy, Krogh, and Mitchison's book Biological Sequence Analysis. The book An Introduction to Bioinformatics Algorithms by Neil C. Jones and Pavel A. Pevzner is more tutorial, but doesn't cover affine gap costs as well. Chapter 6 is the relevant chapter there.

The version of dynamic programming with affine gap costs that we presented in class used only nearest neighbors, but needed 3 matrices:

Mi, j = max { Ici – 1, j – 1 + subst(Ai, Bj)
Mi – 1, j – 1 + subst(Ai, Bj)
Iri – 1, j – 1 + subst(Ai, Bj)

Ici, j = max { Ici , j – 1 – extend
Mi , j – 1 – open
Iri , j – 1 – double_gap

Iri, j = max { Mi–1 , j – open
Iri–1 , j – extend
Ici–1 , j – double_gap

In Python, one can handle this just as well with one array, containing triples of values—for example, (M,Ir,Ic). Some students find handling one subscripting and a triple of information easier.

Be careful of the boundary conditions. For global alignment, you should have "start"-"start" aligned with 0 cost in the MATCH state and –infinity in Ir and Ic, so that a gap opening penalty is incurred if you start with an insertion or deletion. Before start, all rows and columns are –infinity, and the rest of the start row and column can be inferred from the recurrence relation. Make sure that your initial conditions are consistent with the recurrence relations—people sometimes get the indices of their matrices reversed, or confuse the Ir and Ic matrices.

For local alignment, the recurrence relations and the boundary conditions are different:
Mi, j = subst(Ai, Bj) + max { 0
Ici – 1, j – 1
Mi – 1, j – 1
Iri – 1, j – 1

Ici, j = max { Ici , j – 1 – extend
Mi , j – 1 – open
Iri , j – 1 – double_gap

Iri, j = max { Mi–1 , j – open
Iri–1 , j – extend
Ici–1 , j – double_gap

Since local alignment is allowed to start at any match state (with the 0 choice for previous alignment cost), all rows and columns before the first letters of the sequences can be minus infinity.

You'll need a substitution matrix for this assignment. I've provided BLOSUM62 in the format provided by matblas and used by BLAST.

Q1 (hand computation):

Compute the score for the following global alignment using the BLOSUM62 scoring matrix with a gap-opening cost of 12 and a gap-extension cost of 1 and a double-open cost of 2 (all are in 1/2 bit units to be compatible with the BLOSUM62 matrix). Show what you added up to get the score, since the final number alone could easily be incorrect due to a transcription or addition error.

    --CTNIDDQSA--DK-NR-
    VLDAM-EENIGEGDDIKMV
    
Note: this is neither the global nor the local optimal alignment.

Q2 (align module):

Create a module (align.py) that contains two classes: local_aligner and global_aligner. Both classes should have the same methods:

aligner= local_aligner(subst,open=,extend=,double=)
aligner= global_aligner(subst,open=,extend=,double=)
(the __init__ method) records the substitution matrix and gap costs for the alignments to follow. I found it useful to store the substitution matrix as a dict, mapping strings of 2 letters to the substitution matrix value for that pair. The __init__ method should assign reasonable default values for the three gap costs.
score_a2m(self,s1,s2)
Given two strings that represent an alignment in A2M format, score the alignment with the parameters used when the aligner was created. Use this method to check that the alignments you output score what the dynamic programming algorithm thinks that they should have scored.
score=align(row_seq,col_seq)

Creates and fills in alignment matrices and stores row_seq and col_seq.

I found it convenient (though slow) to use dict objects for the matrices, mapping (i,j) 2-tuples to the value of the matrix at that coordinate pair. Note that this allows negative indices. You can make things faster by using numpy arrays with normal 0-based indexing of arrays, but you have to be careful about the bounds.

seq=traceback_col_seq()
Treating the letters of the (cached) row_seq as alignment columns, return the col_seq sequence with appropriate lower-case and "-" characters so that it would correctly provide the alignment in an A2M file (col_seq residues that do not align should be lower-case, and there should be a "-" every time a character of row_seq is skipped).

You may implement the recursive definition of the computation of the alignment matrix in any way that works, as long as subalignments are not recomputed. Although dynamic programming is the standard implementation (being fast and fairly simple), some people find it easier to do memoized computation. For example, you might have three methods m(i,j), ir(i,j), and ic(i,j) that look up the value in the matrices if it is already there, and otherwise do recursion (or base cases) necessary to compute it. There are many ways to implement the recursive definition, and memoization is not the fastest of them.

Write a program test_align
Use your aligner module to write this master-slave alignment program. You should have options
--subst_matrix BLOSUM62
Giving the file or URL to read the substitution matrix from (integer values in BLOSUM format).
--align local and --align local
Specify either local or global alignment.
--open 10
Gap opening penalty.
--extend 2
Gap extension penalty.
--double_gap 3
Penalty for moving from gap in one sequence to gap in the other.

The program should read a fasta file from stdin with multiple sequences and output an alignment as an A2M file to stdout with the same sequences in the same order, and with the same names and comments on the id lines.

The fasta files may have alignment characters in them, like spaces, periods, and hyphens. These are not part of the sequences and should be removed or ignored. The alphabet for the fasta input is taken to be the characters used as indices in the substitution matrix. Whether letters in the fasta file are upper- or lower-case should also be ignored. This means that you should be able to read an A2M file and recover the unaligned sequences in a uniform case. As a sanity check, if you feed the output of your program back into the program with the same alignment options, the output should be identical to the input.

The first sequence of the fasta file is taken to be the "master" sequence and is echoed to stdout. All subsequent sequences are aligned to the first sequence and output to stdout. The score for each alignment should be output to stderr.

The gap penalties should (by default) be open=12, extend=1, double_gap=3. (The open and extend penalties are reasonable, but I'm not so sure about the double_gap penalty.)

The complete output of the program should be a legal A2M file.

>test1
CCCCCCCCCCCCCCCADEFGCCCCCCCCCCCCCCCC
>test2
eeeeeeeeee---------------ADEFG----------------eeeeeeeeeeeeeeeee
or
>test1
CCCCCCCCCCCCCCCADEFGCCCCCCCCCCCCCCCC
>test2
---------------eeeeeeeeeeADEFGeeeeeeeeeeeeeeeee----------------
See http://compbio.soe.ucsc.edu/a2m-desc.html for more details of a2m format. You can check that the file is formatted ok, by running it through the SAM program checkseq, which unfortunately does not use stdin. To check the file x.a2m, run
checkseq foo -d x.a2m
You can also use the SAM program prettyalign to convert the a2m into a more visually appealing format. The SAM programs are installed in /projects/compbio/bin/i686 and /projects/compbio/bin/x86_64

Note: Python is not the language of choice for computation-intensive applications like sequence alignment. It may be better than Perl, but it is still not very efficient. My memoized implementation on 1001 alignments to a 302-long sequence (input from 1bqcA.fa.gz) which requires about 98.2 million cell updates took 1380 seconds on my laptop, or about 71,000 cell-updates-per-second (CUPS). A reasonable software implementation in an efficient language should do about 30 million CUPS on a somewhat newer machine, and a highly optimized implementation taking advantage of all the hardware on fast modern desktops may manage 46 billion CUPS (http://www.clcbio.com/wp-content/uploads/2012/09/Cell_A413.pdf). Python is required for this assignment only for pedagogic purposes, and my implementation was chosen for simplicity, not speed. Even in Python, one could go much faster---the memoized one is ludicrously slow.

Name your program "test_align", so that I can test everyone's program with a script that does test_align --align=local < test.fasta and test_align --align=global < test.fasta . Turn in your program after testing it on the following fasta files:

(I may use other alignment test cases, if I can come up with some good ones for teasing out bugs in traceback.)

Don't waste paper turning in printouts of the alignments, as I'll be recomputing them with your program anyway. I'll be grading based on whether your alignment scores the same as an optimal alignment that my program finds, so it is essential that your output be a legal A2M file that correctly represents the alignment, but it need not be the same optimal alignment that I find.

To help you with your debugging, I provide possible output for local alignment of align-3, with g=12, e=1, d=2:

>fat
ISAYGARFIELDISAFATCAT
>cat
----GARFIELDwasthelASTFATCAT
>cat2
----GARFIELDISAFATCAT
>cat3
----GARFIELDISTheFATCAT
>cat4
----GARFIELDISAFAsTCAT
>cat4
----GARFIELDISTHefaSTCAT
>cat5 
----GARFIELDIS-FAT---
>say
-SAYFAt---------------
with scores 58, 85, 68, 73, 57, 51, and 16. Possible global alignments are
>fat
ISAYGARFIELDISAFATCAT
>cat
----GARFIELDwasthelASTFATCAT
>cat2
----GARFIELDISAFATCAT
>cat3
----GARFIELDISTheFATCAT
>cat4
----GARFIELDISAFAsTCAT
>cat4
----GARFIELDISTHefaSTCAT
>cat5
----GARFIELDIS----FAT
>say
-SAY--------------FAT
with scores 43, 70, 53, 58, 42, 25, and -15.

Bonus point—linear gap costs
Implement linear gap costs and compare running times for affine and linear algorithms. You may need to use longer test sequences to get repeatable timing results, as the expected 3-fold difference in the inner loop may be buried in overhead.

Things learned for next year

I was thinking of adding a requirement that score_a2m be symmetric in s1 and s2, not be assuming that one is a master sequence. Unfortunately, that requires figuring out how to score inserts or deletes that occur simultaneously in both sequences, which our simplified notion of affine gap costs does not cover. So perhaps I should make score_a2m explicitly allowing master/slave alignment scoring. Add another test case to align-3 that does not have a match in the final position for global alignment (currently only align-4 tests that). Add more test cases for local alignment that start in the middle of both sequences.

baskin-icon
SoE home
sketch of Kevin Karplus by Abe
Kevin Karplus's home page
BME-slug-icon
Biomolecular Engineering Department
BME 205 home page UCSC Bioinformatics research

Questions about page content should be directed to Kevin Karplus
Biomolecular Engineering
University of California, Santa Cruz
Santa Cruz, CA 95064
USA
karplus@soe.ucsc.edu
1-831-459-4250
318 Physical Sciences Building

Locations of visitors to pages with this footer (started 3 Nov 2008)