[Biococoa-dev] testing alignements

Philipp Seibel biococoa at bioworxx.com
Sat Mar 19 03:49:40 EST 2005


Sorry me writing so late, but i was in cologne for two days.
I will comment the code soon, i promise.
Your code doesn't work because you used an old version.... , there is a 
new version in the cvs since wednesday.
If the new piece of code doesn't work, please tell.

Phil

Am 17.03.2005 um 06:22 schrieb 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
> _______________________________________________
> Biococoa-dev mailing list
> Biococoa-dev at bioinformatics.org
> https://bioinformatics.org/mailman/listinfo/biococoa-dev
>
>




More information about the Biococoa-dev mailing list