# [Biococoa-dev] testing alignements

Charles PARNOT charles.parnot at stanford.edu
Thu Mar 17 00:22:03 EST 2005

```>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

```