Function alignOverlap

Align/smith_waterman.cpp:179–292

This is the Smith-Waterman algorithm (a variation of * Needleman-Wunsch algorithm), finds one optimal local alignment * Modified to find overlap of seq_a and seq_b (alignment that is * pinned at the end of seq_a and beginning of seq_b). * Actually, this should be a variation of needleman algorithm, that * looks for a global alignment, but without penalizing overhangs... * and make sure the a

 * of seqB).
 */
void alignOverlap(const string& seq_a, const string& seq_b, unsigned seq_a_start_pos,
	vector<overlap_align>& overlaps, bool multi_align, bool verbose)
{
	// get the actual lengths of the sequences
	int N_a = seq_a.length();
	int N_b = seq_b.length();

	// initialize H
	int i, j;
	double** H;
	int **I_i, **I_j;
	H = new double*[N_a+1];
	I_i = new int*[N_a+1];
	I_j = new int*[N_a+1];
	bool** V = new bool*[N_a+1];

	for(i=0;i<=N_a;i++){
		H[i] = new double[N_b+1];
		I_i[i] = new int[N_b+1];
		I_j[i] = new int[N_b+1];
		H[i][0]=0; //only need to initialize first row and first column
		I_i[i][0] = i-1;
		V[i] = new bool[N_b+1];
		V[i][0] = true; //valid start
	}

	for (j = 0; j <= N_b; j++) {
		H[0][j] = 0; //initialize first column
		I_j[0][j] = j-1;
		V[0][j] = false; //wrong start, not overlap
	}
	V[0][0] = true;

	for(i=1;i<=N_a;i++){
		for(j=1;j<=N_b;j++){
			char a = seq_a[i-1], b = seq_b[j-1];
			double scores[3] = {
				V[i-1][j-1] ? H[i-1][j-1] + matchScore(a, b)
					: -DBL_MAX, // match or mismatch
				V[i-1][j] ? H[i-1][j] + gapScore(I_j[i-1][j] == j)
					: -DBL_MAX, // deletion in sequence A
				V[i][j-1] ? H[i][j-1] + gapScore(I_i[i][j-1] == i)
					: -DBL_MAX // deletion in sequence B
			};
			double* pMax = max_element(scores, scores + 3);
			H[i][j] = *pMax;
			switch (pMax - scores) {
			  case 0: // match or mismatch
				I_i[i][j] = i-1;
				I_j[i][j] = j-1;
				break;
			  case 1: // deletion in sequence A
				I_i[i][j] = i-1;
				I_j[i][j] = j;
				break;
			  case 2: // deletion in sequence B
				I_i[i][j] = i;
				I_j[i][j] = j-1;