7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
14 //Values as provided in Probcons V1.1
15 static float EXP_UNDERFLOW_THRESHOLD = -4.60f;
16 static float LOG_UNDERFLOW_THRESHOLD = 7.50f;
17 //static float LOG_ZERO = -FLT_MAX;
18 static float LOG_ZERO=-200000004008175468544.000000;
19 static float LOG_ONE = 0.0f;
21 //Probabilistic Alignment DNA
26 static float DNAinitDistrib2Default[] = { 0.9615409374f, 0.0000004538f, 0.0000004538f, 0.0192291681f, 0.0192291681f };
27 static float DNAgapOpen2Default[] = { 0.0082473317f, 0.0082473317f, 0.0107844425f, 0.0107844425f };
28 static float DNAgapExtend2Default[] = { 0.3210460842f, 0.3210460842f, 0.3298229277f, 0.3298229277f };
30 static char DNAalphabetDefault[] = "ACGUTN";
31 static float DNAemitSingleDefault[6] = {0.2270790040f, 0.2422080040f, 0.2839320004f, 0.2464679927f, 0.2464679927f, 0.0003124650f};
33 static float DNAemitPairsDefault[6][6] = {
34 { 0.1487240046f, 0.0184142999f, 0.0361397006f, 0.0238473993f, 0.0238473993f, 0.0000375308f },
35 { 0.0184142999f, 0.1583919972f, 0.0275536999f, 0.0389291011f, 0.0389291011f, 0.0000815823f },
36 { 0.0361397006f, 0.0275536999f, 0.1979320049f, 0.0244289003f, 0.0244289003f, 0.0000824765f },
37 { 0.0238473993f, 0.0389291011f, 0.0244289003f, 0.1557479948f, 0.1557479948f, 0.0000743985f },
38 { 0.0238473993f, 0.0389291011f, 0.0244289003f, 0.1557479948f, 0.1557479948f, 0.0000743985f },
39 { 0.0000375308f, 0.0000815823f, 0.0000824765f, 0.0000743985f, 0.0000743985f, 0.0000263252f }
42 //Probabilistic ALignment Blosum62mt
47 static float initDistrib2Default[] = { 0.6814756989f, 8.615339902e-05f, 8.615339902e-05f, 0.1591759622f, 0.1591759622 };
48 static float gapOpen2Default[] = { 0.0119511066f, 0.0119511066f, 0.008008334786f, 0.008008334786 };
49 static float gapExtend2Default[] = { 0.3965826333f, 0.3965826333f, 0.8988758326f, 0.8988758326 };
51 static char alphabetDefault[] = "ARNDCQEGHILKMFPSTWYV";
52 static float emitSingleDefault[20] = {
53 0.07831005f, 0.05246024f, 0.04433257f, 0.05130349f, 0.02189704f,
54 0.03585766f, 0.05615771f, 0.07783433f, 0.02601093f, 0.06511648f,
55 0.09716489f, 0.05877077f, 0.02438117f, 0.04463228f, 0.03940142f,
56 0.05849916f, 0.05115306f, 0.01203523f, 0.03124726f, 0.07343426f
59 static float emitPairsDefault[20][20] = {
60 {0.02373072f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
61 {0.00244502f, 0.01775118f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
62 {0.00210228f, 0.00207782f, 0.01281864f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
63 {0.00223549f, 0.00161657f, 0.00353540f, 0.01911178f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
64 {0.00145515f, 0.00044701f, 0.00042479f, 0.00036798f, 0.01013470f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
65 {0.00219102f, 0.00253532f, 0.00158223f, 0.00176784f, 0.00032102f, 0.00756604f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
66 {0.00332218f, 0.00268865f, 0.00224738f, 0.00496800f, 0.00037956f, 0.00345128f, 0.01676565f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
67 {0.00597898f, 0.00194865f, 0.00288882f, 0.00235249f, 0.00071206f, 0.00142432f, 0.00214860f, 0.04062876f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
68 {0.00114353f, 0.00132105f, 0.00141205f, 0.00097077f, 0.00026421f, 0.00113901f, 0.00131767f, 0.00103704f, 0.00867996f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
69 {0.00318853f, 0.00138145f, 0.00104273f, 0.00105355f, 0.00094040f, 0.00100883f, 0.00124207f, 0.00142520f, 0.00059716f, 0.01778263f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
70 {0.00449576f, 0.00246811f, 0.00160275f, 0.00161966f, 0.00138494f, 0.00180553f, 0.00222063f, 0.00212853f, 0.00111754f, 0.01071834f, 0.03583921f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
71 {0.00331693f, 0.00595650f, 0.00257310f, 0.00252518f, 0.00046951f, 0.00312308f, 0.00428420f, 0.00259311f, 0.00121376f, 0.00157852f, 0.00259626f, 0.01612228f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
72 {0.00148878f, 0.00076734f, 0.00063401f, 0.00047808f, 0.00037421f, 0.00075546f, 0.00076105f, 0.00066504f, 0.00042237f, 0.00224097f, 0.00461939f, 0.00096120f, 0.00409522f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
73 {0.00165004f, 0.00090768f, 0.00084658f, 0.00069041f, 0.00052274f, 0.00059248f, 0.00078814f, 0.00115204f, 0.00072545f, 0.00279948f, 0.00533369f, 0.00087222f, 0.00116111f, 0.01661038f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
74 {0.00230618f, 0.00106268f, 0.00100282f, 0.00125381f, 0.00034766f, 0.00090111f, 0.00151550f, 0.00155601f, 0.00049078f, 0.00103767f, 0.00157310f, 0.00154836f, 0.00046718f, 0.00060701f, 0.01846071f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
75 {0.00631752f, 0.00224540f, 0.00301397f, 0.00285226f, 0.00094867f, 0.00191155f, 0.00293898f, 0.00381962f, 0.00116422f, 0.00173565f, 0.00250962f, 0.00312633f, 0.00087787f, 0.00119036f, 0.00180037f, 0.01346609f, 0.0f, 0.0f, 0.0f, 0.0f},
76 {0.00389995f, 0.00186053f, 0.00220144f, 0.00180488f, 0.00073798f, 0.00154526f, 0.00216760f, 0.00214841f, 0.00077747f, 0.00248968f, 0.00302273f, 0.00250862f, 0.00093371f, 0.00107595f, 0.00147982f, 0.00487295f, 0.01299436f, 0.0f, 0.0f, 0.0f},
77 {0.00039119f, 0.00029139f, 0.00021006f, 0.00016015f, 0.00010666f, 0.00020592f, 0.00023815f, 0.00038786f, 0.00019097f, 0.00039549f, 0.00076736f, 0.00028448f, 0.00016253f, 0.00085751f, 0.00015674f, 0.00026525f, 0.00024961f, 0.00563625f, 0.0f, 0.0f},
78 {0.00131840f, 0.00099430f, 0.00074960f, 0.00066005f, 0.00036626f, 0.00070192f, 0.00092548f, 0.00089301f, 0.00131038f, 0.00127857f, 0.00219713f, 0.00100817f, 0.00054105f, 0.00368739f, 0.00047608f, 0.00102648f, 0.00094759f, 0.00069226f, 0.00999315f, 0.0f},
79 {0.00533241f, 0.00169359f, 0.00136609f, 0.00127915f, 0.00119152f, 0.00132844f, 0.00178697f, 0.00194579f, 0.00071553f, 0.01117956f, 0.00914460f, 0.00210897f, 0.00197461f, 0.00256159f, 0.00135781f, 0.00241601f, 0.00343452f, 0.00038538f, 0.00148001f, 0.02075171f}
83 static int suboptimal_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int mode);
84 static int *** forward_so_dp ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
85 static int *** backward_so_dp ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
86 static int *** forward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
87 static int *** backward_so_dp_biphasic ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
88 static int *** forward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
89 static int *** backward_so_dp_glocal ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
97 float ** get_emitPairs (char *mat, char *alp, float **p, float *s);
98 int subop1_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
100 return suboptimal_pair_wise ( A, ns, ls, CL, 1);
103 int subop2_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
105 return suboptimal_pair_wise ( A, ns, ls, CL, 3);
110 int suboptimal_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int mode)
115 int gop, gep,gop2, gep2;
116 int i, I, j, J, n, s1, s2;
120 float opt, min, score, nscore, thres;
127 /*gop2=CL->gop*10*SCORE_K;*/
128 gop2=CL->gop*2*SCORE_K;
131 //Values Adapted from Probcons 1.1
138 ungap(A->seq_al[ls[0][0]]);
139 ungap(A->seq_al[ls[1][0]]);
141 seqI=A->seq_al[ls[0][0]];
142 seqJ=A->seq_al[ls[1][0]];
144 I=strlen (seqI); J=strlen (seqJ);
145 pos0=aln2pos_simple ( A,-1, ns, ls);
146 l1=strlen (A->seq_al[ls[0][0]]);
147 l2=strlen (A->seq_al[ls[1][0]]);
151 F=forward_so_dp (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
152 B=backward_so_dp (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
156 F=forward_so_dp_glocal (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
157 B=backward_so_dp_glocal (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
161 F=forward_so_dp_biphasic (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
162 B=backward_so_dp_biphasic (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
164 if ( MAX5(F[match][l1][l2], F[ins][l1][l2], F[del][l1][l2],F[ins2][l1][l2], F[del2][l1][l2] )!=MAX5( B[match][1][1], B[ins][1][1], B[del][1][1], B[ins2][1][1], B[del2][1][1]))
166 HERE ("ERROR in subop_pair");
167 fprintf ( stdout, "\nForward: %d", MAX3(F[match][l1][l2], F[ins][l1][l2], F[del][l1][l2]));
168 fprintf ( stdout, "\nBackWard: %d \n\n",MAX3( B[match][1][1], B[ins][1][1], B[del][1][1]));
172 for (opt=0,min=0, set=0, i=1; i<=I; i++)
175 if ( F[match][i][j]==UNDEFINED)continue;
176 F[match][i][j]+=B[match][i][j]-(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
178 {set=1; opt=F[match][i][j];min=F[match][i][j];}
179 opt=MAX(F[match][i][j],opt);
180 min=MIN(F[match][i][j],min);
184 s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
185 s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
187 id=idscore_pairseq(seqI,seqJ,-12, -1, CL->M, "idmat");
189 entry=vcalloc ( CL->entry_len, CL->el_size);
190 entry[SEQ1]=s1;entry[SEQ2]=s2;
193 for ( n=0,i=1; i<=I; i++)
198 nscore=((score-min))/(opt-min);
203 entry[R1]=i;entry[R2]=j;
206 add_entry2list (entry,A->CL);
218 /************************************************************************************************************************/
224 /************************************************************************************************************************/
225 int *** forward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop, int gep,int gop2, int gep2,Constraint_list *CL)
231 int match=0, del=1, ins=2;
233 M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
235 for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<5; c++)M[c][i][j]=-999999;
239 for (i=1; i<=I; i++){M[del] [i][0]=i*gep;M[umatch][i][0]=i*gep2+gop2;}
240 for (j=1; j<=J; j++){M[ins] [0][j]=j*gep;M[umatch][0][j]=j*gep2+gop2;}
245 for ( j=1; j<=J; j++)
247 sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
249 M[match][i][j] =MAX4 (M[match][i-1][j-1],M[del][i-1][j-1], M[ins][i-1][j-1],M[umatch][i-1][j-1])+sub;
250 M[del][i][j] =MAX2 ((M[match][i-1][j]+gop), M[del][i-1][j])+gep;
251 M[ins][i][j] =MAX2 ((M[match][i][j-1]+gop), M[ins][i][j-1])+gep;
252 M[umatch][i][j]=MAX6 (M[match][i-1][j-1]+gop2, M[match][i][j-1]+gop2, M[match][i-1][j]+gop2,M[umatch][i-1][j-1], M[umatch][i-1][j], M[umatch][i][j-1])+gep2;
257 int *** backward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
266 M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
267 for ( i=I+1; i>=0; i--)for (j=J+1; j>=0; j--)for (c=0; c<5; c++)M[c][i][j]=-999999;
268 M[match][I+1][J+1]=0;
270 for (i=I; i>0; i--){M[ins] [i][J+1]=i*gep;M[umatch] [i][J+1]=i*gep2+gop2;}
271 for (j=J; j>0; j--){M[del] [I+1][j]=j*gep;M[umatch] [I+1][j]=j*gep2+gop2;}
277 sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
279 M[match ][i][j] =MAX4 ((M[del][i+1][j+1]+gop), (M[ins][i+1][j+1]+gop), M[match][i+1][j+1], M[umatch][i+1][j+1]+gop2)+sub;
280 M[del ][i][j] =MAX2 (M[match][i+1][j], M[del][i+1][j])+gep;
281 M[ins ][i][j] =MAX2 (M[match][i][j+1], M[ins][i][j+1])+gep;
282 M[umatch][i][j] =MAX6 (M[match][i+1][j+1], M[match][i+1][j],M[match][i][j+1], M[umatch][i+1][j+1], M[umatch][i+1][j], M[umatch][i][j+1])+gep2;
292 /************************************************************************************************************************/
298 /************************************************************************************************************************/
300 int *** forward_so_dp ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop, int gep,int gop2, int gep2,Constraint_list *CL)
310 M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
311 for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<3; c++)M[c][i][j]=-999999;
314 for (i=1; i<=I; i++){M[del] [i][0]=i*gep;}
315 for (j=1; j<=J; j++){M[ins] [0][j]=j*gep;}
321 for ( j=1; j<=J; j++)
323 lgop=(i==I || j==J)?0:gop;
324 sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
326 M[match][i][j]=MAX3 (M[del][i-1][j-1], M[ins][i-1][j-1], M[match][i-1][j-1])+sub;
327 M[del][i][j] =MAX ((M[match][i-1][j]+lgop),M[del][i-1][j])+gep;
328 M[ins][i][j] =MAX ((M[match][i][j-1]+lgop), M[ins][i][j-1])+gep;
335 int *** backward_so_dp ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
344 for (b=0; b<ns[a]; b++)
346 invert_string2(A->seq_al[ls[a][b]]);
347 invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
349 T=forward_so_dp(A,ns,ls,pos0, I, J, gop, gep, gop2, gep2, CL);
351 for (b=0; b<ns[a]; b++)
353 invert_string2(A->seq_al[ls[a][b]]);
354 invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
357 M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
363 M[match][i+1][j+1]=T[match][I-i][J-j];
364 M[ins][i+1][j+1]=T[ins][I-i][J-j];
365 M[del][i+1][j+1]=T[del][I-i][J-j];
370 /************************************************************************************************************************/
376 /************************************************************************************************************************/
377 int biphasic_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
383 int M1, D1, D2, I1, I2, LEN;
386 char **al, **aln, *char_buf;
387 int gop1, gop2, gep1, gep2;
389 int score, trace, ntrace;
390 M1=n++; D1=n++; D2=n++; I1=n++, I2=n++;
392 I=strlen (A->seq_al[ls[0][0]]);
393 J=strlen (A->seq_al[ls[1][0]]);
394 m=declare_arrayN (3, sizeof (int),n, I+1, J+1);
395 t=declare_arrayN (3, sizeof (int),n, I+1, J+1);
396 pos0=aln2pos_simple ( A,-1, ns, ls);
397 al=declare_char (2, I+J+1);
398 for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<n; c++)m[c][i][j]=-999999;
400 gop1=CL->gop*SCORE_K*2;
401 gep1=CL->gep*SCORE_K/2;
403 gop2=CL->gop*SCORE_K/2;
404 gep2=CL->gep*SCORE_K*2;
407 for (i=1; i<=I; i++){m[I1][i][0]=gep1*i;}
408 for (j=1; j<=J; j++){m[D1][0][j]=gep1*j;}
410 for (i=1; i<=I; i++){m[I2] [i][0]=gep2*i;}
411 for (j=1; j<=J; j++){m[D2] [0][j]=gep2*j;}
415 for ( j=1; j<=J; j++)
417 sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
418 m[M1][i][j]=max_int (&t[M1][i][j],D1,m[D1][i-1][j-1],I1,m[I1][i-1][j-1], M1, m[M1][i-1][j-1],D2,m[D2][i-1][j-1],I2,m[I2][i-1][j-1], -1)+sub;
420 m[D1][i][j]=max_int (&t[D1][i][j],M1,(m[M1][i][j-1]+gop1),D1,m[D1][i][j-1], -1)+gep1;
421 m[I1][i][j]=max_int (&t[I1][i][j],M1,(m[M1][i-1][j]+gop1),I1,m[I1][i-1][j], -1)+gep1;
423 m[D2][i][j]=max_int (&t[D2][i][j],M1,(m[M1][i][j-1]+gop2),D2,m[D2][i][j-1], -1)+gep2;
424 m[I2][i][j]=max_int (&t[I2][i][j],M1,(m[M1][i-1][j]+gop2),I2,m[I2][i-1][j], -1)+gep2;
428 score=max_int (&trace,M1,m[M1][I][J],D1,m[D1][I][J],I1, m[I1][I][J],D2,m[D2][I][J],I2,m[I2][I][J], -1);
432 trace=t[trace][i][j];
433 while (!(i==0 &&j==0))
436 ntrace=t[trace][i][j];
459 else if ( trace==D1 || trace==D2)
466 else if ( trace == I1 || trace==I2)
477 invert_list_char ( al[0], LEN);
478 invert_list_char ( al[1], LEN);
479 if ( A->declared_len<=LEN)A=realloc_aln2 ( A,A->max_n_seq, 2*LEN);
482 char_buf= vcalloc (LEN+1, sizeof (char));
483 for ( c=0; c< 2; c++)
485 for ( a=0; a< ns[c]; a++)
488 for ( b=0; b< LEN; b++)
491 char_buf[b]=aln[ls[c][a]][ch++];
496 sprintf (aln[ls[c][a]],"%s", char_buf);
503 free_arrayN((void *)m, 3);
504 free_arrayN((void *)t, 3);
509 int *** forward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop1, int gep1,int gop2, int gep2,Constraint_list *CL)
515 int match=0, del=1, ins=2;
516 int lgop1, lgop2, lgep1, lgep2;
518 M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
520 for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<5; c++)M[c][i][j]=-999999;
524 for (i=1; i<=I; i++){M[del] [i][0]=gep1*i+gop1;}
525 for (j=1; j<=J; j++){M[ins] [0][j]=gep1*j+gop1;}
527 for (i=1; i<=I; i++){M[del2] [i][0]=gep2*i+gop2;}
528 for (j=1; j<=J; j++){M[ins2] [0][j]=gep2*j+gop2;}
532 for ( j=1; j<=J; j++)
534 lgop1=(i==I || j==J)?gop1:gop1;
535 lgop2=(i==I || j==J)?gop2:gop2;
539 sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
540 M[match][i][j]=MAX5 (M[del][i-1][j-1], M[ins][i-1][j-1], M[match][i-1][j-1], M[ins2][i-1][j-1], M[del2][i-1][j-1])+sub;
542 M[del ][i][j] =MAX2 ((M[match][i-1][j]+lgop1), M[del ][i-1][j])+lgep1;
543 M[del2][i][j] =MAX2 ((M[match][i-1][j]+lgop2), M[del2][i-1][j])+lgep2;
545 M[ins ][i][j] =MAX2 ((M[match][i][j-1]+lgop1), M[ins ][i][j-1] )+lgep1;
546 M[ins2][i][j] =MAX2 ((M[match][i][j-1]+lgop2), M[ins2][i][j-1] )+lgep2;
551 int *** backward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
560 for (b=0; b<ns[a]; b++)
562 invert_string2(A->seq_al[ls[a][b]]);
563 invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
565 T=forward_so_dp_biphasic(A,ns,ls,pos0, I, J, gop, gep, gop2, gep2, CL);
567 for (b=0; b<ns[a]; b++)
569 invert_string2(A->seq_al[ls[a][b]]);
570 invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
573 M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
579 M[match][i+1][j+1]=T[match][I-i][J-j];
580 M[ins][i+1][j+1]=T[ins][I-i][J-j];
581 M[del][i+1][j+1]=T[del][I-i][J-j];
582 M[ins2][i+1][j+1]=T[ins2][I-i][J-j];
583 M[del2][i+1][j+1]=T[del2][I-i][J-j];
590 int get_tot_prob (Alignment *A1,Alignment *A2, int *ns, int **ls, int nstates, float **matchProb, float **insProb, float *TmatchProb, float ***TinsProb, Constraint_list *CL);
592 float * forward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *TmatchProb, float ***TinsProb, float **transProb);
593 float * backward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *TmatchProb, float ***TinsProb,float **transProb);
594 float ComputeTotalProbability (int seq1Length, int seq2Length,int NumMatrixTypes, int NumInsertStates,float *forward, float *backward) ;
595 int ProbabilisticModel (int NumMatrixTypes, int NumInsertStates,float *initDistribMat,float *emitSingle, float** emitPairs, float *gapOpen, float *gapExtend, float **transMat, float *initialDistribution, float **matchProb, float **insProb, float **transProb);
597 Constraint_list *ProbaMatrix2CL (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL);
600 void free_proba_pair_wise ()
602 proba_pair_wise (NULL, NULL, NULL, NULL);
604 int proba_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
606 int NumMatrixTypes=5;
607 int NumInsertStates=2;
608 static float **transMat, **insProb, **matchProb, *initialDistribution, **transProb, **emitPairs, *emitSingle, ***TinsProb, *TmatchProb;
609 static int TinsProb_ml, TmatchProb_ml;
614 float thr=0.01;//ProbCons Default
619 //Free all the memory
622 free_float (transMat, -1);transMat=NULL;
623 free_float (insProb, -1);insProb=NULL;
624 free_float (matchProb, -1);matchProb=NULL;
625 vfree (initialDistribution); initialDistribution=NULL;
626 free_float (transProb, -1);transProb=NULL;
627 free_float (emitPairs, -1);emitPairs=NULL;
628 vfree (emitSingle);emitSingle=NULL;
631 free_arrayN((void***)TinsProb, 3);TinsProb=NULL;
632 vfree (TmatchProb);TmatchProb=NULL;
633 TinsProb_ml=0; TmatchProb_ml=0;
635 forward_proba_pair_wise (NULL, NULL, 0,0,NULL,NULL,NULL,NULL,NULL);
636 backward_proba_pair_wise (NULL, NULL, 0,0,NULL,NULL,NULL,NULL,NULL);
637 ProbaMatrix2CL(NULL, NULL, NULL, 0, 0, NULL, NULL, 0, NULL);
641 if (!transMat && (strm (retrieve_seq_type(), "DNA") ||strm (retrieve_seq_type(), "RNA")) )
648 l=strlen (DNAalphabetDefault);
649 p=declare_float (l,l);
650 s=vcalloc (l, sizeof (float));
653 s[a]=DNAemitSingleDefault[a];
655 p[a][b]=DNAemitPairsDefault[a][b];
658 p=get_emitPairs (CL->method_matrix, DNAalphabetDefault,p,s);
659 alphabet=DNAalphabetDefault;
660 emitPairs=declare_float (256, 256);
661 emitSingle=vcalloc (256, sizeof (float));
662 for (i=0; i<256; i++)
665 for (j=0; j<256; j++)
666 emitPairs[i][j]=1e-10;
673 c1=tolower(alphabet[i]);
674 C1=toupper(alphabet[i]);
679 c2=tolower(alphabet[j]);
680 C2=toupper(alphabet[j]);
682 emitPairs[c1][c2]=p[i][j];
683 emitPairs[C1][c2]=p[i][j];
684 emitPairs[C1][C2]=p[i][j];
685 emitPairs[c1][C2]=p[i][j];
686 emitPairs[c2][c1]=p[i][j];
687 emitPairs[C2][c1]=p[i][j];
688 emitPairs[C2][C1]=p[i][j];
689 emitPairs[c2][C1]=p[i][j];
694 transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
695 transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
696 insProb=declare_float (256,NumMatrixTypes);
697 matchProb=declare_float (256, 256);
698 initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
700 ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,DNAgapOpen2Default,DNAgapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
702 else if ( !transMat && strm (retrieve_seq_type(), "PROTEIN"))
709 l=strlen (alphabetDefault);
710 p=declare_float (l,l);
711 s=vcalloc (l, sizeof (float));
714 s[a]=emitSingleDefault[a];
716 p[a][b]=emitPairsDefault[a][b];
719 p=get_emitPairs (CL->method_matrix, alphabetDefault,p,s);
720 alphabet=alphabetDefault;
721 emitPairs=declare_float (256, 256);
722 emitSingle=vcalloc (256, sizeof (float));
723 for (i=0; i<256; i++)
725 //emitSingle[i]=1e-5;
727 for (j=0; j<256; j++)
728 //emitPairs[i][j]=1e-10;
737 c1=tolower(alphabet[i]);
738 C1=toupper(alphabet[i]);
743 c2=tolower(alphabet[j]);
744 C2=toupper(alphabet[j]);
746 emitPairs[c1][c2]=p[i][j];
747 emitPairs[C1][c2]=p[i][j];
748 emitPairs[C1][C2]=p[i][j];
749 emitPairs[c1][C2]=p[i][j];
750 emitPairs[c2][c1]=p[i][j];
751 emitPairs[C2][c1]=p[i][j];
752 emitPairs[C2][C1]=p[i][j];
753 emitPairs[c2][C1]=p[i][j];
759 transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
760 transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
761 insProb=declare_float (256,NumMatrixTypes);
762 matchProb=declare_float (256, 256);
763 initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
765 ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,gapOpen2Default,gapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
768 I=strlen (A->seq_al[ls[0][0]]);
769 J=strlen (A->seq_al[ls[1][0]]);
770 //TmatchProb=vcalloc ((I+1)*(J+1), sizeof (float));
771 //TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,MAX(I,J)+1);
777 if (TmatchProb)TmatchProb=vrealloc(TmatchProb,TmatchProb_ml*sizeof (float));
778 else TmatchProb=vcalloc ( l, sizeof (float));
784 if (TinsProb)free_arrayN (TinsProb, 3);
785 TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,TinsProb_ml);
788 get_tot_prob (A,A, ns,ls,NumMatrixTypes, matchProb, insProb,TmatchProb,TinsProb, CL);
790 F=forward_proba_pair_wise (A->seq_al[ls[0][0]], A->seq_al[ls[1][0]], NumMatrixTypes,NumInsertStates,transMat, initialDistribution,TmatchProb,TinsProb, transProb);
791 B=backward_proba_pair_wise (A->seq_al[ls[0][0]], A->seq_al[ls[1][0]], NumMatrixTypes,NumInsertStates,transMat, initialDistribution,TmatchProb,TinsProb, transProb);
792 A->CL=ProbaMatrix2CL(A,ns, ls,NumMatrixTypes,NumInsertStates, F, B, thr,CL);
794 //free_proba_pair_wise();
798 int get_tot_prob (Alignment *A1,Alignment *A2, int *ns, int **ls, int nstates, float **matchProb, float **insProb, float *TmatchProb, float ***TinsProb, Constraint_list *CL)
800 int i, j, a, b, c,d, k, n,n1,n2, ij;
803 int ***VA1,***VA2, *observed, index;
804 //Pre-computation of the pairwise scores in order to use potential profiles
805 //The profiles are vectorized AND Compressed so that the actual alphabet size (proteins/DNA) does not need to be considered
808 if (ns[0]==1 && ns[1]==1 )
812 Alignment *NA1, *NA2;
814 nns=vcalloc ( 2, sizeof (int));
815 nls=vcalloc (2, sizeof (int*));
817 s1=A1->order[ls[0][0]][0];
818 s2=A2->order[ls[1][0]][0];
819 NA1=seq2R_template_profile (CL->S,s1);
820 NA2=seq2R_template_profile (CL->S,s2);
827 nls[0]=vcalloc (NA1->nseq, sizeof (int));
828 for (a=0; a<NA1->nseq; a++)
835 nls[0]=vcalloc (ns[0], sizeof (int));
836 for (a=0; a<ns[0]; a++)
843 nls[1]=vcalloc (NA2->nseq, sizeof (int));
844 for (a=0; a<NA2->nseq; a++)
851 nls[1]=vcalloc (ns[1], sizeof (int));
852 for (a=0; a<ns[1]; a++)
856 get_tot_prob (NA1, NA2, nns, nls, nstates, matchProb, insProb, TmatchProb, TinsProb, CL);
857 vfree (nns); free_int (nls,-1);
862 I=strlen (A1->seq_al[ls[0][0]]);
863 J=strlen (A2->seq_al[ls[1][0]]);
870 for (k=0; k<nstates; k++)
873 for (n=0,b=0; b<ns[0]; b++)
875 c1=A1->seq_al[ls[0][b]][i-1];
878 TinsProb[0][k][i]+=insProb[c1][k];
883 if (n)TinsProb[0][k][i]/=n;
889 for (k=0; k<nstates; k++)
892 for (n=0,b=0; b<ns[1]; b++)
894 c2=A2->seq_al[ls[1][b]][j-1];
897 TinsProb[1][k][j]+=insProb[c2][k];
902 if (n)TinsProb[1][k][j]/=n;
906 observed=vcalloc ( 26, sizeof (int));
907 VA1=declare_arrayN (3, sizeof (int),2,26,I);
910 for (index=0, b=0; b<ns[0]; b++)
913 c1=tolower(A1->seq_al[ls[0][b]][i]);
914 if ( c1=='-')continue;
917 if (!(in=observed[c1])){in=observed[c1]=++index;}
924 for (b=0; b<26; b++)observed[b]=0;
927 VA2=declare_arrayN (3, sizeof (int),2,26,J);
930 for (index=0, b=0; b<ns[1]; b++)
934 c1=tolower(A2->seq_al[ls[1][b]][i]);
935 if ( c1=='-')continue;
938 if (!(in=observed[c1])){in=observed[c1]=++index;}
944 for (b=0; b<26; b++)observed[b]=0;
948 for ( ij=0,i=0; i<=I; i++)
950 for ( j=0; j<=J ; j++, ij++)
958 while (VA1[0][c][i-1]!=-1)
960 c1=VA1[0][c][i-1]+'a';
963 while (VA2[0][d][j-1]!=-1)
965 c2=VA2[0][d][j-1]+'a';
967 TmatchProb[ij]+=matchProb[c1][c2]*(double)n1*(double)n2;
974 if (n)TmatchProb[ij]/=n;
978 free_arrayN ((void **)VA1, 3);
979 free_arrayN ((void **)VA2, 3);
986 Constraint_list *ProbaMatrix2CL (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL)
989 int ij, i, j,k, I, J, s1, s2;
998 static float F=4; //potential number of full suboptimal alignmnents incorporated in the library
999 static int tot_old, tot_new;
1003 free_int (list, -1);list=NULL;
1006 vfree(entry); entry=NULL;
1010 I=strlen (A->seq_al[ls[0][0]]);
1011 J=strlen (A->seq_al[ls[1][0]]);
1012 s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
1013 s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
1017 if ( list_max<list_size)
1019 free_int (list, -1);
1021 list=declare_int (list_max, 3);
1025 totalProb = ComputeTotalProbability (I,J,NumMatrixTypes, NumInsertStates,forward, backward);
1028 for (list_n=0,ij=0,i =0; i <= I; i++)
1030 for (j =0; j <= J; j++, ij+=NumMatrixTypes)
1032 v= EXP (MIN(LOG_ONE,(forward[ij] + backward[ij] - totalProb)));
1033 if (v>thr)//Conservative reduction of the list size to speed up the sorting
1037 list[list_n][2]=(int)((float)v*(float)NORM_F);
1044 sort_int_inv (list, 3, 2, 0, list_n-1);
1045 if (!entry)entry=vcalloc ( CL->entry_len, CL->el_size);
1046 list_n=MIN(list_n,(F*MIN(I,J)));
1047 for (i=0; i<list_n; i++)
1051 entry[R1] =list[i][0];
1052 entry[R2] =list[i][1];
1053 entry[WE] =list[i][2];
1055 add_entry2list (entry,A->CL);
1059 // HERE ("LIB_SIZE NEW: %d (new) %d (old) [%.2f]", list_n, old_n, (float)tot_new/(float)tot_old);
1063 Constraint_list *ProbaMatrix2CL_old (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL)
1066 int ij, i, j,k, I, J, s1, s2;
1071 I=strlen (A->seq_al[ls[0][0]]);
1072 J=strlen (A->seq_al[ls[1][0]]);
1073 s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
1074 s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
1076 totalProb = ComputeTotalProbability (I,J,NumMatrixTypes, NumInsertStates,forward, backward);
1077 if (!entry)entry=vcalloc ( CL->entry_len, CL->el_size);
1082 for (ij=0,i =0; i <= I; i++)
1084 for (j =0; j <= J; j++)
1086 v= EXP (MIN(LOG_ONE,(forward[ij] + backward[ij] - totalProb)));
1087 if (i && j && v>=thr)
1090 entry[SEQ1]=s1;entry[SEQ2]=s2;
1091 entry[R1]=i;entry[R2]=j;
1092 entry[WE]=(int)((float)v*(float)NORM_F);
1094 add_entry2list (entry,A->CL);
1097 ij += NumMatrixTypes;
1100 HERE ("LIB_SIZE_OLD: %d", lib_size);
1104 float ComputeTotalProbability (int seq1Length, int seq2Length,int NumMatrixTypes, int NumInsertStates,float *forward, float *backward)
1107 float totalForwardProb = LOG_ZERO;
1108 float totalBackwardProb = LOG_ZERO;
1111 for (k = 0; k < NumMatrixTypes; k++)
1113 LOG_PLUS_EQUALS (&totalForwardProb,forward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] + backward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
1116 totalBackwardProb =forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] +backward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)];
1118 for (k = 0; k < NumInsertStates; k++)
1120 LOG_PLUS_EQUALS (&totalBackwardProb,forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] +backward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)]);
1121 LOG_PLUS_EQUALS (&totalBackwardProb,forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] +backward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)]);
1123 return (totalForwardProb + totalBackwardProb) / 2;
1127 float * backward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *matchProb, float ***insProb, float **transProb)
1129 static float *backward;
1133 int k, i, j,ij, i1j1, i1j, ij1,a, l, seq1Length, seq2Length, m;
1135 char *iter1, *iter2;
1140 backward=NULL; max_l=0;
1146 seq1Length=strlen (seq1);
1147 seq2Length=strlen (seq2);
1148 l=(seq1Length+1)*(seq2Length+1)*NumMatrixTypes;
1152 backward=vcalloc (l, sizeof (float));
1157 backward=vrealloc (backward, l*sizeof(float));
1161 for (a=0; a<l; a++)backward[a]=LOG_ZERO;
1163 for (k = 0; k < NumMatrixTypes; k++)
1164 backward[NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1) + k] = initialDistribution[k];
1166 // remember offset for each index combination
1167 ij = (seq1Length+1) * (seq2Length+1) - 1;
1169 i1j = ij + seq2Length + 1;
1171 i1j1 = ij + seq2Length + 2;
1172 ij *= NumMatrixTypes;
1173 i1j *= NumMatrixTypes;
1174 ij1 *= NumMatrixTypes;
1175 i1j1 *= NumMatrixTypes;
1177 // compute backward scores
1178 for (i = seq1Length; i >= 0; i--)
1180 c1 = (i == seq1Length) ? '~' : (unsigned char) iter1[i+1];
1181 for (j = seq2Length; j >= 0; j--)
1183 c2 = (j == seq2Length) ? '~' : (unsigned char) iter2[j+1];
1185 if (i < seq1Length && j < seq2Length)
1187 m=((i+1)*(seq2Length+1))+j+1;//The backward and the forward are offset by 1
1188 float ProbXY = backward[0 + i1j1] + matchProb[m];
1191 for (k = 0; k < NumMatrixTypes; k++)
1193 LOG_PLUS_EQUALS (&backward[k + ij], ProbXY + transProb[k][0]);
1198 for (k = 0; k < NumInsertStates; k++)
1200 LOG_PLUS_EQUALS (&backward[0 + ij], backward[2*k+1 + i1j] + insProb[0][k][i+1] + transProb[0][2*k+1]);
1201 LOG_PLUS_EQUALS (&backward[2*k+1 + ij], backward[2*k+1 + i1j] + insProb[0][k][i+1] + transProb[2*k+1][2*k+1]);
1206 for (k = 0; k < NumInsertStates; k++)
1208 //+1 because the backward and the forward are offset by 1
1209 LOG_PLUS_EQUALS (&backward[0 + ij], backward[2*k+2 + ij1] + insProb[1][k][j+1] + transProb[0][2*k+2]);
1210 LOG_PLUS_EQUALS (&backward[2*k+2 + ij], backward[2*k+2 + ij1] + insProb[1][k][j+1] + transProb[2*k+2][2*k+2]);
1214 ij -= NumMatrixTypes;
1215 i1j -= NumMatrixTypes;
1216 ij1 -= NumMatrixTypes;
1217 i1j1 -= NumMatrixTypes;
1223 float * forward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *matchProb, float ***insProb, float **transProb)
1225 static float *forward;
1227 int k, i, j,ij, i1j1, i1j, ij1, seq1Length, seq2Length, m;
1228 char *iter1, *iter2;
1234 forward=NULL; max_l=0;
1239 seq1Length=strlen (seq1);
1240 seq2Length=strlen (seq2);
1241 l=(seq1Length+1)*(seq2Length+1)*NumMatrixTypes;
1245 forward=vcalloc (l, sizeof (float));
1250 forward=vrealloc (forward, l*sizeof(float));
1253 for (a=0; a<l; a++)forward[a]=LOG_ZERO;
1256 forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] = initialDistribution[0] + matchProb[seq2Length+2];
1258 for (k = 0; k < NumInsertStates; k++)
1260 forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] = initialDistribution[2*k+1] + insProb[0][k][1];
1261 forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] = initialDistribution[2*k+2] + insProb[1][k][1];
1264 // remember offset for each index combination
1266 i1j = -seq2Length - 1;
1268 i1j1 = -seq2Length - 2;
1270 ij *= NumMatrixTypes;
1271 i1j *= NumMatrixTypes;
1272 ij1 *= NumMatrixTypes;
1273 i1j1 *= NumMatrixTypes;
1276 // compute forward scores
1277 for (m=0,i = 0; i <= seq1Length; i++)
1279 for (j = 0; j <= seq2Length; j++, m++)
1285 //Sum over all possible alignments
1286 forward[0 + ij] = forward[0 + i1j1] + transProb[0][0];
1287 for (k = 1; k < NumMatrixTypes; k++)
1289 LOG_PLUS_EQUALS (&forward[0 + ij], forward[k + i1j1] + transProb[k][0]);
1291 forward[0 + ij] += matchProb[m];
1295 for (k = 0; k < NumInsertStates; k++)
1297 forward[2*k+1 + ij] = insProb[0][k][i] + LOG_ADD (forward[0 + i1j] + transProb[0][2*k+1],forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1]);
1302 for (k = 0; k < NumInsertStates; k++)
1304 forward[2*k+2 + ij] = insProb[1][k][j] +LOG_ADD (forward[0 + ij1] + transProb[0][2*k+2],forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2]);
1309 ij += NumMatrixTypes;
1310 i1j += NumMatrixTypes;
1311 ij1 += NumMatrixTypes;
1312 i1j1 += NumMatrixTypes;
1318 int ProbabilisticModel (int NumMatrixTypes, int NumInsertStates,float *initDistribMat,float *emitSingle, float **emitPairs, float *gapOpen, float *gapExtend, float **transMat, float *initialDistribution, float **matchProb, float **insProb, float **transProb)
1322 // build transition matrix
1327 for (i = 0; i < NumInsertStates; i++)
1329 transMat[0][2*i+1] = gapOpen[2*i];
1330 transMat[0][2*i+2] = gapOpen[2*i+1];
1331 transMat[0][0] -= (gapOpen[2*i] + gapOpen[2*i+1]);
1333 transMat[2*i+1][2*i+1] = gapExtend[2*i];
1334 transMat[2*i+2][2*i+2] = gapExtend[2*i+1];
1335 transMat[2*i+1][2*i+2] = 0;
1336 transMat[2*i+2][2*i+1] = 0;
1337 transMat[2*i+1][0] = 1 - gapExtend[2*i];
1338 transMat[2*i+2][0] = 1 - gapExtend[2*i+1];
1343 // create initial and transition probability matrices
1344 for (i = 0; i < NumMatrixTypes; i++){
1345 initialDistribution[i] = (float)log ((float)initDistribMat[i]);
1346 for (j = 0; j < NumMatrixTypes; j++)
1347 transProb[i][j] = (float)log ((float)transMat[i][j]);
1350 // create insertion and match probability matrices
1351 for (i = 0; i < 256; i++)
1353 for (j = 0; j < NumMatrixTypes; j++)
1355 insProb[i][j] = (float)log((float)emitSingle[i]);
1357 for (j = 0; j < 256; j++)
1359 matchProb[i][j] = (float)log((float)emitPairs[i][j]);
1366 int viterbi_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1369 char *alphabet, *char_buf;
1371 int seq1Length, seq2Length, I, J;
1372 int i, j,ij, i1j1, i1j, ij1, k, a, b,l, LEN, r, c, m, state;
1373 int NumMatrixTypes=5;
1374 int NumInsertStates=2;
1377 static float **transMat, **insProb, **matchProb, *initialDistribution, **transProb, **emitPairs, *emitSingle, *TmatchProb, ***TinsProb;
1380 ungap_sub_aln (A, ns[0],ls[0]);
1381 ungap_sub_aln (A, ns[1],ls[1]);
1383 seq1Length=I=strlen (A->seq_al[ls[0][0]]);
1384 seq2Length=J=strlen (A->seq_al[ls[1][0]]);
1388 alphabet=alphabetDefault;
1389 emitPairs=declare_float (256, 256);
1390 emitSingle=vcalloc (256, sizeof (float));
1391 for (i=0; i<256; i++)
1394 for (j=0; j<256; j++)
1395 emitPairs[i][j]=1e-10;
1397 l=strlen (alphabet);
1402 c1=tolower(alphabet[i]);
1403 C1=toupper(alphabet[i]);
1404 emitSingle[c1]=emitSingleDefault[i];
1405 emitSingle[C1]=emitSingleDefault[i];
1406 for (j=0; j<=i; j++)
1408 c2=tolower(alphabet[j]);
1409 C2=toupper(alphabet[j]);
1411 emitPairs[c1][c2]=emitPairsDefault[i][j];
1412 emitPairs[C1][c2]=emitPairsDefault[i][j];
1413 emitPairs[C1][C2]=emitPairsDefault[i][j];
1414 emitPairs[c1][C2]=emitPairsDefault[i][j];
1415 emitPairs[c2][c1]=emitPairsDefault[i][j];
1416 emitPairs[C2][c1]=emitPairsDefault[i][j];
1417 emitPairs[C2][C1]=emitPairsDefault[i][j];
1418 emitPairs[c2][C1]=emitPairsDefault[i][j];
1423 transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
1424 transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
1425 insProb=declare_float (256,NumMatrixTypes);
1426 matchProb=declare_float (256, 256);
1427 initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
1429 ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,gapOpen2Default,gapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
1433 TmatchProb=vcalloc ((I+1)*(J+1), sizeof (float));
1434 TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,MAX(I,J)+1);
1435 get_tot_prob (A,A, ns,ls,NumMatrixTypes, matchProb, insProb,TmatchProb,TinsProb, CL);
1437 // create viterbi matrix
1438 l=NumMatrixTypes * (seq1Length+1) * (seq2Length+1);
1439 viterbi =vcalloc (l, sizeof (float));
1440 for (a=0; a<l; a++)viterbi[a]=LOG_ZERO;
1441 traceback=vcalloc (l, sizeof (int));
1442 for (a=0; a<l; a++)traceback[a]=-1;
1444 // initialization condition
1445 for (k = 0; k < NumMatrixTypes; k++)
1446 viterbi[k] = initialDistribution[k];
1448 // remember offset for each index combination
1450 i1j = -seq2Length - 1;
1452 i1j1 = -seq2Length - 2;
1454 ij *= NumMatrixTypes;
1455 i1j *= NumMatrixTypes;
1456 ij1 *= NumMatrixTypes;
1457 i1j1 *= NumMatrixTypes;
1459 // compute viterbi scores
1460 for (m=0,i = 0; i <= seq1Length; i++)
1462 for ( j = 0; j <= seq2Length; j++, m++)
1466 for (k = 0; k < NumMatrixTypes; k++)
1468 float newVal = viterbi[k + i1j1] + transProb[k][0] + TmatchProb[m];
1469 if (viterbi[0 + ij] < newVal)
1471 viterbi[0 + ij] = newVal;
1472 traceback[0 + ij] = k;
1478 for (k = 0; k < NumInsertStates; k++)
1480 float valFromMatch = TinsProb[0][k][i] + viterbi[0 + i1j] + transProb[0][2*k+1];
1481 float valFromIns = TinsProb[0][k][i] + viterbi[2*k+1 + i1j] + transProb[2*k+1][2*k+1];
1482 if (valFromMatch >= valFromIns){
1483 viterbi[2*k+1 + ij] = valFromMatch;
1484 traceback[2*k+1 + ij] = 0;
1487 viterbi[2*k+1 + ij] = valFromIns;
1488 traceback[2*k+1 + ij] = 2*k+1;
1494 for (k = 0; k < NumInsertStates; k++){
1495 float valFromMatch = TinsProb[1][k][j] + viterbi[0 + ij1] + transProb[0][2*k+2];
1496 float valFromIns = TinsProb[1][k][j] + viterbi[2*k+2 + ij1] + transProb[2*k+2][2*k+2];
1497 if (valFromMatch >= valFromIns){
1498 viterbi[2*k+2 + ij] = valFromMatch;
1499 traceback[2*k+2 + ij] = 0;
1503 viterbi[2*k+2 + ij] = valFromIns;
1504 traceback[2*k+2 + ij] = 2*k+2;
1509 ij += NumMatrixTypes;
1510 i1j += NumMatrixTypes;
1511 ij1 += NumMatrixTypes;
1512 i1j1 += NumMatrixTypes;
1516 // figure out best terminating cell
1517 bestProb = LOG_ZERO;
1519 for (k = 0; k < NumMatrixTypes; k++)
1521 float thisProb = viterbi[k + NumMatrixTypes * ((seq1Length+1)*(seq2Length+1) - 1)] + initialDistribution[k];
1522 if (bestProb < thisProb)
1524 bestProb = thisProb;
1531 // compute traceback
1532 al=declare_char(2,seq1Length+seq2Length);
1534 r = seq1Length, c = seq2Length;
1535 while (r != 0 || c != 0)
1537 int newState = traceback[state + NumMatrixTypes * (r * (seq2Length+1) + c)];
1539 if (state == 0){ c--; r--; al[0][LEN]=1;al[1][LEN]=1;}
1540 else if (state % 2 == 1) {r--; al[0][LEN]=1;al[1][LEN]=0;}
1541 else { c--; al[0][LEN]=0;al[1][LEN]=1;}
1547 invert_list_char ( al[0], LEN);
1548 invert_list_char ( al[1], LEN);
1549 if ( A->declared_len<=LEN)A=realloc_aln2 ( A,A->max_n_seq, 2*LEN);
1551 char_buf= vcalloc (LEN+1, sizeof (char));
1552 for ( c=0; c< 2; c++)
1554 for ( a=0; a< ns[c]; a++)
1557 for ( b=0; b< LEN; b++)
1560 char_buf[b]=aln[ls[c][a]][ch++];
1565 sprintf (aln[ls[c][a]],"%s", char_buf);
1571 A->nseq=ns[0]+ns[1];
1579 return (int)(bestProb*(float)1000);
1582 float ** get_emitPairs (char *mat, char *alp, float **p, float *s)
1589 if (!rmat)rmat=vcalloc (100, sizeof (char));
1591 if (!mat || !mat[0] || strm (mat, "default"))return p;
1592 else if (strm (rmat, mat))return p;
1594 sprintf (rmat,"%s", mat);
1596 M=read_matrice (mat);
1606 sc=M[alp[a]-'A'][alp[b]-'A'];
1607 p[a][b]=e*exp ((double)sc*k);
1627 /*********************************COPYRIGHT NOTICE**********************************/
1628 /*© Centro de Regulacio Genomica */
1630 /*Cedric Notredame */
1631 /*Tue Oct 27 10:12:26 WEST 2009. */
1632 /*All rights reserved.*/
1633 /*This file is part of T-COFFEE.*/
1635 /* T-COFFEE is free software; you can redistribute it and/or modify*/
1636 /* it under the terms of the GNU General Public License as published by*/
1637 /* the Free Software Foundation; either version 2 of the License, or*/
1638 /* (at your option) any later version.*/
1640 /* T-COFFEE is distributed in the hope that it will be useful,*/
1641 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1642 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
1643 /* GNU General Public License for more details.*/
1645 /* You should have received a copy of the GNU General Public License*/
1646 /* along with Foobar; if not, write to the Free Software*/
1647 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
1648 /*............................................... |*/
1649 /* If you need some more information*/
1650 /* cedric.notredame@europe.com*/
1651 /*............................................... |*/
1655 /*********************************COPYRIGHT NOTICE**********************************/