[Biococoa-dev] testing alignements

[Charles PARNOT]

>I actually started writing the program by copying and pasting the code, but the alignement does not work. I will post the code on a separate message for Phil and you to have a look, because I have no idea how the algorithm works (Koen, you are not alone!).

Here is the code I used for testing purpose, but the alignement does not work, it seems


/*****************STARTING HERE***************/

/*
 *  TestAlignment.c
 *
 */

#define DIAG (idxB - 1) * lenA + (idxA - 1)
#define LEFT idxB * lenA + (idxA - 1)
#define UP  (idxB -1) * lenA + idxA

typedef enum {
    kNone = 0,
	kDiagonal,
    kLeft,
    kUp
} Pointers;


void alignSequences(char *seqA, char *seqB, int lenA, int lenB)
{

	/* set up score matrix */
	int *scoreMatrix = (int *)malloc( sizeof( int ) * 128 * 128 );
	int match = 1;
	int mismatch = -1;
	int ii,jj;
	for ( ii = 0; ii < 128 ; ii++ )
		for ( jj = 0; jj < 128 ; jj++ )
			if (ii == jj )
				scoreMatrix[ii+128*jj] = match;
			else
				scoreMatrix[ii+128*jj] = mismatch;
	
	int gapCosts = -1;
	
	int *backtracking = (int *)malloc( sizeof( int ) * lenA * lenB );
	int *dynMatrix    = (int *)malloc( sizeof( int ) * lenA * lenB );
	
	unsigned int idxA;
	unsigned int idxB;
	
	dynMatrix[ 0 ] = scoreMatrix[seqA[0]*128+seqB[0]];
	backtracking[ 0 ] = kNone;
	
	for ( idxA = 1; idxA < lenA; idxA++ ) {
	    backtracking[ idxA ] = kLeft;
		dynMatrix[ idxA ] = idxA * gapCosts;
	}
	
	for ( idxB = 1; idxB < lenB; idxB++ ) {
		backtracking[ idxB * lenA ] = kUp;
		dynMatrix[ idxB * lenA ] = idxB * gapCosts;
	}
	
	for ( idxA = 1; idxA < lenA; idxA++ ) {
		for ( idxB = 1; idxB < lenB; idxB++ ) {
			unsigned int currPos = idxB * lenA + idxA;
			
			int substitutionScore = scoreMatrix[seqA[idxA]*128+seqB[idxB]];
			int diagScore = dynMatrix[ DIAG ] + substitutionScore;
			int rightScore = dynMatrix[ LEFT ] + gapCosts;
			int downScore = dynMatrix[ UP ] + gapCosts;
			
			if ( diagScore >= rightScore ) {
				if ( diagScore > downScore ) {
				    backtracking[ currPos ] = kDiagonal;
					dynMatrix[ currPos ] = diagScore;
				}
				else {
					backtracking[ currPos ] = kUp;
					dynMatrix[ currPos ] = downScore;
				}
			}
			else {
				if ( rightScore > downScore ) {
					backtracking[ currPos ] = kLeft;
					dynMatrix[ currPos ] = rightScore;
				}
				else {
					backtracking[ currPos ] = kUp;
					dynMatrix[ currPos ] = downScore;
				}
			}
		}
	}
	
	int i = lenA;
	int j = lenB;
	int k = 0;
	char *a = ( char * ) malloc( (lenA + lenB) * sizeof(char));
	char *b = ( char * ) malloc( (lenA + lenB) * sizeof(char));
	
	while ( 1 ) {
		// escape when origin is reached
		if(backtracking[i * lenA + j ] == kNone) break;
		
		switch(backtracking[i * lenA + j ]){
			case kDiagonal :
				a[k] = seqA[i - 1];
				b[k] = seqB[j - 1];
				i--;
				j--;
				k++;
				break;
				
			case kLeft :
				a[k] = seqA[i - 1];
				b[k] = '-';
				i--;
				k++;
				break;
				
			case kUp :
				a[k] = '-';
				b[k] = seqB[j - 1];
				j--;
				k++;
				break;
		}
	}

	for(i=k-1;i>=0;i--) printf("%c",a[i]);
	printf("\n");
	for(j=k-1;j>=0;j--) printf("%c",b[j]);
	printf("\n");
	
}

void alignment1()
{
	char *seqA = "ATGTAGTCTGATGATGAGATGACGT";
	char *seqB = "ATGTCAGTCTGATGATGAGATGACGAT";
	alignSequences(seqA,seqB,25,27);
}


int main (int argc, const char * argv[]) {
	while(1<2) {
		alignment1();
	}
    return 0;
}

/*****************ENDING HERE***************/

-- 
Help science go fast forward:
http://cmgm.stanford.edu/~cparnot/xgrid-stanford/

Charles Parnot
charles.parnot at stanford.edu

Room  B157 in Beckman Center
279, Campus Drive
Stanford University
Stanford, CA 94305 (USA)

Tel +1 650 725 7754
Fax +1 650 725 8021