5 * Created by Pat Schloss on 12/15/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
8 * This class is an Alignment child class that implements the Gotoh pairwise alignment algorithm as described in:
10 * Gotoh O. 1982. An improved algorithm for matching biological sequences. J. Mol. Biol. 162:705-8.
11 * Myers, EW & Miller, W. 1988. Optimal alignments in linear space. Comput Appl Biosci. 4:11-7.
13 * This method is nice because it allows for an affine gap penalty to be assessed, which is analogous to what is used
14 * in blast and is an alternative to Needleman-Wunsch, which only charges the same penalty for each gap position.
15 * Because this method typically has problems at the ends when two sequences do not full overlap, we employ a separate
16 * method to fix the ends (see Overlap class documentation)
22 #include "alignmentcell.hpp"
23 #include "overlap.hpp"
24 #include "alignment.hpp"
25 #include "gotohoverlap.hpp"
27 /**************************************************************************************************/
29 GotohOverlap::GotohOverlap(float gO, float gE, float m, float mm, int r) :
30 gapOpen(gO), gapExtend(gE), match(m), mismatch(mm), Alignment(r) {
32 for(int i=1;i<nCols;i++){ // we initialize the dynamic programming matrix by setting the pointers in
33 alignment[0][i].prevCell = 'l'; // the first row to the left
34 alignment[0][i].cValue = 0;
35 alignment[0][i].dValue = 0;
38 for(int i=1;i<nRows;i++){ // we initialize the dynamic programming matrix by setting the pointers in
39 alignment[i][0].prevCell = 'u'; // the first column upward
40 alignment[i][0].cValue = 0;
41 alignment[i][0].iValue = 0;
45 /**************************************************************************************************/
47 void GotohOverlap::align(string A, string B){
49 seqA = ' ' + A; lA = seqA.length(); // the algorithm requires that the first character be a dummy value
50 seqB = ' ' + B; lB = seqB.length(); // the algorithm requires that the first character be a dummy value
52 for(int i=1;i<lB;i++){ // the recursion here is shown in Webb and Miller, Fig. 1A. Note that
53 for(int j=1;j<lA;j++){ // if we need to conserve on space we should see Fig. 1B, which is linear
54 // in space, which I think is unnecessary
56 if(seqB[i] == seqA[j]) { diagonal = alignment[i-1][j-1].cValue + match; }
57 else { diagonal = alignment[i-1][j-1].cValue + mismatch; }
59 alignment[i][j].iValue = max(alignment[i][j-1].iValue, alignment[i][j-1].cValue + gapOpen) + gapExtend;
60 alignment[i][j].dValue = max(alignment[i-1][j].dValue, alignment[i-1][j].cValue + gapOpen) + gapExtend;
62 if(alignment[i][j].iValue > alignment[i][j].dValue){
63 if(alignment[i][j].iValue > diagonal){
64 alignment[i][j].cValue = alignment[i][j].iValue;
65 alignment[i][j].prevCell = 'l';
68 alignment[i][j].cValue = diagonal;
69 alignment[i][j].prevCell = 'd';
73 if(alignment[i][j].dValue > diagonal){
74 alignment[i][j].cValue = alignment[i][j].dValue;
75 alignment[i][j].prevCell = 'u';
78 alignment[i][j].cValue = diagonal;
79 alignment[i][j].prevCell = 'd';
86 over.setOverlap(alignment, lA, lB, 0); // Fix the gaps at the ends of the sequences
87 traceBack(); // Construct the alignment and set seqAaln and seqBaln
90 /**************************************************************************************************/