1 #include "mltaln.h"
2 #include "dp.h"
3
4 #define DEBUG 0
5 #define DEBUG2 0
6 #define XXXXXXX 0
7 #define USE_PENALTY_EX 1
8
9 static TLS int localstop;
10
11 #if 1
match_calc_mtx(double ** mtx,double * match,char ** s1,char ** s2,int i1,int lgth2)12 static void match_calc_mtx( double **mtx, double *match, char **s1, char **s2, int i1, int lgth2 )
13 {
14 char *seq2 = s2[0];
15 double *doubleptr = mtx[(int)s1[0][i1]];
16
17 while( lgth2-- )
18 *match++ = doubleptr[(int)*seq2++];
19 }
20 #else
match_calc(double * match,char ** s1,char ** s2,int i1,int lgth2)21 static void match_calc( double *match, char **s1, char **s2, int i1, int lgth2 )
22 {
23 int j;
24
25 for( j=0; j<lgth2; j++ )
26 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
27 }
28 #endif
29
gentracking(double * lasthorizontalw,double * lastverticalw,char ** seq1,char ** seq2,char ** mseq1,char ** mseq2,double ** cpmx1,double ** cpmx2,int ** ijpi,int ** ijpj,int * off1pt,int * off2pt,int endi,int endj)30 static double gentracking( double *lasthorizontalw, double *lastverticalw,
31 char **seq1, char **seq2,
32 char **mseq1, char **mseq2,
33 double **cpmx1, double **cpmx2,
34 int **ijpi, int **ijpj, int *off1pt, int *off2pt, int endi, int endj )
35 {
36 int i, j, l, iin, jin, lgth1, lgth2, k, limk;
37 int ifi=0, jfi=0; // by D.Mathog
38 // char gap[] = "-";
39 char *gap;
40 gap = newgapstr;
41 lgth1 = strlen( seq1[0] );
42 lgth2 = strlen( seq2[0] );
43
44 #if 0
45 for( i=0; i<lgth1; i++ )
46 {
47 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
48 }
49 #endif
50
51 for( i=0; i<lgth1+1; i++ )
52 {
53 ijpi[i][0] = localstop;
54 ijpj[i][0] = localstop;
55 }
56 for( j=0; j<lgth2+1; j++ )
57 {
58 ijpi[0][j] = localstop;
59 ijpj[0][j] = localstop;
60 }
61
62 mseq1[0] += lgth1+lgth2;
63 *mseq1[0] = 0;
64 mseq2[0] += lgth1+lgth2;
65 *mseq2[0] = 0;
66 iin = endi; jin = endj;
67 limk = lgth1+lgth2;
68 for( k=0; k<=limk; k++ )
69 {
70
71 ifi = ( ijpi[iin][jin] );
72 jfi = ( ijpj[iin][jin] );
73 l = iin - ifi;
74 // if( ijpi[iin][jin] < 0 || ijpj[iin][jin] < 0 )
75 // {
76 // fprintf( stderr, "skip! %d-%d\n", ijpi[iin][jin], ijpj[iin][jin] );
77 // fprintf( stderr, "1: %c-%c\n", seq1[0][iin], seq1[0][ifi] );
78 // fprintf( stderr, "2: %c-%c\n", seq2[0][jin], seq2[0][jfi] );
79 // }
80 while( --l )
81 {
82 *--mseq1[0] = seq1[0][ifi+l];
83 *--mseq2[0] = *gap;
84 k++;
85 }
86 l= jin - jfi;
87 while( --l )
88 {
89 *--mseq1[0] = *gap;
90 *--mseq2[0] = seq2[0][jfi+l];
91 k++;
92 }
93
94 if( iin <= 0 || jin <= 0 ) break;
95 *--mseq1[0] = seq1[0][ifi];
96 *--mseq2[0] = seq2[0][jfi];
97
98 if( ijpi[ifi][jfi] == localstop ) break;
99 if( ijpj[ifi][jfi] == localstop ) break;
100 k++;
101 iin = ifi; jin = jfi;
102 }
103 if( ifi == -1 ) *off1pt = 0; else *off1pt = ifi;
104 if( jfi == -1 ) *off2pt = 0; else *off2pt = jfi;
105
106 // fprintf( stderr, "ifn = %d, jfn = %d\n", ifi, jfi );
107
108
109 return( 0.0 );
110 }
111
112
genL__align11(double ** n_dynamicmtx,char ** seq1,char ** seq2,int alloclen,int * off1pt,int * off2pt)113 double genL__align11( double **n_dynamicmtx, char **seq1, char **seq2, int alloclen, int *off1pt, int *off2pt )
114 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
115 {
116 // int k;
117 register int i, j;
118 int lasti, lastj;
119 int lgth1, lgth2;
120 int resultlen;
121 double wm = 0.0; /* int ?????? */
122 double g;
123 double *currentw, *previousw;
124 #if 1
125 double *wtmp;
126 int *ijpipt;
127 int *ijpjpt;
128 double *mjpt, *Mjpt, *prept, *curpt;
129 int *mpjpt, *Mpjpt;
130 #endif
131 static TLS double mi, *m;
132 static TLS double Mi, *largeM;
133 static TLS int **ijpi;
134 static TLS int **ijpj;
135 static TLS int mpi, *mp;
136 static TLS int Mpi, *Mp;
137 static TLS double *w1, *w2;
138 static TLS double *match;
139 static TLS double *initverticalw; /* kufuu sureba iranai */
140 static TLS double *lastverticalw; /* kufuu sureba iranai */
141 static TLS char **mseq1;
142 static TLS char **mseq2;
143 static TLS char **mseq;
144 static TLS double **cpmx1;
145 static TLS double **cpmx2;
146 static TLS int **intwork;
147 static TLS double **doublework;
148 static TLS int orlgth1 = 0, orlgth2 = 0;
149 static TLS double **amino_dynamicmtx = NULL; // ??
150 double maxwm;
151 double tbk;
152 int tbki, tbkj;
153 int endali, endalj;
154 // double localthr = 0.0;
155 // double localthr2 = 0.0;
156 double fpenalty = (double)penalty;
157 double fpenalty_OP = (double)penalty_OP;
158 double fpenalty_ex = (double)penalty_ex;
159 // double fpenalty_EX = (double)penalty_EX;
160 double foffset = (double)offset;
161 double localthr = -foffset;
162 double localthr2 = -foffset;
163
164 if( seq1 == NULL )
165 {
166 if( orlgth1 > 0 && orlgth2 > 0 )
167 {
168 orlgth1 = 0;
169 orlgth2 = 0;
170 free( mseq1 );
171 free( mseq2 );
172 FreeFloatVec( w1 );
173 FreeFloatVec( w2 );
174 FreeFloatVec( match );
175 FreeFloatVec( initverticalw );
176 FreeFloatVec( lastverticalw );
177
178 FreeFloatVec( m );
179 FreeIntVec( mp );
180 free( largeM );
181 free( Mp );
182
183 FreeCharMtx( mseq );
184
185 FreeFloatMtx( cpmx1 );
186 FreeFloatMtx( cpmx2 );
187
188 FreeFloatMtx( doublework );
189 FreeIntMtx( intwork );
190 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
191
192 }
193 return( 0.0 );
194 }
195
196
197
198 // fprintf( stderr, "@@@@@@@@@@@@@ penalty_OP = %f, penalty_EX = %f, pelanty = %f\n", fpenalty_OP, fpenalty_EX, fpenalty );
199
200 if( orlgth1 == 0 )
201 {
202 mseq1 = AllocateCharMtx( njob, 0 );
203 mseq2 = AllocateCharMtx( njob, 0 );
204 }
205
206
207 lgth1 = strlen( seq1[0] );
208 lgth2 = strlen( seq2[0] );
209
210 if( lgth1 > orlgth1 || lgth2 > orlgth2 )
211 {
212 int ll1, ll2;
213
214 if( orlgth1 > 0 && orlgth2 > 0 )
215 {
216 FreeFloatVec( w1 );
217 FreeFloatVec( w2 );
218 FreeFloatVec( match );
219 FreeFloatVec( initverticalw );
220 FreeFloatVec( lastverticalw );
221
222 FreeFloatVec( m );
223 FreeIntVec( mp );
224 FreeFloatVec( largeM );
225 FreeIntVec( Mp );
226
227 FreeCharMtx( mseq );
228
229 FreeFloatMtx( cpmx1 );
230 FreeFloatMtx( cpmx2 );
231
232 FreeFloatMtx( doublework );
233 FreeIntMtx( intwork );
234 if( amino_dynamicmtx ) FreeDoubleMtx( amino_dynamicmtx ); amino_dynamicmtx = NULL;
235 }
236
237 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
238 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
239
240 #if DEBUG
241 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
242 #endif
243
244 w1 = AllocateFloatVec( ll2+2 );
245 w2 = AllocateFloatVec( ll2+2 );
246 match = AllocateFloatVec( ll2+2 );
247
248 initverticalw = AllocateFloatVec( ll1+2 );
249 lastverticalw = AllocateFloatVec( ll1+2 );
250
251 m = AllocateFloatVec( ll2+2 );
252 mp = AllocateIntVec( ll2+2 );
253 largeM = AllocateFloatVec( ll2+2 );
254 Mp = AllocateIntVec( ll2+2 );
255
256 mseq = AllocateCharMtx( njob, ll1+ll2 );
257
258 cpmx1 = AllocateFloatMtx( nalphabets, ll1+2 );
259 cpmx2 = AllocateFloatMtx( nalphabets, ll2+2 );
260
261 doublework = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 );
262 intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 );
263
264 amino_dynamicmtx = AllocateDoubleMtx( 0x80, 0x80 );
265
266 #if DEBUG
267 fprintf( stderr, "succeeded\n" );
268 #endif
269
270 orlgth1 = ll1 - 100;
271 orlgth2 = ll2 - 100;
272 }
273
274 for( i=0; i<nalphabets; i++) for( j=0; j<nalphabets; j++ )
275 amino_dynamicmtx[(int)amino[i]][(int)amino[j]] = (double)n_dynamicmtx[i][j];
276
277 mseq1[0] = mseq[0];
278 mseq2[0] = mseq[1];
279
280
281 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
282 {
283 int ll1, ll2;
284
285 if( commonAlloc1 && commonAlloc2 )
286 {
287 FreeIntMtx( commonIP );
288 FreeIntMtx( commonJP );
289 }
290
291 ll1 = MAX( orlgth1, commonAlloc1 );
292 ll2 = MAX( orlgth2, commonAlloc2 );
293
294 #if DEBUG
295 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
296 #endif
297
298 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
299 commonJP = AllocateIntMtx( ll1+10, ll2+10 );
300
301 #if DEBUG
302 fprintf( stderr, "succeeded\n\n" );
303 #endif
304
305 commonAlloc1 = ll1;
306 commonAlloc2 = ll2;
307 }
308 ijpi = commonIP;
309 ijpj = commonJP;
310
311
312 #if 0
313 for( i=0; i<lgth1; i++ )
314 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
315 #endif
316
317 currentw = w1;
318 previousw = w2;
319
320 match_calc_mtx( amino_dynamicmtx, initverticalw, seq2, seq1, 0, lgth1 );
321
322 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, 0, lgth2 );
323
324
325 lasti = lgth2+1;
326 for( j=1; j<lasti; ++j )
327 {
328 m[j] = currentw[j-1]; mp[j] = 0;
329 largeM[j] = currentw[j-1]; Mp[j] = 0;
330 }
331
332 lastverticalw[0] = currentw[lgth2-1];
333
334
335 #if 0
336 fprintf( stderr, "currentw = \n" );
337 for( i=0; i<lgth1+1; i++ )
338 {
339 fprintf( stderr, "%5.2f ", currentw[i] );
340 }
341 fprintf( stderr, "\n" );
342 fprintf( stderr, "initverticalw = \n" );
343 for( i=0; i<lgth2+1; i++ )
344 {
345 fprintf( stderr, "%5.2f ", initverticalw[i] );
346 }
347 fprintf( stderr, "\n" );
348 #endif
349 #if DEBUG2
350 fprintf( stderr, "\n" );
351 fprintf( stderr, " " );
352 for( j=0; j<lgth2+1; j++ )
353 fprintf( stderr, "%c ", seq2[0][j] );
354 fprintf( stderr, "\n" );
355 #endif
356
357 localstop = lgth1+lgth2+1;
358 maxwm = -999999999.9;
359 endali = endalj = 0;
360 #if DEBUG2
361 fprintf( stderr, "\n" );
362 fprintf( stderr, "%c ", seq1[0][0] );
363
364 for( j=0; j<lgth2+1; j++ )
365 fprintf( stderr, "%5.0f ", currentw[j] );
366 fprintf( stderr, "\n" );
367 #endif
368
369 lasti = lgth1+1;
370 for( i=1; i<lasti; i++ )
371 {
372 wtmp = previousw;
373 previousw = currentw;
374 currentw = wtmp;
375
376 previousw[0] = initverticalw[i-1];
377
378 match_calc_mtx( amino_dynamicmtx, currentw, seq1, seq2, i, lgth2 );
379 #if DEBUG2
380 fprintf( stderr, "%c ", seq1[0][i] );
381 fprintf( stderr, "%5.0f ", currentw[0] );
382 #endif
383
384 #if XXXXXXX
385 fprintf( stderr, "\n" );
386 fprintf( stderr, "i=%d\n", i );
387 fprintf( stderr, "currentw = \n" );
388 for( j=0; j<lgth2; j++ )
389 {
390 fprintf( stderr, "%5.2f ", currentw[j] );
391 }
392 fprintf( stderr, "\n" );
393 #endif
394 #if XXXXXXX
395 fprintf( stderr, "\n" );
396 fprintf( stderr, "i=%d\n", i );
397 fprintf( stderr, "currentw = \n" );
398 for( j=0; j<lgth2; j++ )
399 {
400 fprintf( stderr, "%5.2f ", currentw[j] );
401 }
402 fprintf( stderr, "\n" );
403 #endif
404 currentw[0] = initverticalw[i];
405
406 mi = previousw[0]; mpi = 0;
407 Mi = previousw[0]; Mpi = 0;
408
409 #if 0
410 if( mi < localthr ) mi = localthr2;
411 #endif
412
413 ijpipt = ijpi[i] + 1;
414 ijpjpt = ijpj[i] + 1;
415 mjpt = m + 1;
416 Mjpt = largeM + 1;
417 prept = previousw;
418 curpt = currentw + 1;
419 mpjpt = mp + 1;
420 Mpjpt = Mp + 1;
421 tbk = -999999.9;
422 tbki = 0;
423 tbkj = 0;
424 lastj = lgth2+1;
425 for( j=1; j<lastj; j++ )
426 {
427 wm = *prept;
428 *ijpipt = i-1;
429 *ijpjpt = j-1;
430
431
432 // fprintf( stderr, "i,j=%d,%d %c-%c\n", i, j, seq1[0][i], seq2[0][j] );
433 // fprintf( stderr, "wm=%f\n", wm );
434 #if 0
435 fprintf( stderr, "%5.0f->", wm );
436 #endif
437 g = mi + fpenalty;
438 #if 0
439 fprintf( stderr, "%5.0f?", g );
440 #endif
441 if( g > wm )
442 {
443 wm = g;
444 // *ijpipt = i - 1;
445 *ijpjpt = mpi;
446 }
447 g = *prept;
448 if( g > mi )
449 {
450 mi = g;
451 mpi = j-1;
452 }
453
454 #if USE_PENALTY_EX
455 mi += fpenalty_ex;
456 #endif
457
458 #if 0
459 fprintf( stderr, "%5.0f->", wm );
460 #endif
461 g = *mjpt + fpenalty;
462 #if 0
463 fprintf( stderr, "m%5.0f?", g );
464 #endif
465 if( g > wm )
466 {
467 wm = g;
468 *ijpipt = *mpjpt;
469 *ijpjpt = j - 1; //IRU!
470 }
471 g = *prept;
472 if( g > *mjpt )
473 {
474 *mjpt = g;
475 *mpjpt = i-1;
476 }
477 #if USE_PENALTY_EX
478 *mjpt += fpenalty_ex;
479 #endif
480
481
482 g = tbk + fpenalty_OP;
483 // g = tbk;
484 if( g > wm )
485 {
486 wm = g;
487 *ijpipt = tbki;
488 *ijpjpt = tbkj;
489 // fprintf( stderr, "hit! i%d, j%d, ijpi = %d, ijpj = %d\n", i, j, *ijpipt, *ijpjpt );
490 }
491 // g = Mi;
492 if( Mi > tbk )
493 {
494 tbk = Mi; //error desu.
495 tbki = i-1;
496 tbkj = Mpi;
497 }
498 // g = *Mjpt;
499 if( *Mjpt > tbk )
500 {
501 tbk = *Mjpt;
502 tbki = *Mpjpt;
503 tbkj = j-1;
504 }
505 // tbk += fpenalty_EX;// + foffset;
506
507 // g = *prept;
508 if( *prept > *Mjpt )
509 {
510 *Mjpt = *prept;
511 *Mpjpt = i-1;
512 }
513 // *Mjpt += fpenalty_EX;// + foffset;
514
515 // g = *prept;
516 if( *prept > Mi )
517 {
518 Mi = *prept;
519 Mpi = j-1;
520 }
521 // Mi += fpenalty_EX;// + foffset;
522
523
524 // fprintf( stderr, "wm=%f, tbk=%f(%c-%c), mi=%f, *mjpt=%f\n", wm, tbk, seq1[0][tbki], seq2[0][tbkj], mi, *mjpt );
525 // fprintf( stderr, "ijp = %c,%c\n", seq1[0][abs(*ijpipt)], seq2[0][abs(*ijpjpt)] );
526
527
528 if( maxwm < wm )
529 {
530 maxwm = wm;
531 endali = i;
532 endalj = j;
533 }
534 #if 1
535 if( wm < localthr )
536 {
537 // fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
538 *ijpipt = localstop;
539 // *ijpjpt = localstop;
540 wm = localthr2;
541 }
542 #endif
543 #if 0
544 fprintf( stderr, "%5.0f ", *curpt );
545 #endif
546 #if DEBUG2
547 fprintf( stderr, "%5.0f ", wm );
548 // fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
549 #endif
550
551 *curpt += wm;
552 ijpipt++;
553 ijpjpt++;
554 mjpt++;
555 Mjpt++;
556 prept++;
557 mpjpt++;
558 Mpjpt++;
559 curpt++;
560 }
561 #if DEBUG2
562 fprintf( stderr, "\n" );
563 #endif
564
565 lastverticalw[i] = currentw[lgth2-1];
566 }
567
568
569 #if DEBUG2
570 fprintf( stderr, "maxwm = %f\n", maxwm );
571 fprintf( stderr, "endali = %d\n", endali );
572 fprintf( stderr, "endalj = %d\n", endalj );
573 #endif
574
575 if( ijpi[endali][endalj] == localstop ) // && ijpj[endali][endalj] == localstop )
576 {
577 strcpy( seq1[0], "" );
578 strcpy( seq2[0], "" );
579 *off1pt = *off2pt = 0;
580 return( 0.0 );
581 }
582
583
584 gentracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijpi, ijpj, off1pt, off2pt, endali, endalj );
585
586 // fprintf( stderr, "### impmatch = %f\n", *impmatch );
587
588 resultlen = strlen( mseq1[0] );
589 if( alloclen < resultlen || resultlen > N )
590 {
591 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
592 ErrorExit( "LENGTH OVER!\n" );
593 }
594
595
596 strcpy( seq1[0], mseq1[0] );
597 strcpy( seq2[0], mseq2[0] );
598
599 #if 0
600 fprintf( stderr, "\n" );
601 fprintf( stderr, ">\n%s\n", mseq1[0] );
602 fprintf( stderr, ">\n%s\n", mseq2[0] );
603 #endif
604
605
606 return( maxwm );
607 }
608
609