1 /***********************************************
2 * # Copyright 2009-2010. Liu Yongchao
3 * # Contact: Liu Yongchao, School of Computer Engineering,
4 * # Nanyang Technological University.
5 * # Emails: liuy0039@ntu.edu.sg; nkcslyc@hotmail.com
7 * # GPL version 3.0 applies.
9 * ************************************************/
10 #include "SafeVector.h"
19 #include "MultiSequence.h"
20 #include "ScoreType.h"
22 #define TRACE 0 // 0: NOTRACE 1: TRACE
24 #define endgaps 1 // 1: engap penaties enabled 0: disabled
25 #define PART_FULL_MEMORY 0 //0: LOW MEM OPTION
26 #define REVPART_FULL_MEMORY 0 //0: LOW MEM OPTION
30 #define OS_HUGE_VALL HUGE_VAL
32 #define OS_HUGE_VALL HUGE_VALL
41 char opt; //can be 'P' or 'M'
46 typedef struct sequence {
52 typedef struct alignment {
58 ////////////////////////////////////////////////////////
59 //externs related to scoring matrix and input arguments
60 ///////////////////////////////////////////////////////////
61 extern float g_gap_open1, g_gap_open2, g_gap_ext1, g_gap_ext2;
62 extern char aminos[26], matrixtype[20], bases[26];
64 extern double sub_matrix[26][26];
65 extern double normalized_matrix[26][26]; // add by YE Yongtao
66 extern int subst_index[26];
68 extern float TEMPERATURE;
69 extern int MATRIXTYPE;
73 extern argument_decl argument;
75 //////////////////////////////////////////////////////////////////////////////
76 //calculates reverse partition function values based on z matrices
77 //and also simulaneously calculates the propability of each basepair
78 //or aminoacid residue pair i,j
79 //////////////////////////////////////////////////////////////////////////////
81 VF *revers_partf(fasta sequences[2], const double termgapopen,
82 const double termgapextend, long double **Zfm, const double d,
84 // printf("revpart\n");
85 //rest of the declarations
87 long double **Zm = NULL;
88 long double **Ze = NULL;
89 long double **Zf = NULL;
94 double endgapopen, endgapextend;
97 //Init lengths of sequences
98 len0 = strlen(sequences[0].text);
99 len1 = strlen(sequences[1].text);
101 //Safe vector declared
102 VF *posteriorPtr = new VF((len0 + 1) * (len1 + 1));
103 VF & posterior = *posteriorPtr;
104 VF::iterator ptr = posterior.begin();
106 if (TRACE) //open the trace file
107 fo = fopen("revpartdump", "a");
110 endgapopen = termgapopen;
111 endgapextend = termgapextend;
113 //instantiate the z matrix
114 if (REVPART_FULL_MEMORY) {
116 Ze = new long double *[sequences[1].length + 1];
117 Zf = new long double *[sequences[1].length + 1];
118 Zm = new long double *[sequences[1].length + 1];
121 printf("\n\n %e %e\n", d, e);
123 //DYNAMICALLY GROW 2D Zm Zf Ze MARICES (long double)
124 for (i = 0; i <= sequences[1].length; i++) {
125 Ze[i] = new long double[sequences[0].length + 1];
126 Zf[i] = new long double[sequences[0].length + 1];
127 Zm[i] = new long double[sequences[0].length + 1];
130 Zm = new long double *[2];
131 Ze = new long double *[2];
132 Zf = new long double *[2];
133 for (i = 0; i <= 1; i++) {
134 Zm[i] = new long double[sequences[0].length + 1];
135 Ze[i] = new long double[sequences[0].length + 1];
136 Zf[i] = new long double[sequences[0].length + 1];
142 printf("in rev partf---");
146 if (REVPART_FULL_MEMORY) {
147 for (i = 0; i <= len1; i++)
148 for (j = 0; j <= len0; j++) {
155 for (j = 0; j <= len0; j++) {
165 //fill the probability matrix with 0s
166 for (i = 0; i <= len1; i++)
167 for (j = 0; j <= len0; j++)
168 ptr[j * (len1 + 1) + i] = 0;
172 Ze[len1][len0] = Zf[len1][len0] = 0;
173 Zf[len1 - 1][len0] = Zm[len1][len0] * d;
174 Ze[len1][len0 - 1] = Zm[len1][len0] * d;
177 if (REVPART_FULL_MEMORY) {
178 for (i = len1 - 2; i >= 0; i--) {
179 Zf[i][len0] = Zf[i + 1][len0] * e;
184 if (REVPART_FULL_MEMORY) {
185 for (j = len0 - 2; j >= 0; j--) {
186 Ze[len1][j] = Ze[len1][j + 1] * e;
189 for (j = len0 - 2; j >= 0; j--) {
190 Ze[0][j] = Ze[0][j + 1] * e;
195 if (REVPART_FULL_MEMORY) {
198 Ze[len1][len0] = Zf[len1][len0] = 0;
199 Zf[len1 - 1][len0] = Zm[len1][len0] * endgapopen;
200 Ze[len1][len0 - 1] = Zm[len1][len0] * endgapopen;
203 for (i = len1 - 2; i >= 0; i--) {
204 Zf[i][len0] = Zf[i + 1][len0] * endgapextend;
210 for (j = len0 - 2; j >= 0; j--) {
211 Ze[len1][j] = Ze[len1][j + 1] * endgapextend;
217 // Zm(0) be the current row being filled/computed
218 // Zm(1) be the previous row
221 Ze[0][len0] = Zf[0][len0] = 0;
222 Zf[1][len0] = Zm[1][len0] * endgapopen;
223 Ze[0][len0 - 1] = Zm[1][len0] * endgapopen;
226 for (j = len0 - 2; j >= 0; j--) {
227 Ze[0][j] = Ze[0][j + 1] * endgapextend;
232 } //END FULL MEMORY and GAP enablement IF STATEMENT
234 double scorez, zz = 0;
236 for (i = len1 - 1; i >= 0; i--) {
238 for (j = len0 - 1; j >= 0; j--) {
239 Si = subst_index[sequences[1].text[i] - 'A'];
240 Tj = subst_index[sequences[0].text[j] - 'A'];
241 scorez = sub_matrix[Si][Tj];
243 //endgaps modification aug 10
244 double open0, extend0, open1, extend1;
247 extend0 = extend1 = e;
251 //check to see if one of the 2 sequences or both reach the end
255 extend0 = endgapextend;
261 extend1 = endgapextend;
266 if (REVPART_FULL_MEMORY) {
269 Ze[i][j] = Zm[i][j + 1] * open0 + Ze[i][j + 1] * extend0;
270 Zf[i][j] = Zm[i + 1][j] * open1 + Zf[i + 1][j] * extend1;
271 Zm[i][j] = (Zm[i + 1][j + 1] + Zf[i + 1][j + 1]
272 + Ze[i + 1][j + 1]) * scorez;
273 zz = Zm[i][j] + Zf[i][j] + Ze[i][j];
277 //2 ROW zE zF ALGORITHM GOES...:
278 //Ze[1][j] =Zm[i][j + 1] * exp(beta * open0) + Ze[1][j + 1] *exp(beta * extend0);
279 //Zf[1][j] = Zm[i + 1][j] * exp(beta * open1) + Zf[0][j] * exp(beta * extend1);
280 //Zm[i][j] = (Zm[i + 1][j + 1] + Zf[0][j + 1] + Ze[0][j + 1]) * exp(beta * scorez);
281 //zz = Zm[0][j] + Zf[1][j] + Ze[1][j];
283 //lowmem code for merging probability calculating module
284 //Here we make use of Zm as a 2 row matrix
286 Zf[1][j] = Zm[1][j] * open1 + Zf[0][j] * extend1;
287 Ze[1][j] = Zm[0][j + 1] * open0 + Ze[1][j + 1] * extend0;
288 Zm[0][j] = (Zm[1][j + 1] + Zf[0][j + 1] + Ze[0][j + 1])
291 tempvar = Zfm[i + 1][j + 1] * Zm[0][j];
292 //divide P(i,j) i.e. pairwise probability by denominator
293 tempvar /= (scorez * Zfm[0][0]);
294 probability = (float) tempvar;
296 //store only noticable probabilities
297 //if (probability <= 1 && probability >= 0.001) {
299 //validprob[i + 1][j + 1] = probability;
300 ptr[(j + 1) * (len1 + 1) + (i + 1)] = probability;
302 //lowmem code ends here
308 if (REVPART_FULL_MEMORY == 0) {
309 for (int t = 0; t <= sequences[0].length; t++) {
327 printf("\n\nrM:....\n\n");
328 if (REVPART_FULL_MEMORY) {
329 for (i = 0; i <= len1; i++) {
330 for (j = 0; j <= len0; j++)
331 printf("%.2Le ", Zm[i][j]);
335 printf("\n\nrE:....\n\n");
336 for (i = 0; i <= len1; i++) {
337 for (j = 0; j <= len0; j++)
338 printf("%.2Le ", Ze[i][j]);
343 printf("\n\nrF:....\n\n");
344 for (i = 0; i <= len1; i++) {
345 for (j = 0; j <= len0; j++)
346 printf("%.2Le ", Zf[i][j]);
360 //delete unused memory
362 if (REVPART_FULL_MEMORY) {
363 for (i = 0; i <= len1; i++) {
378 for (i = 0; i <= len1; i++) {
395 return (posteriorPtr);
399 //////////////////////////////////////////////////////////////
400 //forward partition function
401 /////////////////////////////////////////////////////////////
403 long double **partf(fasta sequences[2], const double termgapopen,
404 const double termgapextend, const double d, const double e) {
406 int i, j, len1, len0;
407 long double **Zm = NULL, **Zf = NULL, **Ze = NULL, zz = 0;
408 double endgapopen, endgapextend;
411 endgapopen = termgapopen;
412 endgapextend = termgapextend;
414 //the flag endgaps is set at the #define section
415 if (PART_FULL_MEMORY) {
417 Zf = new long double *[sequences[1].length + 1];
418 Ze = new long double *[sequences[1].length + 1];
419 Zm = new long double *[sequences[1].length + 1];
423 printf("\nPARTF:====\n");
425 //DYNAMICALLY GROW 2D M,IX,IY,PIX,PIY MARICES
426 for (i = 0; i <= sequences[1].length; i++) {
427 Zf[i] = new long double[sequences[0].length + 1];
428 Ze[i] = new long double[sequences[0].length + 1];
429 Zm[i] = new long double[sequences[0].length + 1];
432 Zm = new long double *[sequences[1].length + 1];
433 Ze = new long double *[2];
434 Zf = new long double *[2];
435 for (i = 0; i <= sequences[1].length; i++) {
436 Zm[i] = new long double[sequences[0].length + 1];
438 Ze[0] = new long double[sequences[0].length + 1];
439 Zf[0] = new long double[sequences[0].length + 1];
440 Ze[1] = new long double[sequences[0].length + 1];
441 Zf[1] = new long double[sequences[0].length + 1];
444 len0 = strlen(sequences[0].text);
445 len1 = strlen(sequences[1].text);
447 if (PART_FULL_MEMORY) {
448 for (i = 0; i <= sequences[1].length; i++)
449 for (j = 0; j <= sequences[0].length; j++) {
455 for (i = 0; i <= len1; i++) {
456 for (j = 0; j <= len0; j++) {
460 for (j = 0; j <= len0; j++) {
473 Zf[0][0] = Ze[0][0] = 0;
474 Zf[1][0] = Zm[0][0] * d;
475 Ze[0][1] = Zm[0][0] * d;
478 if (PART_FULL_MEMORY) {
479 for (i = 2; i <= sequences[1].length; i++) {
480 Zf[i][0] = Zf[i - 1][0] * e;
485 for (j = 2; j <= sequences[0].length; j++) {
486 Ze[0][j] = Ze[0][j - 1] * e;
491 Zf[0][0] = Ze[0][0] = 0;
492 Zf[1][0] = Zm[0][0] * endgapopen;
493 Ze[0][1] = Zm[0][0] * endgapopen;
496 if (PART_FULL_MEMORY) {
497 for (i = 2; i <= sequences[1].length; i++) {
498 Zf[i][0] = Zf[i - 1][0] * endgapextend;
503 for (j = 2; j <= sequences[0].length; j++) {
504 Ze[0][j] = Ze[0][j - 1] * endgapextend;
513 for (i = 1; i <= sequences[1].length; i++) {
515 for (j = 1; j <= sequences[0].length; j++) {
517 Si = subst_index[sequences[1].text[i - 1] - 'A'];
518 Tj = subst_index[sequences[0].text[j - 1] - 'A'];
520 score = sub_matrix[Si][Tj];
522 double open0, extend0, open1, extend1;
525 extend0 = extend1 = e;
528 //check to see if one of the 2 sequences or both reach the end
530 if (i == sequences[1].length) {
532 extend0 = endgapextend;
536 if (j == sequences[0].length) {
538 extend1 = endgapextend;
543 //z computation using open and extend temp vars
544 //open0 is gap open in seq0 and open1 is gap open in seq1
545 //entend0 is gap extend in seq0 and extend1 is gap extend in seq1
547 if (PART_FULL_MEMORY) {
548 Ze[i][j] = Zm[i][j - 1] * open0 + Ze[i][j - 1] * extend0;
550 if (Ze[i][j] >= OS_HUGE_VALL) {
551 printf("ERROR: huge val error for Ze\n");
555 Zf[i][j] = Zm[i - 1][j] * open1 + Zf[i - 1][j] * extend1;
557 if (Zf[i][j] >= OS_HUGE_VALL) {
558 printf("ERROR: huge val error for Zf\n");
562 Zm[i][j] = (Zm[i - 1][j - 1] + Ze[i - 1][j - 1]
563 + Zf[i - 1][j - 1]) * score;
565 if (Zm[i][j] >= OS_HUGE_VALL) {
566 printf("ERROR: huge val error for Zm\n");
570 zz = Zm[i][j] + Ze[i][j] + Zf[i][j];
572 Ze[1][j] = Zm[i][j - 1] * open0 + Ze[1][j - 1] * extend0;
574 if (Ze[1][j] >= OS_HUGE_VALL) {
575 printf("ERROR: huge val error for zE\n");
579 Zf[1][j] = Zm[i - 1][j] * open1 + Zf[0][j] * extend1;
581 if (Zf[1][j] >= OS_HUGE_VALL) {
582 printf("ERROR: huge val error for zF\n");
586 Zm[i][j] = (Zm[i - 1][j - 1] + Ze[0][j - 1] + Zf[0][j - 1])
589 if (Zm[i][j] >= OS_HUGE_VALL) {
590 printf("ERROR: huge val error for zM\n");
594 zz = Zm[i][j] + Ze[1][j] + Zf[1][j];
599 if (!PART_FULL_MEMORY) {
600 for (int t = 0; t <= sequences[0].length; t++) {
614 //store the sum of zm zf ze (m,n)s in zm's 0,0 th position
619 //print the 3 Z matrices namely Zm Zf and Ze
621 printf("\n\nFINAL Zm:\n");
622 for (i = 0; i <= sequences[1].length; i++) {
623 for (j = 0; j <= sequences[0].length; j++)
624 printf("%.2Le ", Zm[i][j]);
628 printf("FINAL Zf \n");
629 for (i = 0; i <= sequences[1].length; i++) {
630 for (j = 0; j <= sequences[0].length; j++)
631 printf("%.2Le ", Zf[i][j]);
635 printf("FINAL Ze \n");
636 for (i = 0; i <= sequences[1].length; i++) {
637 for (j = 0; j <= sequences[0].length; j++)
638 printf("%.2Le ", Ze[i][j]);
642 //end debug dump code
646 if (PART_FULL_MEMORY) {
647 for (i = 0; i <= sequences[1].length; i++) {
663 } //end of forward partition function
665 /////////////////////////////////////////////////////////////////////////////////////////
666 //entry point (was the main function) , returns the posterior probability safe vector
667 ////////////////////////////////////////////////////////////////////////////////////////
668 VF *ComputePostProbs(int a, int b, string seq1, string seq2) {
669 //printf("probamod\n");
670 double gap_open = -22, gap_ext = -1, beta = 0.2;//T = 5, beta = 1/T = 0.2, by default
673 double termgapopen = 1.0f; //exp(0)
674 double termgapextend = 1.0f; //exp(0)
676 //initialize the sequence structure
679 sequences[0].length = strlen((char *) seq1.c_str());
680 sequences[0].text = (char *) seq1.c_str();
681 sequences[0].title = new char[10];
682 strcpy(sequences[0].title, "seq0");
683 sequences[1].length = strlen((char *) seq2.c_str());
684 sequences[1].text = (char *) seq2.c_str();
685 sequences[1].title = new char[10];
686 strcpy(sequences[1].title, "seq1");
691 printf("%d %d %s\n%d %d %s\n--\n", a, sequences[0].length,
692 sequences[0].text, b, sequences[1].length, sequences[1].text);
693 printf("after init\n");
695 FILE *dump1 = fopen("dump1", "a");
696 fprintf(dump1, "%d %d %s\n%d %d %s\n--\n", a, sequences[0].length,
697 sequences[0].text, b, sequences[1].length, sequences[1].text);
701 gap_open = argument.gapopen;
702 gap_ext = argument.gapext;
703 beta = argument.beta;
705 stock_loop = argument.N;
706 le = argument.matrix;
708 //compute the values of exp(beta * ?)
709 termgapopen = exp(beta * 0.0);
710 termgapextend = exp(beta * 0.0);
711 gap_open = exp(beta * gap_open);
712 gap_ext = exp(beta * gap_ext);
715 printf("%f %f %f %d\n", gap_open, gap_ext, beta, le);
717 //call for calculating the posterior probabilities
718 // 1. call partition function partf
719 // 2. calculate revpartition using revers_parf
720 // 3. calculate probabilities
721 /// MODIFICATION... POPULATE SAFE VECTOR
725 MAT1 = partf(sequences, termgapopen, termgapextend, gap_open, gap_ext);
727 return revers_partf(sequences, termgapopen, termgapextend, MAT1, gap_open,
732 //////////////////////////////////////////////////////////////
733 //Compute Viterbi Alignment
734 // Added by YE Yongtao
735 /////////////////////////////////////////////////////////////
737 pair<SafeVector<char> *, float> partViterbi(string seq1, string seq2) {
740 double gap_open = -12, gap_ext = -1, beta = 0.2;//T = 5, beta = 1/T = 0.2, by default
743 //double termgapopen = 1.0f; //exp(0)
744 //double termgapextend = 1.0f; //exp(0)
746 //initialize the sequence structure
748 sequences[0].length = strlen((char *) seq1.c_str());
749 sequences[0].text = (char *) seq1.c_str();
750 sequences[0].title = new char[10];
751 strcpy(sequences[0].title, "seq0");
752 sequences[1].length = strlen((char *) seq2.c_str());
753 sequences[1].text = (char *) seq2.c_str();
754 sequences[1].title = new char[10];
755 strcpy(sequences[1].title, "seq1");
757 gap_open = argument.gapopen;
758 gap_ext = argument.gapext;
759 beta = argument.beta;
761 stock_loop = argument.N;
762 le = argument.matrix;
764 //compute the values of exp(beta * ?)
765 double endgapopen = exp(beta * 0.0);
766 double endgapextend = exp(beta * 0.0);
767 double d = exp(beta * gap_open);
768 double e = exp(beta * gap_ext);
770 int i, j, len1, len0;
771 long double **Zm = NULL, **Zf = NULL, **Ze = NULL;
772 int **traceZm = NULL, **traceZf = NULL, **traceZe = NULL;
774 //the flag endgaps is set at the #define section
775 Zf = new long double *[sequences[1].length + 1];
776 Ze = new long double *[sequences[1].length + 1];
777 Zm = new long double *[sequences[1].length + 1];
779 traceZf = new int *[sequences[1].length + 1];
780 traceZe = new int *[sequences[1].length + 1];
781 traceZm = new int *[sequences[1].length + 1];
783 //DYNAMICALLY GROW 2D M,IX,IY,PIX,PIY MARICES
784 for (i = 0; i <= sequences[1].length; i++) {
785 Zf[i] = new long double[sequences[0].length + 1];
786 Ze[i] = new long double[sequences[0].length + 1];
787 Zm[i] = new long double[sequences[0].length + 1];
789 traceZf[i] = new int[sequences[0].length + 1];
790 traceZe[i] = new int[sequences[0].length + 1];
791 traceZm[i] = new int[sequences[0].length + 1];
794 len0 = strlen(sequences[0].text);
795 len1 = strlen(sequences[1].text);
798 for (i = 0; i <= sequences[1].length; i++)
799 for (j = 0; j <= sequences[0].length; j++) {
814 Zf[0][0] = Ze[0][0] = 0;
815 Zf[1][0] = Zm[0][0] * d;
816 Ze[0][1] = Zm[0][0] * d;
820 for (i = 2; i <= sequences[1].length; i++) {
821 Zf[i][0] = Zf[i - 1][0] * e;
827 for (j = 2; j <= sequences[0].length; j++) {
828 Ze[0][j] = Ze[0][j - 1] * e;
834 Zf[0][0] = Ze[0][0] = 0;
835 Zf[1][0] = Zm[0][0] * endgapopen;
836 Ze[0][1] = Zm[0][0] * endgapopen;
840 for (i = 2; i <= sequences[1].length; i++) {
841 Zf[i][0] = Zf[i - 1][0] * endgapextend;
845 for (j = 2; j <= sequences[0].length; j++) {
846 Ze[0][j] = Ze[0][j - 1] * endgapextend;
856 for (i = 1; i <= sequences[1].length; i++) {
858 for (j = 1; j <= sequences[0].length; j++) {
860 Si = subst_index[sequences[1].text[i - 1] - 'A'];
861 Tj = subst_index[sequences[0].text[j - 1] - 'A'];
863 score = sub_matrix[Si][Tj];
865 double open0, extend0, open1, extend1;
868 extend0 = extend1 = e;
871 //check to see if one of the 2 sequences or both reach the end
873 if (i == sequences[1].length) {
875 extend0 = endgapextend;
879 if (j == sequences[0].length) {
881 extend1 = endgapextend;
886 //z computation using open and extend temp vars
887 //open0 is gap open in seq0 and open1 is gap open in seq1
888 //entend0 is gap extend in seq0 and extend1 is gap extend in seq1
889 Zf[i][j] = Zf[i - 1][j] * extend1;
892 if(Zm[i - 1][j] * open1 > Zf[i][j]){
893 Zf[i][j] = Zm[i - 1][j] * open1;
896 if (Zf[i][j] >= OS_HUGE_VALL) {
897 printf("ERROR: huge val error for Zf\n");
900 Ze[i][j] = Ze[i][j - 1] * extend0;
902 if(Zm[i][j - 1] * open0 > Ze[i][j]){
903 Ze[i][j] = Zm[i][j - 1] * open0;
907 if (Ze[i][j] >= OS_HUGE_VALL) {
908 printf("ERROR: huge val error for Ze\n");
912 Zm[i][j] = Zm[i - 1][j - 1] * score;
914 if(Zf[i - 1][j - 1] * score > Zm[i][j]){
915 Zm[i][j] = Zf[i - 1][j - 1] * score;
918 if(Ze[i - 1][j - 1] * score > Zm[i][j]){
919 Zm[i][j] = Ze[i - 1][j - 1] * score;
922 if (Zm[i][j] >= OS_HUGE_VALL) {
923 printf("ERROR: huge val error for Zm\n");
929 // figure out best terminating cell
931 float bestProb = Zm[sequences[1].length][sequences[0].length];
933 if( bestProb < Zf[sequences[1].length][sequences[0].length]){
934 bestProb = Zf[sequences[1].length][sequences[0].length];
937 if( bestProb < Ze[sequences[1].length][sequences[0].length]){
938 bestProb = Ze[sequences[1].length][sequences[0].length];
941 assert (state != -1);
944 SafeVector<char> *alignment = new SafeVector<char>; assert (alignment);
945 int c = sequences[1].length, r = sequences[0].length;
946 while (r != 0 || c != 0){
949 newState = traceZm[c][r];
950 c--; r--; alignment->push_back ('B');
953 newState = traceZe[c][r];
954 r--; alignment->push_back ('X');
957 newState = traceZf[c][r];
958 c--; alignment->push_back ('Y');
963 reverse (alignment->begin(), alignment->end());
965 for (i = 0; i <= sequences[1].length; i++) {
981 return make_pair(alignment, bestProb);
984 //////////////////////////////////////////////////////////////
985 // Compute two sequences' similarity defined as the normalized alignment score without gap penalties
986 // Added by YE Yongtao
987 /////////////////////////////////////////////////////////////
989 float computeSimilarity(string seq1, string seq2, SafeVector<char> * alignment) {
991 //initialize the sequence structure
993 sequences[0].length = strlen((char *) seq1.c_str());
994 sequences[0].text = (char *) seq1.c_str();
995 sequences[0].title = new char[10];
996 strcpy(sequences[0].title, "seq0");
997 sequences[1].length = strlen((char *) seq2.c_str());
998 sequences[1].text = (char *) seq2.c_str();
999 sequences[1].title = new char[10];
1000 strcpy(sequences[1].title, "seq1");
1005 int i = 1;int j = 1;
1006 for (SafeVector<char>::iterator iter = alignment->begin();
1007 iter != alignment->end(); ++iter){
1009 Si = subst_index[sequences[1].text[j - 1] - 'A'];
1010 Tj = subst_index[sequences[0].text[i - 1] - 'A'];
1011 score = normalized_matrix[Si][Tj];
1015 else if(*iter == 'X') i++;
1016 else if(*iter == 'Y') j++;
1018 if(i!= sequences[0].length + 1 || j!= sequences[1].length + 1 ) cerr << "similarity error"<< endl;
1019 bestProb /= alignment->size();
1020 //bestProb /= min(sequences[0].length, sequences[1].length);
1023 //end of posterior probability module