3 void profilealignment( int n0, int n1, int n2, char **aln0, char **aln1, char **aln2, int alloclen, char alg ) // n1 ha allgap
6 double *effarr0, *effarr2;
9 effarr0 = AllocateDoubleVec( n0 );
10 effarr2 = AllocateDoubleVec( n2 );
12 commongappick( n0, aln0 );
13 commongappick( n2, aln2 );
16 eff = 1.0 / (double)n0; for( i=0; i<n0; i++ ) effarr0[i] = eff;
17 eff = 1.0 / (double)n2; for( i=0; i<n2; i++ ) effarr2[i] = eff;
20 MSalignmm( aln0, aln2, effarr0, effarr2, n0, n2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); //outgap=1??
22 A__align( aln0, aln2, effarr0, effarr2, n0, n2, alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); //outgap=1??
24 newlen = strlen( aln0[0] );
26 for( i=0; i<newlen; i++ ) aln1[0][i] = '-';
28 for( i=1; i<n1; i++ ) strcpy( aln1[i], aln1[0] );
35 void eq2dash( char *s )
39 if( *s == '=' ) *s = '-';
44 void findnewgaps( int n, char **seq, int *gaplen )
47 len = strlen( seq[0] );
49 // fprintf( stderr, "seq[0] = %s\n", seq[0] );
50 for( i=0; i<len; i++ ) gaplen[i] = 0;
53 for( i=0; i<len; i++ )
55 if( seq[0][i] == '=' )
57 // fprintf( stderr, "Newgap! pos = %d\n", pos );
65 void findcommongaps( int n, char **seq, int *gaplen )
67 int i, j, pos, len, len1;
68 len = strlen( seq[0] );
71 // fprintf( stderr, "seq[0] = %s\n", seq[0] );
72 for( i=0; i<len1; i++ ) gaplen[i] = 0;
75 for( i=0; i<len; i++ )
78 if( seq[j][i] != '-' ) break;
80 if( j == n ) gaplen[pos]++;
85 for( i=0; i<pos; i++ )
87 fprintf( stderr, "vec[%d] = %d\n", i, gaplen[i] );
92 void adjustgapmap( int newlen, int *gapmap, char *seq )
96 int newlen1 = newlen+1;
99 tmpmap = AllocateIntVec( newlen+2 );
104 // fprintf( stderr, "j=%d *seq = %c\n", j, *seq );
109 tmpmap[j++] = gapmap[pos++];
112 tmpmap[j++] = gapmap[pos];
114 for(j=0; j<newlen1; j++)
115 gapmap[j] = tmpmap[j];
120 void insertnewgaps( int njob, int *alreadyaligned, char **seq, int *ex1, int *ex2, int *gaplen, int *gapmap, int alloclen, char alg )
125 int i, j, k, len, rep, len0;
126 char **mseq2, **mseq0, **mseq1;
128 int ngroup2, ngroup0, ngroup1;
129 int *list0, *list1, *list2;
130 int posin12, gapshift, newpos;
132 mar = calloc( njob, sizeof( int ) );
133 list0 = calloc( njob, sizeof( int ) );
134 list1 = calloc( njob, sizeof( int ) );
135 list2 = calloc( njob, sizeof( int ) );
137 for( i=0; i<njob; i++ ) mar[i] = 0;
138 for( i=0; i<njob; i++ )
140 if( alreadyaligned[i]==0 ) mar[i] = 3;
142 for( i=0; (k=ex1[i])>-1; i++ )
145 // fprintf( stderr, "excluding %d\n", ex1[i] );
147 for( i=0; (k=ex2[i])>-1; i++ )
150 // fprintf( stderr, "excluding %d\n", ex2[i] );
153 ngroup2 = ngroup1 = ngroup0 = 0;
154 for( i=0; i<njob; i++ )
172 list0[ngroup0] = list1[ngroup1] = list2[ngroup2] = -1;
175 fprintf( stderr, "Nothing to do\n" );
183 for( i=0; i<njob; i++ ) if( mar[i] == 0 ) break;
185 len = strlen( seq[rep] );
191 //// fprintf( stderr, "Nothing to do\n" );
196 mseq2 = AllocateCharMtx( ngroup2, alloclen );
197 mseq1 = AllocateCharMtx( ngroup1, alloclen );
198 mseq0 = AllocateCharMtx( ngroup0, alloclen );
199 aseq = AllocateCharMtx( njob, alloclen );
200 gaps = calloc( alloclen, sizeof( char ) );
202 for( i=0; i<njob; i++ ) aseq[i][0] = 0;
204 for( j=0; j<len0; j++ )
208 for( i=0; i<ngroup0; i++ ) mseq0[i][0] = 0;
209 for( i=0; i<ngroup1; i++ ) mseq1[i][0] = 0;
210 for( i=0; i<ngroup2; i++ ) mseq2[i][0] = 0;
212 gapshift = gaplen[j];
214 while( gapshift-- ) *cptr++ = '-';
216 gapshift = gaplen[j];
218 for( i=0; i<ngroup0; i++ ) strncat( mseq0[i], gaps, gapshift );
219 for( i=0; i<ngroup1; i++ ) strncat( mseq1[i], seq[list1[i]]+posin12, gapshift );
220 for( i=0; i<ngroup2; i++ ) strncat( mseq2[i], seq[list2[i]]+posin12, gapshift );
223 gapshift = gapmap[posin12];
224 // fprintf( stderr, "gapmap[%d] kouho = %d\n", posin12, gapmap[posin12] );
227 for( i=0; i<ngroup0; i++ ) strncat( mseq0[i], seq[list0[i]]+j, gapshift );
228 for( i=0; i<ngroup1; i++ ) strncat( mseq1[i], seq[list1[i]]+posin12, gapshift );
229 for( i=0; i<ngroup2; i++ ) strncat( mseq2[i], seq[list2[i]]+posin12, gapshift );
231 for( i=0; i<ngroup0; i++ ) fprintf( stderr, "### mseq0[%d] = %s\n", i, mseq0[i] );
232 for( i=0; i<ngroup1; i++ ) fprintf( stderr, "### mseq1[%d] = %s\n", i, mseq1[i] );
233 for( i=0; i<ngroup2; i++ ) fprintf( stderr, "### mseq2[%d] = %s\n", i, mseq2[i] );
236 if( gapshift ) profilealignment( ngroup0, ngroup1, ngroup2, mseq0, mseq1, mseq2, alloclen, alg );
241 for( i=0; i<ngroup0; i++ ) strcat( aseq[list0[i]], mseq0[i] );
242 for( i=0; i<ngroup1; i++ ) strcat( aseq[list1[i]], mseq1[i] );
243 for( i=0; i<ngroup2; i++ ) strcat( aseq[list2[i]], mseq2[i] );
246 newpos = strlen( aseq[rep] );
247 for( i=0; i<ngroup0; i++ ) aseq[list0[i]][newpos] = seq[list0[i]][j];
248 for( i=0; i<ngroup1; i++ ) aseq[list1[i]][newpos] = seq[list1[i]][posin12];
249 for( i=0; i<ngroup2; i++ ) aseq[list2[i]][newpos] = seq[list2[i]][posin12];
250 for( i=0; i<ngroup0; i++ ) aseq[list0[i]][newpos+1] = 0;
251 for( i=0; i<ngroup1; i++ ) aseq[list1[i]][newpos+1] = 0;
252 for( i=0; i<ngroup2; i++ ) aseq[list2[i]][newpos+1] = 0;
257 // for( i=0; i<njob; i++ ) if( mar[i] != 3 ) strcpy( seq[i], aseq[i] );
258 for( i=0; i<ngroup0; i++ ) strcpy( seq[list0[i]], aseq[list0[i]] );
259 for( i=0; i<ngroup1; i++ ) strcpy( seq[list1[i]], aseq[list1[i]] );
260 for( i=0; i<ngroup2; i++ ) strcpy( seq[list2[i]], aseq[list2[i]] );
267 FreeCharMtx( mseq2 );
268 FreeCharMtx( mseq0 );
272 void restorecommongaps( int njob, char **seq, int *ex1, int *ex2, int *gaplen, int alloclen )
279 int i, j, k, len, rep, len1;
281 mar = calloc( njob, sizeof( int ) );
282 tmpseq = calloc( alloclen, sizeof( char ) );
283 tmpgaplen = calloc( alloclen, sizeof( int ) );
284 // tmpseq = calloc( alloclen+2, sizeof( char ) );
285 // tmpgaplen = calloc( alloclen+2, sizeof( int ) );
288 for( i=0; i<njob; i++ ) mar[i] = 0;
289 for( i=0; (k=ex1[i])>-1; i++ )
292 // fprintf( stderr, "excluding %d\n", ex1[i] );
294 for( i=0; (k=ex2[i])>-1; i++ )
297 // fprintf( stderr, "excluding %d\n", ex2[i] );
300 for( i=0; i<njob; i++ )
305 // fprintf( stderr, "Nothing to do\n" );
312 len = strlen( seq[rep] );
316 for( i=0; i<njob; i++ )
318 if( mar[i] == 0 ) continue;
320 for( j=0; j<len1; j++ )
322 for( k=0; k<gaplen[j]; k++ )
324 *(cptr++) = seq[i][j];
327 strcpy( seq[i], tmpseq );
331 for( j=0; j<len1; j++ )
333 *(iptr++) = gaplen[j];
334 for( k=0; k<gaplen[j]; k++ )
340 while( *iptr != -1 ) *gaplen++ = *iptr++;
348 int samemember( int *mem, int *cand )
353 fprintf( stderr, "mem = " );
354 for( i=0; mem[i]>-1; i++ ) fprintf( stderr, "%d ", mem[i] );
355 fprintf( stderr, "\n" );
357 fprintf( stderr, "cand = " );
358 for( i=0; cand[i]>-1; i++ ) fprintf( stderr, "%d ", cand[i] );
359 fprintf( stderr, "\n" );
362 for( i=0, j=0; mem[i]>-1; )
364 if( mem[i++] != cand[j++] ) return( 0 );
377 int samemember( int *mem, int *cand )
382 fprintf( stderr, "mem = " );
383 for( i=0; mem[i]>-1; i++ ) fprintf( stderr, "%d ", mem[i] );
384 fprintf( stderr, "\n" );
386 fprintf( stderr, "cand = " );
387 for( i=0; cand[i]>-1; i++ ) fprintf( stderr, "%d ", cand[i] );
388 fprintf( stderr, "\n" );
391 nm = 0; for( i=0; mem[i]>-1; i++ ) nm++;
392 nc = 0; for( i=0; cand[i]>-1; i++ ) nc++;
394 if( nm != nc ) return( 0 );
396 for( i=0; mem[i]>-1; i++ )
398 for( j=0; cand[j]>-1; j++ )
399 if( mem[i] == cand[j] ) break;
400 if( cand[j] == -1 ) return( 0 );
415 int includemember( int *mem, int *cand ) // mem in cand
420 fprintf( stderr, "mem = " );
421 for( i=0; mem[i]>-1; i++ ) fprintf( stderr, "%d ", mem[i] );
422 fprintf( stderr, "\n" );
424 fprintf( stderr, "cand = " );
425 for( i=0; cand[i]>-1; i++ ) fprintf( stderr, "%d ", cand[i] );
426 fprintf( stderr, "\n" );
429 for( i=0; mem[i]>-1; i++ )
431 for( j=0; cand[j]>-1; j++ )
432 if( mem[i] == cand[j] ) break;
433 if( cand[j] == -1 ) return( 0 );
435 // fprintf( stderr, "INCLUDED! mem[0]=%d\n", mem[0] );