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 "alignment.hpp"
24 #include "overlap.hpp"
25 #include "needlemanoverlap.hpp"
27 /**************************************************************************************************/
29 NeedlemanOverlap::NeedlemanOverlap(float gO, float m, float mm, int r) :// note that we don't have a gap extend
30 gap(gO), match(m), mismatch(mm), Alignment(r) { // the gap openning penalty is assessed for
31 // every gapped position
32 for(int i=1;i<nCols;i++){
33 alignment[0][i].prevCell = 'l'; // initialize first row by pointing all poiters to the left
34 alignment[0][i].cValue = 0; // and the score to zero
37 for(int i=1;i<nRows;i++){
38 alignment[i][0].prevCell = 'u'; // initialize first column by pointing all poiters upwards
39 alignment[i][0].cValue = 0; // and the score to zero
43 /**************************************************************************************************/
45 NeedlemanOverlap::~NeedlemanOverlap(){ /* do nothing */ }
47 /**************************************************************************************************/
49 void NeedlemanOverlap::align(string A, string B){
51 seqA = ' ' + A; lA = seqA.length(); // algorithm requires a dummy space at the beginning of each string
52 seqB = ' ' + B; lB = seqB.length(); // algorithm requires a dummy space at the beginning of each string
54 for(int i=1;i<lB;i++){ // This code was largely translated from Perl code provided in Ex 3.1
55 for(int j=1;j<lA;j++){ // of the O'Reilly BLAST book. I found that the example output had a
58 if(seqB[i] == seqA[j]) { diagonal = alignment[i-1][j-1].cValue + match; }
59 else { diagonal = alignment[i-1][j-1].cValue + mismatch; }
61 float up = alignment[i-1][j].cValue + gap;
62 float left = alignment[i][j-1].cValue + gap;
66 alignment[i][j].cValue = diagonal;
67 alignment[i][j].prevCell = 'd';
70 alignment[i][j].cValue = left;
71 alignment[i][j].prevCell = 'l';
76 alignment[i][j].cValue = up;
77 alignment[i][j].prevCell = 'u';
80 alignment[i][j].cValue = left;
81 alignment[i][j].prevCell = 'l';
87 over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences
88 traceBack(); // Traceback the alignment to populate seqAaln and seqBaln
91 /**************************************************************************************************/