Mercurial > repos > clustalomega > clustalomega
view clustalomega/clustal-omega-0.2.0/src/hhalign/hhhit-C.h @ 0:ea6d0e588642 default tip
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
| author | clustalomega |
|---|---|
| date | Tue, 07 Jun 2011 16:13:02 -0400 |
| parents | |
| children |
line wrap: on
line source
/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */ /********************************************************************* * Clustal Omega - Multiple sequence alignment * * Copyright (C) 2010 University College Dublin * * Clustal-Omega is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation; either version 2 of the * License, or (at your option) any later version. * * This file is part of Clustal-Omega. * ********************************************************************/ /* * RCS $Id: hhhit-C.h 237 2011-04-14 15:01:48Z fabian $ */ // hhhit.C #ifndef MAIN #define MAIN #include <iostream> // cin, cout, cerr #include <fstream> // ofstream, ifstream #include <stdio.h> // printf using std::cout; using std::cerr; using std::endl; using std::ios; using std::ifstream; using std::ofstream; #include <stdlib.h> // exit #include <string> // strcmp, strstr #include <math.h> // sqrt, pow #include <limits.h> // INT_MIN #include <float.h> // FLT_MIN #include <time.h> // clock #include <ctype.h> // islower, isdigit etc #include "util-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc. #include "list.h" // list data structure #include "hash.h" // hash data structure #include "hhdecl-C.h" // constants, class #include "hhutil-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc. #include "hhhmm.h" // class HMM #include "hhalignment.h" // class Alignment #include "hhhitlist.h" // class HitList #endif #define CALCULATE_MAX6(max, var1, var2, var3, var4, var5, var6, varb) \ if (var1>var2) { max=var1; varb=STOP;} \ else { max=var2; varb=MM;}; \ if (var3>max) { max=var3; varb=GD;}; \ if (var4>max) { max=var4; varb=IM;}; \ if (var5>max) { max=var5; varb=DG;}; \ if (var6>max) { max=var6; varb=MI;}; #define CALCULATE_SUM6(sum, var1, var2, var3, var4, var5, var6, varb) \ if (var1>var2) { sum=var1; varb=STOP;} \ else { sum=var2; varb=MM;}; \ if (var3>sum) { sum=var3; varb=GD;}; \ if (var4>sum) { sum=var4; varb=IM;}; \ if (var5>sum) { sum=var5; varb=DG;}; \ if (var6>sum) { sum=var6; varb=MI;}; \ sum = var1 + var2 + var3 + var4 + var5 + var6; #define CALCULATE_MAX4(max, var1, var2, var3, var4, varb) \ if (var1>var2) { max=var1; varb=STOP;} \ else { max=var2; varb=MM;}; \ if (var3>max) { max=var3; varb=MI;}; \ if (var4>max) { max=var4; varb=IM;}; // Generate random number in [0,1[ #define frand() ((float) rand()/(RAND_MAX+1.0)) // Function declarations inline float Score(float* qi, float* tj); inline float ProbFwd(float* qi, float* tj); inline float max2(const float& xMM, const float& xX, char& b); inline int pickprob2(const double& xMM, const double& xX, const int& state); inline int pickprob3_GD(const double& xMM, const double& xDG, const double& xGD); inline int pickprob3_IM(const double& xMM, const double& xMI, const double& xIM); inline int pickprob6(const double& x0, const double& xMM, const double& xGD, const double& xIM, const double& xDG, const double& xMI); inline int pickmax2(const double& xMM, const double& xX, const int& state); inline int pickmax3_GD(const double& xMM, const double& xDG, const double& xGD); inline int pickmax3_IM(const double& xMM, const double& xMI, const double& xIM); inline int pickmax6(const double& x0, const double& xMM, const double& xGD, const double& xIM, const double& xDG, const double& xMI); inline double Pvalue(double x, double a[]); inline double Pvalue(float x, float lamda, float mu); inline double logPvalue(float x, float lamda, float mu); inline double logPvalue(float x, double a[]); inline double Probab(Hit& hit); ////////////////////////////////////////////////////////////////////////////// //// Constructor ////////////////////////////////////////////////////////////////////////////// Hit::Hit() { longname = name = file = dbfile = NULL; sname = NULL; seq = NULL; bMM = bGD = bDG = bIM = bMI = NULL; self = 0; i = j = NULL; states = NULL; S = S_ss = P_posterior = NULL; Xcons = NULL; B_MM=B_MI=B_IM=B_DG=B_GD=NULL; F_MM=F_MI=F_IM=F_DG=F_GD=NULL; cell_off = NULL; scale = NULL; sum_of_probs=0.0; Neff_HMM=0.0; score_ss = Pval = logPval = Eval = Probab = Pforward = 0.0; nss_conf = nfirst = i1 = i2 = j1 = j2 = matched_cols = ssm1 = ssm2 = 0; } ////////////////////////////////////////////////////////////////////////////// /** * @brief Free all allocated memory (to delete list of hits) */ void Hit::Delete() { if (i){ delete[] i; i = NULL; } if (j){ delete[] j; j = NULL; } if (states){ delete[] states; states = NULL; } if (S){ delete[] S; S = NULL; } if (S_ss){ delete[] S_ss; S_ss = NULL; } if (P_posterior){ delete[] P_posterior; P_posterior = NULL; } if (Xcons){ delete[] Xcons; Xcons = NULL; } // delete[] l; l = NULL; i = j = NULL; states = NULL; S = S_ss = P_posterior = NULL; Xcons = NULL; if (irep==1) // if irep>1 then longname etc point to the same memory locations as the first repeat. { // but these have already been deleted. // printf("Delete name = %s\n",name);////////////////// delete[] longname; longname = NULL; delete[] name; name = NULL; delete[] file; file = NULL; delete[] dbfile; dbfile = NULL; for (int k=0; k<n_display; k++) { //delete[] sname[k]; sname[k] = NULL; delete[] seq[k]; seq[k] = NULL; } //delete[] sname; sname = NULL; delete[] seq; seq = NULL; } } /////////////////////////////////////////////////////////////////////////// /** * @brief Allocate/delete memory for dynamic programming matrix */ void Hit::AllocateBacktraceMatrix(int Nq, int Nt) { int i; bMM=new(char*[Nq]); bMI=new(char*[Nq]); bIM=new(char*[Nq]); bDG=new(char*[Nq]); bGD=new(char*[Nq]); cell_off=new(char*[Nq]); for (i=0; i<Nq; i++) { bMM[i]=new(char[Nt]); bMI[i]=new(char[Nt]); bIM[i]=new(char[Nt]); bGD[i]=new(char[Nt]); bDG[i]=new(char[Nt]); cell_off[i]=new(char[Nt]); if (!bMM[i] || !bMI[i] || !bIM[i] || !bGD[i] || !bDG[i] || !cell_off[i]) { fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq); fprintf(stderr,"Suggestions:\n"); fprintf(stderr,"1. Cut query sequence into shorter segments\n"); fprintf(stderr,"2. Check stack size limit (Linux: ulimit -a)\n"); fprintf(stderr,"3. Run on a computer with bigger memory\n"); exit(3); } } } /** * @brief */ void Hit::DeleteBacktraceMatrix(int Nq) { int i; for (i=0; i<Nq; i++) { delete[] bMM[i]; bMM[i] = NULL; delete[] bMI[i]; bMI[i] = NULL; delete[] bIM[i]; bIM[i] = NULL; delete[] bGD[i]; bGD[i] = NULL; delete[] bDG[i]; bDG[i] = NULL; delete[] cell_off[i]; cell_off[i] = NULL; } delete[] bMM; bMM = NULL; delete[] bMI; bMI = NULL; delete[] bIM; bIM = NULL; delete[] bDG; bDG = NULL; delete[] bGD; bGD = NULL; delete[] cell_off; cell_off = NULL; } /////////////////////////////////////////////////////////////////////////////// /** * @brief Allocate/delete memory for Forward dynamic programming matrix */ void Hit::AllocateForwardMatrix(int Nq, int Nt) { F_MM=new(double*[Nq]); F_MI=new(double*[Nq]); F_DG=new(double*[Nq]); F_IM=new(double*[Nq]); F_GD=new(double*[Nq]); scale=new(double[Nq+1]); // need Nq+3? for (int i=0; i<Nq; i++) { F_MM[i] = new(double[Nt]); F_MI[i] = new(double[Nt]); F_DG[i] = new(double[Nt]); F_IM[i] = new(double[Nt]); F_GD[i] = new(double[Nt]); if (!F_MM[i] || !F_MI[i] || !F_IM[i] || !F_GD[i] || !F_DG[i]) { fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq); fprintf(stderr,"Suggestions:\n"); fprintf(stderr,"1. Cut query sequence into shorter segments\n"); fprintf(stderr,"2. Check stack size limit (Linux: ulimit -a)\n"); fprintf(stderr,"3. Run on a computer with bigger memory\n"); exit(3); } } } /** * @brief */ void Hit::DeleteForwardMatrix(int Nq) { for (int i=0; i<Nq; i++) { delete[] F_MM[i]; F_MM[i] = NULL; delete[] F_MI[i]; F_MI[i] = NULL; delete[] F_IM[i]; F_IM[i] = NULL; delete[] F_GD[i]; F_GD[i] = NULL; delete[] F_DG[i]; F_DG[i] = NULL; } delete[] F_MM; F_MM = NULL; delete[] F_MI; F_MI = NULL; delete[] F_IM; F_IM = NULL; delete[] F_DG; F_DG = NULL; delete[] F_GD; F_GD = NULL; delete[] scale; scale = NULL; } ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Allocate/delete memory for Backward dynamic programming matrix (DO ONLY AFTER FORWARD MATRIX HAS BEEN ALLOCATED) */ void Hit::AllocateBackwardMatrix(int Nq, int Nt) { B_MM=new(double*[Nq]); B_MI=F_MI; B_DG=F_DG; B_IM=F_IM; B_GD=F_GD; for (int i=0; i<Nq; i++) { B_MM[i] = new(double[Nt]); if (!B_MM[i]) { fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq); fprintf(stderr,"Suggestions:\n"); fprintf(stderr,"1. Cut query sequence into shorter segments\n"); fprintf(stderr,"2. Check stack size limit (Linux: ulimit -a)\n"); fprintf(stderr,"3. Run on a computer with bigger memory\n"); exit(3); } } } void Hit::DeleteBackwardMatrix(int Nq) { for (int i=0; i<Nq; i++) { delete[] B_MM[i]; B_MM[i] = NULL; /* is this all? FS */ } delete[] B_MM; B_MM = NULL; B_MM=B_MI=B_IM=B_DG=B_GD=NULL; } ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Compare HMMs with one another and look for sub-optimal alignments that share no pair with previous ones * The function is called with q and t * If q and t are equal (self==1), only the upper right part of the matrix is calculated: j>=i+3 */ void Hit::Viterbi(HMM& q, HMM& t, float** Sstruc) { // Linear topology of query (and template) HMM: // 1. The HMM HMM has L+2 columns. Columns 1 to L contain // a match state, a delete state and an insert state each. // 2. The Start state is M0, the virtual match state in column i=0 (j=0). (Therefore X[k][0]=ANY) // This column has only a match state and it has only a transitions to the next match state. // 3. The End state is M(L+1), the virtual match state in column i=L+1.(j=L+1) (Therefore X[k][L+1]=ANY) // Column L has no transitions to the delete state: tr[L][M2D]=tr[L][D2D]=0. // 4. Transitions I->D and D->I are ignored, since they do not appear in PsiBlast alignments // (as long as the gap opening penalty d is higher than the best match score S(a,b)). // Pairwise alignment of two HMMs: // 1. Pair-states for the alignment of two HMMs are // MM (Q:Match T:Match) , GD (Q:Gap T:Delete), IM (Q:Insert T:Match), DG (Q:Delelte, T:Match) , MI (Q:Match T:Insert) // 2. Transitions are allowed only between the MM-state and each of the four other states. // Saving space: // The best score ending in pair state XY sXY[i][j] is calculated from left to right (j=1->t.L) // and top to bottom (i=1->q.L). To save space, only the last row of scores calculated is kept in memory. // (The backtracing matrices are kept entirely in memory [O(t.L*q.L)]). // When the calculation has proceeded up to the point where the scores for cell (i,j) are caculated, // sXY[i-1][j'] = sXY[j'] for j'>=j (A below) // sXY[i][j'] = sXY[j'] for j'<j (B below) // sXY[i-1][j-1]= sXY_i_1_j_1 (C below) // sXY[i][j] = sXY_i_j (D below) // j-1 // j // i-1: CAAAAAAAAAAAAAAAAAA // i : BBBBBBBBBBBBBD // Variable declarations //float sMM[MAXRES]; // sMM[i][j] = score of best alignment up to indices (i,j) ending in (Match,Match) //float sGD[MAXRES]; // sGD[i][j] = score of best alignment up to indices (i,j) ending in (Gap,Delete) //float sDG[MAXRES]; // sDG[i][j] = score of best alignment up to indices (i,j) ending in (Delete,Gap) //float sIM[MAXRES]; // sIM[i][j] = score of best alignment up to indices (i,j) ending in (Ins,Match) //float sMI[MAXRES]; // sMI[i][j] = score of best alignment up to indices (i,j) ending in (Match,Ins) float *sMM = new(float[par.maxResLen]); // sMM[i][j] = score of best alignment up to indices (i,j) ending in (Match,Match) float *sGD = new(float[par.maxResLen]); // sGD[i][j] = score of best alignment up to indices (i,j) ending in (Gap,Delete) float *sDG = new(float[par.maxResLen]); // sDG[i][j] = score of best alignment up to indices (i,j) ending in (Delete,Gap) float *sIM = new(float[par.maxResLen]); // sIM[i][j] = score of best alignment up to indices (i,j) ending in (Ins,Match) float *sMI = new(float[par.maxResLen]); // sMI[i][j] = score of best alignment up to indices (i,j) ending in (Match,Ins) float smin=(par.loc? 0:-FLT_MAX); //used to distinguish between SW and NW algorithms in maximization int i,j; //query and template match state indices float sMM_i_j=0,sMI_i_j,sIM_i_j,sGD_i_j,sDG_i_j; float sMM_i_1_j_1,sMI_i_1_j_1,sIM_i_1_j_1,sGD_i_1_j_1,sDG_i_1_j_1; int jmin,jmax; // Reset crossed out cells? if(irep==1) InitializeForAlignment(q,t); // Initialization of top row, i.e. cells (0,j) for (j=0; j<=t.L; j++) { sMM[j] = (self? 0 : -j*par.egt); sIM[j] = sMI[j] = sDG[j] = sGD[j] = -FLT_MAX; } score=-INT_MAX; i2=j2=0; bMM[0][0]=STOP; // Viterbi algorithm for (i=1; i<=q.L; i++) // Loop through query positions i { // if (v>=5) printf("\n"); if (self) { // If q is compared to itself, ignore cells below diagonal+SELFEXCL jmin = i+SELFEXCL; jmax = t.L; if (jmin>jmax) continue; } else { // If q is compared to t, exclude regions where overlap of q with t < min_overlap residues jmin=imax( 1, i+min_overlap-q.L); // Lq-i+j>=Ovlap => j>=i+Ovlap-Lq => jmin=max{1, i+Ovlap-Lq} jmax=imin(t.L,i-min_overlap+t.L); // Lt-j+i>=Ovlap => j<=i-Ovlap+Lt => jmax=min{Lt,i-Ovlap+Lt} } // Initialize cells if (jmin==1) { sMM_i_1_j_1 = -(i-1)*par.egq; // initialize at (i-1,0) sMM[0] = -i*par.egq; // initialize at (i,0) sIM_i_1_j_1 = sMI_i_1_j_1 = sDG_i_1_j_1 = sGD_i_1_j_1 = -FLT_MAX; // initialize at (i-1,jmin-1) } else { // Initialize at (i-1,jmin-1) if lower left triagonal is excluded due to min overlap sMM_i_1_j_1 = sMM[jmin-1]; // initialize at (i-1,jmin-1) sIM_i_1_j_1 = sIM[jmin-1]; // initialize at (i-1,jmin-1) sMI_i_1_j_1 = sMI[jmin-1]; // initialize at (i-1,jmin-1) sDG_i_1_j_1 = sDG[jmin-1]; // initialize at (i-1,jmin-1) sGD_i_1_j_1 = sGD[jmin-1]; // initialize at (i-1,jmin-1) sMM[jmin-1] = -FLT_MAX; // initialize at (i,jmin-1) } if (jmax<t.L) // initialize at (i-1,jmmax) if upper right triagonal is excluded due to min overlap sMM[jmax] = sIM[jmax] = sMI[jmax] = sDG[jmax] = sGD[jmax] = -FLT_MAX; sIM[jmin-1] = sMI[jmin-1] = sDG[jmin-1] = sGD[jmin-1] = -FLT_MAX; // initialize at (i,jmin-1) for (j=jmin; j<=jmax; j++) // Loop through template positions j { if (cell_off[i][j]) { sMM_i_1_j_1 = sMM[j]; // sMM_i_1_j_1 (for j->j+1) = sMM(i-1,(j+1)-1) = sMM[j] sGD_i_1_j_1 = sGD[j]; sIM_i_1_j_1 = sIM[j]; sDG_i_1_j_1 = sDG[j]; sMI_i_1_j_1 = sMI[j]; sMM[j]=sMI[j]=sIM[j]=sDG[j]=sGD[j]=-FLT_MAX; // sMM[j] = sMM(i,j) is cell_off } else { // Recursion relations // printf("S[%i][%i]=%4.1f ",i,j,Score(q.p[i],t.p[j])); // DEBUG!! CALCULATE_MAX6( sMM_i_j, smin, sMM_i_1_j_1 + q.tr[i-1][M2M] + t.tr[j-1][M2M], sGD_i_1_j_1 + q.tr[i-1][M2M] + t.tr[j-1][D2M], sIM_i_1_j_1 + q.tr[i-1][I2M] + t.tr[j-1][M2M], sDG_i_1_j_1 + q.tr[i-1][D2M] + t.tr[j-1][M2M], sMI_i_1_j_1 + q.tr[i-1][M2M] + t.tr[j-1][I2M], bMM[i][j] ); sMM_i_j += Score(q.p[i],t.p[j]) + ScoreSS(q,t,i,j) + par.shift + (Sstruc==NULL? 0: Sstruc[i][j]); sGD_i_j = max2 ( sMM[j-1] + t.tr[j-1][M2D], // MM->GD gap opening in query sGD[j-1] + t.tr[j-1][D2D], // GD->GD gap extension in query bGD[i][j] ); sIM_i_j = max2 ( // sMM[j-1] + q.tr[i][M2I] + t.tr[j-1][M2M] , sMM[j-1] + q.tr[i][M2I] + t.tr[j-1][M2M_GAPOPEN], // MM->IM gap opening in query sIM[j-1] + q.tr[i][I2I] + t.tr[j-1][M2M], // IM->IM gap extension in query bIM[i][j] ); sDG_i_j = max2 ( // sMM[j] + q.tr[i-1][M2D], // sDG[j] + q.tr[i-1][D2D], //gap extension (DD) in query sMM[j] + q.tr[i-1][M2D] + t.tr[j][GAPOPEN], // MM->DG gap opening in template sDG[j] + q.tr[i-1][D2D] + t.tr[j][GAPEXTD], // DG->DG gap extension in template bDG[i][j] ); sMI_i_j = max2 ( sMM[j] + q.tr[i-1][M2M] + t.tr[j][M2I], // MM->MI gap opening M2I in template sMI[j] + q.tr[i-1][M2M] + t.tr[j][I2I], // MI->MI gap extension I2I in template bMI[i][j] ); sMM_i_1_j_1 = sMM[j]; sGD_i_1_j_1 = sGD[j]; sIM_i_1_j_1 = sIM[j]; sDG_i_1_j_1 = sDG[j]; sMI_i_1_j_1 = sMI[j]; sMM[j] = sMM_i_j; sGD[j] = sGD_i_j; sIM[j] = sIM_i_j; sDG[j] = sDG_i_j; sMI[j] = sMI_i_j; //if (isnan(sMM_i_j)||isinf(sMM_i_j)){ // printf("."); /* <DEBUG> FS*/ //} // Find maximum score; global alignment: maxize only over last row and last column if(sMM_i_j>score && (par.loc || i==q.L)) { i2=i; j2=j; score=sMM_i_j; } } // end if } //end for j // if global alignment: look for best cell in last column if (!par.loc && sMM_i_j>score) { i2=i; j2=jmax; score=sMM_i_j; } } // end for i state=MM; // state with maximum score is MM state // If local alignment do length correction: -log(length) if (par.loc) { if (self) score=score-log(0.5*t.L*q.L/200.0/200.0)/LAMDA - 11.2; // offset of -11.2 to get approx same mean as for -global else if (par.idummy==0 && q.lamda>0) ////////////////////////////////////////////// score=score-log(t.L*q.L/200.0/200.0)/q.lamda - 11.2; // offset of -11.2 to get approx same mean as for -global else if (par.idummy<=1) ////////////////////////////////////////////// score=score-log(t.L*q.L/200.0/200.0)/LAMDA - 11.2; // offset of -11.2 to get approx same mean as for -global } // printf("Template=%-12.12s i=%-4i j=%-4i score=%6.3f\n",t.name,i2,j2,score); delete[] sMM; sMM = NULL; delete[] sGD; sGD = NULL; delete[] sDG; sDG = NULL; delete[] sIM; sIM = NULL; delete[] sMI; sMI = NULL; return; } /* this is the end of Hit::Viterbi() */ ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Compare two HMMs with Forward Algorithm in lin-space (~ 2x faster than in log-space) */ void Hit::Forward(HMM& q, HMM& t, float** Pstruc) { // Variable declarations int i,j; // query and template match state indices double pmin=(par.loc? 1.0: 0.0); // used to distinguish between SW and NW algorithms in maximization double Cshift = pow(2.0,par.shift); // score offset transformed into factor in lin-space double Pmax_i; // maximum of F_MM in row i double scale_prod=1.0; // Prod_i=1^i (scale[i]) int jmin; // First alignment of this pair of HMMs? if(irep==1) { q.tr[0][M2D] = q.tr[0][M2I] = 0.0; q.tr[0][I2M] = q.tr[0][I2I] = 0.0; q.tr[0][D2M] = q.tr[0][D2D] = 0.0; t.tr[0][M2M] = 1.0; t.tr[0][M2D] = t.tr[0][M2I] = 0.0; t.tr[0][I2M] = t.tr[0][I2I] = 0.0; t.tr[0][D2M] = t.tr[0][D2D] = 0.0; q.tr[q.L][M2M] = 1.0; q.tr[q.L][M2D] = q.tr[q.L][M2I] = 0.0; q.tr[q.L][I2M] = q.tr[q.L][I2I] = 0.0; q.tr[q.L][D2M] = 1.0; q.tr[q.L][D2D] = 0.0; t.tr[t.L][M2M] = 1.0; t.tr[t.L][M2D] = t.tr[t.L][M2I] = 0.0; t.tr[t.L][I2M] = t.tr[t.L][I2I] = 0.0; t.tr[t.L][D2M] = 1.0; t.tr[t.L][D2D] = 0.0; InitializeForAlignment(q,t); } // Initialization of top row, i.e. cells (0,j) F_MM[1][0] = F_IM[1][0] = F_GD[1][0] = F_MM[0][1] = F_MI[0][1] = F_DG[0][1] = 0.0; for (j=1; j<=t.L; j++) { if (cell_off[1][j]) F_MM[1][j] = F_MI[1][j] = F_DG[1][j] = F_IM[1][j] = F_GD[1][j] = 0.0; else { F_MM[1][j] = ProbFwd(q.p[1],t.p[j]) * fpow2(ScoreSS(q,t,1,j)) * Cshift * (Pstruc==NULL? 1: Pstruc[1][j]) ; F_MI[1][j] = F_DG[1][j] = 0.0; F_IM[1][j] = F_MM[1][j-1] * q.tr[1][M2I] * t.tr[j-1][M2M] + F_IM[1][j-1] * q.tr[1][I2I] * t.tr[j-1][M2M]; F_GD[1][j] = F_MM[1][j-1] * t.tr[j-1][M2D] + F_GD[1][j-1] * t.tr[j-1][D2D]; } } scale[0]=scale[1]=scale[2]=1.0; // Forward algorithm for (i=2; i<=q.L; i++) // Loop through query positions i { // if (v>=5) printf("\n"); if (self) jmin = imin(i+SELFEXCL+1,t.L); else jmin=1; if (scale_prod<DBL_MIN*100) scale_prod = 0.0; else scale_prod *= scale[i]; // Initialize cells at (i,0) if (cell_off[i][jmin]) F_MM[i][jmin] = F_MI[i][jmin] = F_DG[i][jmin] = F_IM[i][jmin] = F_GD[i][jmin] = 0.0; else { F_MM[i][jmin] = scale_prod * ProbFwd(q.p[i],t.p[jmin]) * fpow2(ScoreSS(q,t,i,jmin)) * Cshift * (Pstruc==NULL? 1: Pstruc[i][jmin]); F_IM[i][jmin] = F_GD[i][jmin] = 0.0; F_MI[i][jmin] = scale[i] * (F_MM[i-1][jmin] * q.tr[i-1][M2M] * t.tr[jmin][M2I] + F_MI[i-1][jmin] * q.tr[i-1][M2M] * t.tr[jmin][I2I]); F_DG[i][jmin] = scale[i] * (F_MM[i-1][jmin] * q.tr[i-1][M2D] + F_DG[i-1][jmin] * q.tr[i-1][D2D]); } Pmax_i=0; for (j=jmin+1; j<=t.L; j++) // Loop through template positions j { // Recursion relations // printf("S[%i][%i]=%4.1f ",i,j,Score(q.p[i],t.p[j])); if (cell_off[i][j]) F_MM[i][j] = F_MI[i][j] = F_DG[i][j] = F_IM[i][j] = F_GD[i][j] = 0.0; else { F_MM[i][j] = ProbFwd(q.p[i],t.p[j]) * fpow2(ScoreSS(q,t,i,j)) * Cshift * (Pstruc==NULL? 1: Pstruc[i][j]) * scale[i] * ( pmin + F_MM[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][M2M] // BB -> MM (BB = Begin/Begin, for local alignment) + F_GD[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][D2M] // GD -> MM + F_IM[i-1][j-1] * q.tr[i-1][I2M] * t.tr[j-1][M2M] // IM -> MM + F_DG[i-1][j-1] * q.tr[i-1][D2M] * t.tr[j-1][M2M] // DG -> MM + F_MI[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][I2M] // MI -> MM ); F_GD[i][j] = ( F_MM[i][j-1] * t.tr[j-1][M2D] // GD -> MM + F_GD[i][j-1] * t.tr[j-1][D2D] // GD -> GD + (Pstruc==NULL? 0 : F_DG[i][j-1] * t.tr[j-1][M2D] * q.tr[i][D2M] ) // DG -> GD (only when structure scores given) ); F_IM[i][j] = ( F_MM[i][j-1] * q.tr[i][M2I] * t.tr[j-1][M2M] // MM -> IM + F_IM[i][j-1] * q.tr[i][I2I] * t.tr[j-1][M2M] // IM -> IM + (Pstruc==NULL? 0 : F_MI[i][j-1] * q.tr[i][M2I] * t.tr[j-1][I2M] ) // MI -> IM (only when structure scores given) ); F_DG[i][j] = scale[i] * ( F_MM[i-1][j] * q.tr[i-1][M2D] // DG -> MM + F_DG[i-1][j] * q.tr[i-1][D2D] // DG -> DG ) ; F_MI[i][j] = scale[i] * ( F_MM[i-1][j] * q.tr[i-1][M2M] * t.tr[j][M2I] // MI -> MM + F_MI[i-1][j] * q.tr[i-1][M2M] * t.tr[j][I2I] // MI -> MI ); if(F_MM[i][j]>Pmax_i) Pmax_i=F_MM[i][j]; } // end else } //end for j pmin *= scale[i]; scale[i+1] = 1.0/(Pmax_i+1.0); // scale[i+1] = 1.0; } // end for i // Calculate P_forward * Product_{i=1}^{Lq+1}(scale[i]) if (par.loc) { Pforward = 1.0; // alignment contains no residues (see Mueckstein, Stadler et al.) for (i=1; i<=q.L; i++) // Loop through query positions i { for (j=1; j<=t.L; j++) // Loop through template positions j Pforward += F_MM[i][j]; Pforward *= scale[i+1]; } } else // global alignment { Pforward = 0.0; for (i=1; i<q.L; i++) { Pforward = (Pforward + F_MM[i][t.L]) * scale[i+1]; } for (j=1; j<=t.L; j++) { Pforward += F_MM[q.L][j]; } Pforward *= scale[q.L+1]; } // Calculate log2(P_forward) score = log2(Pforward)-10.0f; for (i=1; i<=q.L+1; i++) score -= log2(scale[i]); // state = MM; if (par.loc) { if (self) score=score-log(0.5*t.L*q.L)/LAMDA+14.; // +14.0 to get approx same mean as for -global else score=score-log(t.L*q.L)/LAMDA+14.; // +14.0 to get approx same mean as for -global } // Debugging output if (v>=6) { const int i0=0, i1=q.L; const int j0=0, j1=t.L; scale_prod=1; printf("\nFwd scale "); for (j=j0; j<=j1; j++) printf("%3i ",j); printf("\n"); for (i=i0; i<=i1; i++) { scale_prod *= scale[i]; printf("%3i: %9.3G ",i,1/scale_prod); for (j=j0; j<=j1; j++) printf("%7.4f ",(F_MM[i][j]+F_MI[i][j]+F_IM[i][j]+F_DG[i][j]+F_GD[i][j])); printf("\n"); // printf(" MM %9.5f ",1/scale[i]); // for (j=j0; j<=j1; j++) // printf("%7.4f ",F_MM[i][j]); // printf("\n"); } } // printf("Template=%-12.12s score=%6.3f i2=%i j2=%i \n",t.name,score,i2,j2); return; } /* this is the end of Hit::Forward() */ ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Compare two HMMs with Backward Algorithm (in lin-space, 2x faster), for use in MAC alignment */ void Hit::Backward(HMM& q, HMM& t) { // Variable declarations int i,j; // query and template match state indices double pmin=(par.loc? 1.0: 0.0); // used to distinguish between SW and NW algorithms in maximization double Cshift = pow(2.0,par.shift); // score offset transformed into factor in lin-space double scale_prod=scale[q.L+1]; int jmin; //double dMaxB = -1.0; // Initialization of top row, i.e. cells (0,j) for (j=t.L; j>=1; j--) { if (cell_off[q.L][j]) B_MM[q.L][j] = 0.0; else B_MM[q.L][j] = scale[q.L+1]; //dMaxB = dMaxB>B_MM[q.L][j]?dMaxB:B_MM[q.L][j]; B_IM[q.L][j] = B_MI[q.L][j] = B_DG[q.L][j] = B_GD[q.L][j] = 0.0; } if (par.loc) pmin = scale[q.L+1]; // transform pmin (for local alignment) to scale of present (i'th) row // Backward algorithm for (i=q.L-1; i>=1; i--) // Loop through query positions i { // if (v>=5) printf("\n"); if (self) jmin = imin(i+SELFEXCL,t.L); else jmin=1; // jmin = i+SELFEXCL and not (i+SELFEXCL+1) to set matrix element at boundary to zero // Initialize cells at (i,t.L+1) scale_prod *= scale[i+1]; if (cell_off[i][t.L]) B_MM[i][t.L] = 0.0; else B_MM[i][t.L] = scale_prod; //if (isnan(B_MM[i][t.L])||isinf(B_MM[i][t.L])){ // printf("."); /* <DEBUG> FS*/ //} //dMaxB = dMaxB>B_MM[i][t.L]?dMaxB:B_MM[i][t.L]; B_IM[i][t.L] = B_MI[i][t.L] = B_DG[i][t.L] = B_GD[i][t.L] = 0.0; pmin *= scale[i+1]; // transform pmin (for local alignment) to scale of present (i'th) row for (j=t.L-1; j>=jmin; j--) // Loop through template positions j { // Recursion relations // printf("S[%i][%i]=%4.1f ",i,j,Score(q.p[i],t.p[j])); if (cell_off[i][j]) B_MM[i][j] = B_GD[i][j] = B_IM[i][j] = B_DG[i][j] = B_MI[i][j] = 0.0; else { double pmatch = B_MM[i+1][j+1] * ProbFwd(q.p[i+1],t.p[j+1]) * fpow2(ScoreSS(q,t,i+1,j+1)) * Cshift * scale[i+1]; //if (isnan(pmatch)||isinf(pmatch)){ // printf("."); /* <DEBUG> FS*/ //} B_MM[i][j] = ( + pmin // MM -> EE (End/End, for local alignment) + pmatch * q.tr[i][M2M] * t.tr[j][M2M] // MM -> MM + B_GD[i][j+1] * t.tr[j][M2D] // MM -> GD (q.tr[i][M2M] is already contained in GD->MM) + B_IM[i][j+1] * q.tr[i][M2I] * t.tr[j][M2M] // MM -> IM + B_DG[i+1][j] * q.tr[i][M2D] * scale[i+1] // MM -> DG (t.tr[j][M2M] is already contained in DG->MM) + B_MI[i+1][j] * q.tr[i][M2M] * t.tr[j][M2I] * scale[i+1] // MM -> MI ); //if (isnan(B_MM[i][j])||isinf(B_MM[i][j])){ // printf("."); /* <DEBUG> FS*/ //} //dMaxB = dMaxB>B_MM[i][j]?dMaxB:B_MM[i][j]; B_GD[i][j] = ( + pmatch * q.tr[i][M2M] * t.tr[j][D2M] // GD -> MM + B_GD[i][j+1] * t.tr[j][D2D] // DG -> DG ); B_IM[i][j] = ( + pmatch * q.tr[i][I2M] * t.tr[j][M2M] // IM -> MM + B_IM[i][j+1] * q.tr[i][I2I] * t.tr[j][M2M] // IM -> IM ); B_DG[i][j] = ( + pmatch * q.tr[i][D2M] * t.tr[j][M2M] // DG -> MM + B_DG[i+1][j] * q.tr[i][D2D] * scale[i+1] // DG -> DG // + B_GD[i][j+1] * q.tr[i][D2M] * t.tr[j][M2D] // DG -> GD ); B_MI[i][j] = ( + pmatch * q.tr[i][M2M] * t.tr[j][I2M] // MI -> MM + B_MI[i+1][j] * q.tr[i][M2M] * t.tr[j][I2I] * scale[i+1] // MI -> MI // + B_IM[i][j+1] * q.tr[i][M2I] * t.tr[j][I2M] // MI -> IM ); } // end else } //end for j } // end for i // Debugging output if (v>=6) { const int i0=0, i1=q.L; const int j0=0, j1=t.L; double scale_prod[q.L+2]; scale_prod[q.L] = scale[q.L+1]; for (i=q.L-1; i>=1; i--) scale_prod[i] = scale_prod[i+1] * scale[i+1]; printf("\nBwd scale "); for (j=j0; j<=j1; j++) printf("%3i ",j); printf("\n"); for (i=i0; i<=i1; i++) { printf("%3i: %9.3G ",i,1/scale_prod[i]); for (j=j0; j<=j1; j++) printf("%7.4f ",(B_MM[i][j]+B_MI[i][j]+B_IM[i][j]+B_DG[i][j]+B_GD[i][j]) * (ProbFwd(q.p[i],t.p[j])*fpow2(ScoreSS(q,t,i,j)) * Cshift)); printf("\n"); // printf("MM %9.5f ",1/scale[i]); // for (j=j0; j<=j1; j++) // printf("%7.4f ",B_MM[i][j] * (ProbFwd(q.p[i],t.p[j])*fpow2(ScoreSS(q,t,i,j)) * Cshift)); // printf("\n"); } printf("\nPost scale "); for (j=j0; j<=j1; j++) printf("%3i ",j); printf("\n"); for (i=i0; i<=i1; i++) { printf("%3i: %9.3G ",i,1/scale_prod[i]); for (j=j0; j<=j1; j++) printf("%7.4f ",B_MM[i][j]*F_MM[i][j]/Pforward); printf("\n"); } printf("\n"); } if (v>=4) printf("\nForward total probability ratio: %8.3G\n",Pforward); // Calculate Posterior matrix and overwrite Backward matrix with it for (i=1; i<=q.L; i++) { for (j=1; j<=t.L; j++) { B_MM[i][j] *= F_MM[i][j]/Pforward; //if (isnan(B_MM[i][j]) || isinf(B_MM[i][j])){ // printf("."); /* <DEBUG> FS*/ //} //dMaxB = dMaxB>B_MM[i][j]?dMaxB:B_MM[i][j]; } } //printf("Max-B_MM = %f\n", dMaxB); return; } /* this is the end of Hit::Backward() */ ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Maximum Accuracy alignment */ void Hit::MACAlignment(HMM& q, HMM& t) { // Use Forward and Backward matrices to find that alignment which // maximizes the expected number of correctly aligned pairs of residues (mact=0) // or, more generally, which maximizes the expectation value of the number of // correctly aligned pairs minus (mact x number of aligned pairs) // "Correctly aligned" can be based on posterior probabilities calculated with // a local or a global version of the Forward-Backward algorithm. int i,j; // query and template match state indices int jmin,jmax; // range of dynamic programming for j double** S=F_MI; // define alias for new score matrix double score_MAC; // score of the best MAC alignment // Initialization of top row, i.e. cells (0,j) for (j=0; j<=t.L; j++) S[0][j] = 0.0; score_MAC=-INT_MAX; i2=j2=0; bMM[0][0]=STOP; // Dynamic programming for (i=1; i<=q.L; i++) // Loop through query positions i { if (self) { // If q is compared to itself, ignore cells below diagonal+SELFEXCL jmin = i+SELFEXCL; jmax = t.L; if (jmin>jmax) continue; } else { // If q is compared to t, exclude regions where overlap of q with t < min_overlap residues jmin=imax( 1, i+min_overlap-q.L); // Lq-i+j>=Ovlap => j>=i+Ovlap-Lq => jmin=max{1, i+Ovlap-Lq} jmax=imin(t.L,i-min_overlap+t.L); // Lt-j+i>=Ovlap => j<=i-Ovlap+Lt => jmax=min{Lt,i-Ovlap+Lt} } // Initialize cells S[i][jmin-1] = 0.0; if (jmax<t.L) S[i-1][jmax] = 0.0; // initialize at (i-1,jmax) if upper right triagonal is excluded due to min overlap for (j=jmin; j<=jmax; j++) // Loop through template positions j { if (cell_off[i][j]) S[i][j] = -FLT_MIN; else { // Recursion // NOT the state before the first MM state) CALCULATE_MAX4( S[i][j], B_MM[i][j] - par.mact, // STOP signifies the first MM state, NOT the state before the first MM state (as in Viterbi) S[i-1][j-1] + B_MM[i][j] - par.mact, // B_MM[i][j] contains posterior probability S[i-1][j] - 0.5*par.mact, // gap penalty prevents alignments such as this: XX--xxXX S[i][j-1] - 0.5*par.mact, // YYyy--YY bMM[i][j] // backtracing matrix ); // if (i==6 && j==8) // printf("i=%i j=%i S[i][j]=%8.3f MM:%7.3f MI:%7.3f IM:%7.3f b:%i\n",i,j,S[i][j],S[i-1][j-1]+B_MM[i][j]-par.mact,S[i-1][j],S[i][j-1],bMM[i][j]); // Find maximum score; global alignment: maximize only over last row and last column if(S[i][j]>score_MAC && (par.loc || i==q.L)) { i2=i; j2=j; score_MAC=S[i][j]; } } // end if } //end for j // if global alignment: look for best cell in last column if (!par.loc && S[i][jmax]>score_MAC) { i2=i; j2=jmax; score_MAC=S[i][jmax]; } } // end for i // DEBUG if (v>=5) { printf("\nScore "); for (j=0; j<=t.L; j++) printf("%3i ",j); printf("\n"); for (i=0; i<=q.L; i++) { printf("%2i: ",i); for (j=0; j<=t.L; j++) printf("%5.2f ",S[i][j]); printf("\n"); } printf("\n"); printf("Template=%-12.12s i=%-4i j=%-4i score=%6.3f\n",t.name,i2,j2,score); } return; } /* this is the end of Hit::MACAlignment() */ ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Trace back alignment of two profiles based on matrices bXX[][] */ void Hit::Backtrace(HMM& q, HMM& t) { // Trace back trough the matrices bXY[i][j] until first match state is found (STOP-state) int step; // counts steps in path through 5-layered dynamic programming matrix int i,j; // query and template match state indices InitializeBacktrace(q,t); // Make sure that backtracing stops when t:M1 or q:M1 is reached (Start state), e.g. sMM[i][1], or sIM[i][1] (M:MM, B:IM) for (i=0; i<=q.L; i++) bMM[i][1]=bGD[i][1]=bIM[i][1] = STOP; for (j=1; j<=t.L; j++) bMM[1][j]=bDG[1][j]=bMI[1][j] = STOP; // Back-tracing loop matched_cols=0; // for each MACTH (or STOP) state matched_col is incremented by 1 step=0; // steps through the matrix correspond to alignment columns (from 1 to nsteps) // state=MM; // state with maximum score must be MM state // already set at the end of Viterbi() i=i2; j=j2; // last aligned pair is (i2,j2) while (state) // while (state!=STOP) because STOP=0 { step++; states[step] = state; this->i[step] = i; this->j[step] = j; // Exclude cells in direct neighbourhood from all further alignments for (int ii=imax(i-2,1); ii<=imin(i+2,q.L); ii++) cell_off[ii][j]=1; for (int jj=imax(j-2,1); jj<=imin(j+2,t.L); jj++) cell_off[i][jj]=1; switch (state) { case MM: // current state is MM, previous state is bMM[i][j] matched_cols++; state = bMM[i--][j--]; break; case GD: // current state is GD switch (bGD[i][j--]) { case STOP: state = STOP; break; // current state does not have predecessor case MM: state = MM; break; // previous state is Match state } // default: previous state is same state (GD) break; case IM: switch (bIM[i][j--]) { case STOP: state = STOP; break; // current state does not have predecessor case MM: state = MM; break; // previous state is Match state } // default: previous state is same state (IM) break; case DG: switch (bDG[i--][j]) { case STOP: state = STOP; break; // current state does not have predecessor case MM: state = MM; break; // previous state is Match state } // default: previous state is same state (DG) break; case MI: switch (bMI[i--][j]) { case STOP: state = STOP; break; // current state does not have predecessor case MM: state = MM; break; // previous state is Match state } // default: previous state is same state (MI) break; default: fprintf(stderr,"Error: unallowed state value %i occurred during backtracing at step %i, (i,j)=(%i,%i)\n",state,step,i,j); state=0; v=4; break; } //end switch (state) } //end while (state) i1 = this->i[step]; j1 = this->j[step]; states[step] = MM; // first state (STOP state) is set to MM state nsteps=step; // Allocate new space for alignment scores if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query S = new( float[nsteps+1]); S_ss = new( float[nsteps+1]); if (!S_ss) MemoryError("space for HMM-HMM alignments"); // Add contribution from secondary structure score, record score along alignment, // and record template consensus sequence in master-slave-alignment to query sequence score_ss=0.0f; int ssm=ssm1+ssm2; for (step=1; step<=nsteps; step++) { switch(states[step]) { case MM: i = this->i[step]; j = this->j[step]; S[step] = Score(q.p[i],t.p[j]); S_ss[step] = ScoreSS(q,t,i,j,ssm); score_ss += S_ss[step]; if (Xcons) Xcons[i]=t.Xcons[j]; //record database consensus sequence break; case MI: //if gap in template case DG: if (Xcons) Xcons[this->i[step]]=GAP; //(no break hereafter) default: //if gap in T or Q S[step]=S_ss[step]=0.0f; break; } } if (ssm2>=1) score-=score_ss; // subtract SS score added during alignment!!!! if (Xcons) { for (i=0; i<i1; i++) Xcons[i]=ENDGAP; // set end gap code at beginning and end of template consensus sequence for (i=i2+1; i<=q.L+1; i++) Xcons[i]=ENDGAP; } // Add contribution from correlation of neighboring columns to score float Scorr=0; if (nsteps) { for (step=2; step<=nsteps; step++) Scorr+=S[step]*S[step-1]; for (step=3; step<=nsteps; step++) Scorr+=S[step]*S[step-2]; for (step=4; step<=nsteps; step++) Scorr+=S[step]*S[step-3]; for (step=5; step<=nsteps; step++) Scorr+=S[step]*S[step-4]; score+=par.corr*Scorr; } // Set score, P-value etc. score_sort = score_aass = -score; logPval=0; Pval=1; if (t.mu) { logPvalt=logPvalue(score,t.lamda,t.mu); Pvalt=Pvalue(score,t.lamda,t.mu); } else { logPvalt=0; Pvalt=1;} // printf("%-10.10s lamda=%-9f score=%-9f logPval=%-9g\n",name,t.lamda,score,logPvalt); //DEBUG: Print out Viterbi path if (v>=4) { printf("NAME=%7.7s score=%7.3f score_ss=%7.3f\n",name,score,score_ss); printf("step Q T i j state score T Q cf ss-score\n"); for (step=nsteps; step>=1; step--) { switch(states[step]) { case MM: printf("%4i %1c %1c ",step,q.seq[q.nfirst][this->i[step]],seq[nfirst][this->j[step]]); break; case GD: case IM: printf("%4i - %1c ",step,seq[nfirst][this->j[step]]); break; case DG: case MI: printf("%4i %1c - ",step,q.seq[q.nfirst][this->i[step]]); break; } printf("%4i %4i %2i %7.2f ",this->i[step],this->j[step],(int)states[step],S[step]); printf("%c %c %1i %7.2f\n",i2ss(t.ss_dssp[this->j[step]]),i2ss(q.ss_pred[this->i[step]]),q.ss_conf[this->i[step]]-1,S_ss[step]); } } return; } /* this is the end of Hit::Backtrace() */ ///////////////////////////////////////////////////////////////////////////////////// /** * @brief GLOBAL stochastic trace back through the forward matrix of probability ratios */ void Hit::StochasticBacktrace(HMM& q, HMM& t, char maximize) { int step; // counts steps in path through 5-layered dynamic programming matrix int i,j; // query and template match state indices // float pmin=(par.loc? 1.0: 0.0); // used to distinguish between SW and NW algorithms in maximization const float pmin=0; double* scale_cum = new(double[q.L+2]); scale_cum[0]=1; for (i=1; i<=q.L+1; i++) scale_cum[i] = scale_cum[i-1]*scale[i]; // Select start cell for GLOBAL alignment // (Implementing this in a local version would make this method work for local backtracing as well) if (maximize) { double F_max=0; for (i=q.L-1; i>=1; i--) if (F_MM[i][t.L]/scale_cum[i]>F_max) {i2=i; j2=t.L; F_max=F_MM[i][t.L]/scale_cum[i];} for (j=t.L; j>=1; j--) if (F_MM[q.L][j]/scale_cum[q.L]>F_max) {i2=q.L; j2=j; F_max=F_MM[q.L][j]/scale_cum[q.L];} } else { // float sumF[q.L+t.L]; double* sumF=new(double[q.L+t.L]); sumF[0]=0.0; for (j=1; j<=t.L; j++) sumF[j] = sumF[j-1] + F_MM[q.L][j]/scale_cum[q.L];; for (j=t.L+1; j<t.L+q.L; j++) sumF[j] = sumF[j-1] + F_MM[j-t.L][t.L]/scale_cum[j-t.L];; float x = sumF[t.L+q.L-1]*frand(); // generate random number between 0 and sumF[t.L+q.L-1] for (j=1; j<t.L+q.L; j++) if (x<sumF[j]) break; if (j<=t.L) {i2=q.L; j2=j;} else {i2=j-t.L; j2=t.L;} delete[] sumF; sumF = NULL; } InitializeBacktrace(q,t); int (*pick2)(const double&, const double&, const int&); int (*pick3_GD)(const double&, const double&, const double&); int (*pick3_IM)(const double&, const double&, const double&); int (*pick6)(const double&, const double&, const double&, const double&, const double&, const double&); if (maximize) { pick2 = &pickmax2; pick3_GD = &pickmax3_GD; pick3_IM = &pickmax3_IM; pick6 = &pickmax6; } else { pick2 = &pickprob2; pick3_GD = &pickprob3_GD; pick3_IM = &pickprob3_IM; pick6 = &pickprob6; } // Back-tracing loop matched_cols=0; // for each MACTH (or STOP) state matched_col is incremented by 1 step=0; // steps through the matrix correspond to alignment columns (from 1 to nsteps) state = MM; i=i2; j=j2; // start at end of query and template while (state) // while not reached STOP state or upper or left border { step++; states[step] = state; this->i[step] = i; this->j[step] = j; switch (state) { case MM: // current state is MM, previous state is state // fprintf(stderr,"%4i %1c %1c %4i %4i MM %7.2f\n",step,q.seq[q.nfirst][i],seq[nfirst][j],i,j,Score(q.p[i],t.p[j])); // printf("0:%7.3f MM:%7.3f GD:%7.3f IM:%7.3f DG:%7.3f MI:%7.3f \n", // pmin*scale_cum[i-1], // F_MM[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][M2M], // F_GD[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][D2M], // F_IM[i-1][j-1] * q.tr[i-1][I2M] * t.tr[j-1][M2M], // F_DG[i-1][j-1] * q.tr[i-1][D2M] * t.tr[j-1][M2M], // F_MI[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][I2M]); matched_cols++; if (j>1 && i>1) state = (*pick6)( pmin*scale_cum[i-1], F_MM[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][M2M], F_GD[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][D2M], F_IM[i-1][j-1] * q.tr[i-1][I2M] * t.tr[j-1][M2M], F_DG[i-1][j-1] * q.tr[i-1][D2M] * t.tr[j-1][M2M], F_MI[i-1][j-1] * q.tr[i-1][M2M] * t.tr[j-1][I2M] ); else state=0; i--; j--; break; case GD: // current state is GD // fprintf(stderr,"%4i - %1c %4i %4i GD %7.2f\n",step,q.seq[q.nfirst][j],i,j,Score(q.p[i],t.p[j])); if (j>1) state = (*pick3_GD)( F_MM[i][j-1] * t.tr[j-1][M2D], F_DG[i][j-1] * t.tr[j-1][M2D] * q.tr[i][D2M], // DG -> GD F_GD[i][j-1] * t.tr[j-1][D2D] // gap extension (DD) in template ); else state=0; j--; break; case IM: // fprintf(stderr,"%4i - %1c %4i %4i IM %7.2f\n",step,q.seq[q.nfirst][j],i,j,Score(q.p[i],t.p[j])); if (j>1) state = (*pick3_IM)( F_MM[i][j-1] * q.tr[i][M2I] * t.tr[j-1][M2M_GAPOPEN], F_MI[i][j-1] * q.tr[i][M2I] * t.tr[j-1][I2M], // MI -> IM F_IM[i][j-1] * q.tr[i][I2I] * t.tr[j-1][M2M] // gap extension (II) in query ); else state=0; j--; break; case DG: // fprintf(stderr,"%4i %1c - %4i %4i DG %7.2f\n",step,q.seq[q.nfirst][i],i,j,Score(q.p[i],t.p[j])); if (i>1) state = (*pick2)( F_MM[i-1][j] * q.tr[i-1][M2D] * t.tr[j][GAPOPEN], F_DG[i-1][j] * q.tr[i-1][D2D] * t.tr[j][GAPEXTD], //gap extension (DD) in query DG ); else state=0; i--; break; case MI: // fprintf(stderr,"%4i %1c - %4i %4i MI %7.2f\n",step,q.seq[q.nfirst][i],i,j,Score(q.p[i],t.p[j])); if (i>1) state = (*pick2)( F_MM[i-1][j] * q.tr[i-1][M2M] * t.tr[j][M2I], F_MI[i-1][j] * q.tr[i-1][M2M] * t.tr[j][I2I], //gap extension (II) in template MI ); else state=0; i--; break; } //end switch (state) } //end while (state) i1 = this->i[step]; j1 = this->j[step]; states[step] = MM; // first state (STOP state) is set to MM state nsteps=step; // Allocate new space for alignment scores if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query S = new( float[nsteps+1]); S_ss = new( float[nsteps+1]); if (!S_ss) MemoryError("space for HMM-HMM alignments"); // Add contribution from secondary structure score, record score along alignment, // and record template consensus sequence in master-slave-alignment to query sequence score_ss=0.0f; int ssm=ssm1+ssm2; for (step=1; step<=nsteps; step++) { switch(states[step]) { case MM: i = this->i[step]; j = this->j[step]; S[step] = Score(q.p[i],t.p[j]); S_ss[step] = ScoreSS(q,t,i,j,ssm); score_ss += S_ss[step]; if (Xcons) Xcons[i]=t.Xcons[j]; //record database consensus sequence break; case MI: //if gap in template case DG: if (Xcons) Xcons[this->i[step]]=GAP; //(no break hereafter) default: //if gap in T or Q S[step]=S_ss[step]=0.0f; break; } } if (ssm2>=1) score-=score_ss; // subtract SS score added during alignment!!!! if (Xcons) { for (i=0; i<i1; i++) Xcons[i]=ENDGAP; // set end gap code at beginning and end of template consensus sequence for (i=i2+1; i<=q.L+1; i++) Xcons[i]=ENDGAP; } delete[] scale_cum; scale_cum = NULL; return; } ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Trace back alignment of two profiles based on matrices bXX[][] */ void Hit::BacktraceMAC(HMM& q, HMM& t) { // Trace back trough the matrix b[i][j] until STOP state is found char** b=bMM; // define alias for backtracing matrix int step; // counts steps in path through 5-layered dynamic programming matrix int i,j; // query and template match state indices InitializeBacktrace(q,t); // Make sure that backtracing stops when t:M1 or q:M1 is reached (Start state), e.g. sMM[i][1], or sIM[i][1] (M:MM, B:IM) for (i=0; i<=q.L; i++) b[i][1] = STOP; for (j=1; j<=t.L; j++) b[1][j] = STOP; // Back-tracing loop // In contrast to the Viterbi-Backtracing, STOP signifies the first Match-Match state, NOT the state before the first MM state matched_cols=1; // for each MACTH (or STOP) state matched_col is incremented by 1 state=MM; // lowest state with maximum score must be match-match state step=0; // steps through the matrix correspond to alignment columns (from 1 to nsteps) i=i2; j=j2; // last aligned pair is (i2,j2) while (state!=STOP) { step++; states[step] = state = b[i][j]; this->i[step] = i; this->j[step] = j; // Exclude cells in direct neighbourhood from all further alignments for (int ii=imax(i-2,1); ii<=imin(i+2,q.L); ii++) cell_off[ii][j]=1; for (int jj=imax(j-2,1); jj<=imin(j+2,t.L); jj++) cell_off[i][jj]=1; if (state==MM) matched_cols++; switch (state) { case MM: i--; j--; break; case IM: j--; break; case MI: i--; break; case STOP: break; default: fprintf(stderr,"Error: unallowed state value %i occurred during backtracing at step %i, (i,j)=(%i,%i)\n",state,step,i,j); state=0; v=4; break; } //end switch (state) } //end while (state) i1 = this->i[step]; j1 = this->j[step]; states[step] = MM; // first state (STOP state) is set to MM state nsteps=step; // Allocate new space for alignment scores if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query S = new( float[nsteps+1]); S_ss = new( float[nsteps+1]); P_posterior = new( float[nsteps+1]); if (!P_posterior) MemoryError("space for HMM-HMM alignments"); // Add contribution from secondary structure score, record score along alignment, // and record template consensus sequence in master-slave-alignment to query sequence score_ss=0.0f; sum_of_probs=0.0; // number of identical residues in query and template sequence int ssm=ssm1+ssm2; // printf("Hit=%s\n",name); ///////////////////////////////////////////////////////////// for (step=1; step<=nsteps; step++) { switch(states[step]) { case MM: i = this->i[step]; j = this->j[step]; S[step] = Score(q.p[i],t.p[j]); S_ss[step] = ScoreSS(q,t,i,j,ssm); score_ss += S_ss[step]; P_posterior[step] = B_MM[this->i[step]][this->j[step]]; // Add probability to sum of probs if no dssp states given or dssp states exist and state is resolved in 3D structure if (t.nss_dssp<0 || t.ss_dssp[j]>0) sum_of_probs += P_posterior[step]; // printf("j=%-3i dssp=%1i P=%4.2f sum=%6.2f\n",j,t.ss_dssp[j],P_posterior[step],sum_of_probs); ////////////////////////// if (Xcons) Xcons[i]=t.Xcons[j]; //record database consensus sequence break; case MI: //if gap in template case DG: if (Xcons) Xcons[this->i[step]]=GAP; //(no break hereafter) default: //if gap in T or Q S[step] = S_ss[step] = P_posterior[step] = 0.0; break; } } // printf("\n"); ///////////////////////////////////////////////////////////// if (ssm2>=1) score-=score_ss; // subtract SS score added during alignment!!!! if (Xcons) { for (i=0; i<i1; i++) Xcons[i]=ENDGAP; // set end gap code at beginning and end of template consensus sequence for (i=i2+1; i<=q.L+1; i++) Xcons[i]=ENDGAP; } // Add contribution from correlation of neighboring columns to score float Scorr=0; if (nsteps) { for (step=1; step<=nsteps-1; step++) Scorr+=S[step]*S[step+1]; for (step=1; step<=nsteps-2; step++) Scorr+=S[step]*S[step+2]; for (step=1; step<=nsteps-3; step++) Scorr+=S[step]*S[step+3]; for (step=1; step<=nsteps-4; step++) Scorr+=S[step]*S[step+4]; score+=par.corr*Scorr; } // Set score, P-value etc. score_sort = score_aass = -score; logPval=0; Pval=1; if (t.mu) { logPvalt=logPvalue(score,t.lamda,t.mu); Pvalt=Pvalue(score,t.lamda,t.mu); } else { logPvalt=0; Pvalt=1;} // printf("%-10.10s lamda=%-9f score=%-9f logPval=%-9g\n",name,t.lamda,score,logPvalt); //DEBUG: Print out MAC alignment path if (v>=4) { float sum_post=0.0; printf("NAME=%7.7s score=%7.3f score_ss=%7.3f\n",name,score,score_ss); printf("step Q T i j state score T Q cf ss-score P_post Sum_post\n"); for (step=nsteps; step>=1; step--) { switch(states[step]) { case MM: sum_post+=P_posterior[step]; printf("%4i %1c %1c ",step,q.seq[q.nfirst][this->i[step]],seq[nfirst][this->j[step]]); break; case IM: printf("%4i - %1c ",step,seq[nfirst][this->j[step]]); break; case MI: printf("%4i %1c - ",step,q.seq[q.nfirst][this->i[step]]); break; } printf("%4i %4i %2i %7.1f ",this->i[step],this->j[step],(int)states[step],S[step]); printf("%c %c %1i %7.1f ",i2ss(t.ss_dssp[this->j[step]]),i2ss(q.ss_pred[this->i[step]]),q.ss_conf[this->i[step]]-1,S_ss[step]); printf("%7.5f %7.2f\n",P_posterior[step],sum_post); } } return; } ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Functions that calculate probabilities */ void Hit::InitializeForAlignment(HMM& q, HMM& t) { int i,j; // SS scoring during (ssm2>0) or after (ssm1>0) alignment? Query SS known or Template SS known? switch (par.ssm) { case 0: ssm1=0; ssm2=0; break; case 1: ssm2=0; // SS scoring after alignment if (t.nss_dssp>=0 && q.nss_pred>=0) ssm1=1; else if (q.nss_dssp>=0 && t.nss_pred>=0) ssm1=2; else if (q.nss_pred>=0 && t.nss_pred>=0) ssm1=3; else ssm1=0; break; case 2: ssm1=0; // SS scoring during alignment if (t.nss_dssp>=0 && q.nss_pred>=0) ssm2=1; else if (q.nss_dssp>=0 && t.nss_pred>=0) ssm2=2; else if (q.nss_pred>=0 && t.nss_pred>=0) ssm2=3; else ssm2=0; break; case 3: ssm2=0; // SS scoring after alignment if (q.nss_pred>=0 && t.nss_pred>=0) ssm1=3; else ssm1=0; break; case 4: ssm1=0; // SS scoring during alignment if (q.nss_pred>=0 && t.nss_pred>=0) ssm2=3; else ssm2=0; break; // case 5: // ssm2=0; // SS scoring after alignment // if (q.nss_dssp>=0 && t.nss_dssp>=0) ssm1=4; else ssm1=0; // break; // case 6: // ssm1=0; // SS scoring during alignment // if (q.nss_dssp>=0 && t.nss_dssp>=0) ssm2=4; else ssm2=0; // break; } if (self) { // Cross out cells in lower diagonal for self-comparison? for (i=1; i<=q.L; i++) { int jmax = imin(i+SELFEXCL,t.L); for (j=1; j<=jmax; j++) cell_off[i][j]=1; // cross out cell near diagonal for (j=jmax+1; j<=t.L+1; j++) cell_off[i][j]=0; // no other cells crossed out yet } } else // Compare two different HMMs Q and T { // Activate all cells in dynamic programming matrix for (i=1; i<=q.L; i++) for (j=1; j<=t.L; j++) cell_off[i][j]=0; // no other cells crossed out yet // Cross out cells that are excluded by the minimum-overlap criterion if (par.min_overlap==0) min_overlap = imin(60, (int)(0.333f*imin(q.L,t.L))+1); // automatic minimum overlap else min_overlap = imin(par.min_overlap, (int)(0.8f*imin(q.L,t.L))); for (i=0; i<min_overlap; i++) for (j=i-min_overlap+t.L+1; j<=t.L; j++) // Lt-j+i>=Ovlap => j<=i-Ovlap+Lt => jmax=min{Lt,i-Ovlap+Lt} cell_off[i][j]=1; for (i=q.L-min_overlap+1; i<=q.L; i++) for (j=1; j<i+min_overlap-q.L; j++) // Lq-i+j>=Ovlap => j>=i+Ovlap-Lq => jmin=max{1, i+Ovlap-Lq} cell_off[i][j]=1; } // Cross out rows which are contained in range given by exclstr ("3-57,238-314") if (par.exclstr) { char* ptr=par.exclstr; int i0, i1; while (1) { i0 = abs(strint(ptr)); i1 = abs(strint(ptr)); if (!ptr) break; for (i=i0; i<=imin(i1,q.L); i++) for (j=1; j<=t.L; j++) cell_off[i][j]=1; } } } ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Allocate memory for data of new alignment (sequence names, alignment, scores,...) */ void Hit::InitializeBacktrace(HMM& q, HMM& t) { if (irep==1) //if this is the first single repeat repeat hit with this template { //Copy information about template profile to hit and reset template pointers to avoid destruction longname=new(char[strlen(t.longname)+1]); name =new(char[strlen(t.name)+1]); file =new(char[strlen(t.file)+1]); if (!file) MemoryError("space for alignments with database HMMs. \nNote that all alignments have to be kept in memory"); strcpy(longname,t.longname); strcpy(name,t.name); strcpy(fam ,t.fam); strcpy(sfam ,t.sfam); strcpy(fold ,t.fold); strcpy(cl ,t.cl); strcpy(file,t.file); sname=new(char*[t.n_display]); // Call Compare only once with irep=1 seq =new(char*[t.n_display]); // Call Compare only once with irep=1 if (!sname || !seq) MemoryError("space for alignments with database HMMs.\nNote that all sequences for display have to be kept in memory"); for (int k=0; k<t.n_display; k++) { if (NULL != t.sname){ sname[k]=t.sname[k]; t.sname[k]=NULL; } else { sname[k]=NULL; } seq[k] =t.seq[k]; t.seq[k]=NULL; } n_display=t.n_display; t.n_display=0; ncons = t.ncons; nfirst = t.nfirst; nss_dssp = t.nss_dssp; nsa_dssp = t.nsa_dssp; nss_pred = t.nss_pred; nss_conf = t.nss_conf; L = t.L; Neff_HMM = t.Neff_HMM; Eval = 1.0; Pval = 1.0; Pvalt = 1.0; logPval = 0.0; logPvalt= 0.0; Probab = 1.0; } // Allocate new space this->i = new( int[i2+j2+2]); this->j = new( int[i2+j2+2]); states = new( char[i2+j2+2]); S = S_ss = P_posterior = NULL; // set to NULL to avoid deleting data from irep=1 when hit with irep=2 is removed Xcons = NULL; } ///////////////////////////////////////////////////////////////////////////////////// // Some score functions ///////////////////////////////////////////////////////////////////////////////////// /** * @brief Calculate score between columns i and j of two HMMs (query and template) */ inline float Score(float* qi, float* tj) { // if (par.columnscore==9) // return (tj[0] *qi[0] +tj[1] *qi[1] +tj[2] *qi[2] +tj[3] *qi[3] +tj[4]*qi[4] // +tj[5] *qi[5] +tj[6] *qi[6] +tj[7] *qi[7] +tj[8] *qi[8] +tj[9]*qi[9] // +tj[10]*qi[10]+tj[11]*qi[11]+tj[12]*qi[12]+tj[13]*qi[13]+tj[14]*qi[14] // +tj[15]*qi[15]+tj[16]*qi[16]+tj[17]*qi[17]+tj[18]*qi[18]+tj[19]*qi[19]); // else return fast_log2( tj[0] *qi[0] +tj[1] *qi[1] +tj[2] *qi[2] +tj[3] *qi[3] +tj[4] *qi[4] +tj[5] *qi[5] +tj[6] *qi[6] +tj[7] *qi[7] +tj[8] *qi[8] +tj[9] *qi[9] +tj[10]*qi[10]+tj[11]*qi[11]+tj[12]*qi[12]+tj[13]*qi[13]+tj[14]*qi[14] +tj[15]*qi[15]+tj[16]*qi[16]+tj[17]*qi[17]+tj[18]*qi[18]+tj[19]*qi[19] ); } /** * @brief Calculate score between columns i and j of two HMMs (query and template) */ inline float ProbFwd(float* qi, float* tj) { return tj[0] *qi[0] +tj[1] *qi[1] +tj[2] *qi[2] +tj[3] *qi[3] +tj[4] *qi[4] +tj[5] *qi[5] +tj[6] *qi[6] +tj[7] *qi[7] +tj[8] *qi[8] +tj[9] *qi[9] +tj[10]*qi[10]+tj[11]*qi[11]+tj[12]*qi[12]+tj[13]*qi[13]+tj[14]*qi[14] +tj[15]*qi[15]+tj[16]*qi[16]+tj[17]*qi[17]+tj[18]*qi[18]+tj[19]*qi[19]; } /** * @brief Calculate secondary structure score between columns i and j of two HMMs (query and template) */ inline float Hit::ScoreSS(HMM& q, HMM& t, int i, int j, int ssm) { switch (ssm) //SS scoring during alignment { case 0: // no SS scoring during alignment return 0.0; case 1: // t has dssp information, q has psipred information return par.ssw * S73[ (int)t.ss_dssp[j]][ (int)q.ss_pred[i]][ (int)q.ss_conf[i]]; case 2: // q has dssp information, t has psipred information return par.ssw * S73[ (int)q.ss_dssp[i]][ (int)t.ss_pred[j]][ (int)t.ss_conf[j]]; case 3: // q has dssp information, t has psipred information return par.ssw * S33[ (int)q.ss_pred[i]][ (int)q.ss_conf[i]][ (int)t.ss_pred[j]][ (int)t.ss_conf[j]]; // case 4: // q has dssp information, t has dssp information // return par.ssw*S77[ (int)t.ss_dssp[j]][ (int)t.ss_conf[j]]; } return 0.0; } /** * @brief Calculate secondary structure score between columns i and j of two HMMs (query and template) */ inline float Hit::ScoreSS(HMM& q, HMM& t, int i, int j) { return ScoreSS(q,t,i,j,ssm2); } /** * @brief Calculate score between columns i and j of two HMMs (query and template) */ inline float Hit::ScoreTot(HMM& q, HMM& t, int i, int j) { return Score(q.p[i],t.p[j]) + ScoreSS(q,t,i,j) + par.shift; } /* * Calculate score between columns i and j of two HMMs (query and template) */ inline float Hit::ScoreAA(HMM& q, HMM& t, int i, int j) { return Score(q.p[i],t.p[j]); } ///////////////////////////////////////////////////////////////////////////////////// /* * Function for Viterbi() */ inline float max2(const float& xMM, const float& xX, char& b) { if (xMM>xX) { b=MM; return xMM;} else { b=SAME; return xX;} } ///////////////////////////////////////////////////////////////////////////////////// /* * Functions for StochasticBacktrace() */ inline int pickprob2(const double& xMM, const double& xX, const int& state) { if ( (xMM+xX)*frand() < xMM) return MM; else return state; } inline int pickprob3_GD(const double& xMM, const double& xDG, const double& xGD) { double x = (xMM+xDG+xGD)*frand(); if ( x<xMM) return MM; else if ( x<xMM+xDG) return DG; else return GD; } inline int pickprob3_IM(const double& xMM, const double& xMI, const double& xIM) { double x = (xMM+xMI+xIM)*frand(); if ( x<xMM) return MM; else if ( x<xMM+xMI) return MI; else return IM; } inline int pickprob6(const double& x0, const double& xMM, const double& xGD, const double& xIM, const double& xDG, const double& xMI) { double x = (x0+xMM+xGD+xIM+xDG+xMI)*frand(); x-=xMM; if (x<0) return MM; x-=x0; if (x<0) return STOP; x-=xGD; if (x<0) return GD; x-=xIM; if (x<0) return IM; if (x < xDG) return DG; else return MI; } inline int pickmax2(const double& xMM, const double& xX, const int& state) { if (xMM > xX) return MM; else return state; } inline int pickmax3_GD(const double& xMM, const double& xDG, const double& xGD) { char state; double x; if ( xMM>xDG) {state=MM; x=xMM;} else {state=DG; x=xDG;} if ( xGD>x) {state=GD; x=xGD;} return state; } inline int pickmax3_IM(const double& xMM, const double& xMI, const double& xIM) { char state; double x; if ( xMM>xMI) {state=MM; x=xMM;} else {state=MI; x=xMI;} if ( xIM>x) {state=IM; x=xIM;} return state; } inline int pickmax6(const double& x0, const double& xMM, const double& xGD, const double& xIM, const double& xDG, const double& xMI) { char state; double x; if ( x0 >xMM) {state=STOP; x=x0;} else {state=MM; x=xMM;} if ( xGD>x) {state=GD; x=xGD;} if ( xIM>x) {state=IM; x=xIM;} if ( xDG>x) {state=DG; x=xDG;} if ( xMI>x) {state=MI; x=xMI;} return state; } ///////////////////////////////////////////////////////////////////////////////////// //// Functions that calculate P-values and probabilities ///////////////////////////////////////////////////////////////////////////////////// //// Evaluate the CUMULATIVE extreme value distribution at point x //// p(s)ds = lamda * exp{ -exp[-lamda*(s-mu)] - lamda*(s-mu) } ds = exp( -exp(-x) - x) dx = p(x) dx //// => P(s>S) = integral_-inf^inf {p(x) dx} = 1 - exp{ -exp[-lamda*(S-mu)] } inline double Pvalue(double x, double a[]) { //a[0]=lamda, a[1]=mu double h = a[0]*(x-a[1]); return (h>10)? exp(-h) : double(1.0)-exp( -exp(-h)); } inline double Pvalue(float x, float lamda, float mu) { double h = lamda*(x-mu); return (h>10)? exp(-h) : (double(1.0)-exp( -exp(-h))); } inline double logPvalue(float x, float lamda, float mu) { double h = lamda*(x-mu); return (h>10)? -h : (h<-2.5)? -exp(-exp(-h)): log( ( double(1.0) - exp(-exp(-h)) ) ); } inline double logPvalue(float x, double a[]) { double h = a[0]*(x-a[1]); return (h>10)? -h : (h<-2.5)? -exp(-exp(-h)): log( ( double(1.0) - exp(-exp(-h)) ) ); } // Calculate probability of true positive : p_TP(score)/( p_TP(score)+p_FP(score) ) // TP: same superfamily OR MAXSUB score >=0.1 inline double Probab(Hit& hit) { double s=-hit.score_aass; double t; if (s>200) return 100.0; if (par.loc) { if (par.ssm && (hit.ssm1 || hit.ssm2) && par.ssw>0) { // local with SS const double a=sqrt(6000.0); const double b=2.0*2.5; const double c=sqrt(0.12); const double d=2.0*32.0; t = a*exp(-s/b) + c*exp(-s/d); } else { // local no SS const double a=sqrt(4000.0); const double b=2.0*2.5; const double c=sqrt(0.15); const double d=2.0*34.0; t = a*exp(-s/b) + c*exp(-s/d); } } else { if ( (par.ssm>0) && (par.ssw>0) ) /* FIXME: was '&', should be '&&' (or not?) */ { // global with SS const double a=sqrt(4000.0); const double b=2.0*3.0; const double c=sqrt(0.13); const double d=2.0*34.0; t = a*exp(-s/b) + c*exp(-s/d); } else { // global no SS const double a=sqrt(6000.0); const double b=2.0*2.5; const double c=sqrt(0.10); const double d=2.0*37.0; t = a*exp(-s/b) + c*exp(-s/d); } } return 100.0/(1.0+t*t); } // #define Weff(Neff) (1.0+par.neffa*(Neff-1.0)+(par.neffb-4.0*par.neffa)/16.0*(Neff-1.0)*(Neff-1.0)) // ///////////////////////////////////////////////////////////////////////////////////// // // Merge HMM with next aligned HMM // ///////////////////////////////////////////////////////////////////////////////////// // void Hit::MergeHMM(HMM& Q, HMM& t, float wk[]) // { // int i,j; // position in query and target // int a; // amino acid // int step; // alignment position (step=1 is end) // float Weff_M, Weff_D, Weff_I; // for (step=nsteps; step>=2; step--) // iterate only to one before last alignment column // { // i = this->i[step]; // j = this->j[step]; // switch(states[step]) // { // case MM: // Weff_M = Weff(t.Neff_M[j]-1.0); // Weff_D = Weff(t.Neff_D[j]-1.0); // Weff_I = Weff(t.Neff_I[j]-1.0); // for (a=0; a<20; a++) Q.f[i][a] += t.f[j][a]*wk[j]*Weff_M; // switch(states[step-1]) // { // case MM: // MM->MM // Q.tr_lin[i][M2M]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][M2D]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][M2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2M]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case MI: // MM->MI // Q.tr_lin[i][M2D]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][M2D]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][M2M]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case DG: // MM->DG // Q.tr_lin[i][M2D]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][M2D]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][M2M]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case IM: // MM->IM // Q.tr_lin[i][M2I]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][M2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][M2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2M]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case GD: // MM->GD // Q.tr_lin[i][M2I]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][M2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][M2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2M]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // } // break; // case MI: // if gap in template // Weff_I = Weff(t.Neff_I[j]-1.0); // switch(states[step-1]) // { // case MI: // MI->MI // Q.tr_lin[i][M2M]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case MM: // MI->MM // Q.tr_lin[i][M2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // break; // } // break; // case DG: // Weff_M = Weff(t.Neff_M[j]-1.0); // Weff_D = Weff(t.Neff_D[j]-1.0); // Weff_I = Weff(t.Neff_I[j]-1.0); // switch(states[step-1]) // { // case DG: // DG->DG // Q.tr_lin[i][D2D]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][M2M]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // Q.tr_lin[i][M2D]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // break; // case MM: // DG->MM // Q.tr_lin[i][D2M]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][D2D]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][D2M]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // Q.tr_lin[i][M2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // break; // } // break; // case IM: // if gap in query // Weff_M = Weff(t.Neff_M[j]-1.0); // switch(states[step-1]) // { // case IM: // IM->IM // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // break; // case MM: // IM->MM // Weff_D = Weff(t.Neff_D[j]-1.0); // Weff_I = Weff(t.Neff_I[j]-1.0); // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][I2M]+= t.tr_lin[j][D2M]*wk[j]*Weff_D; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // } // break; // case GD: // Weff_M = Weff(t.Neff_M[j]-1.0); // switch(states[step-1]) // { // case GD: // GD->GD // Weff_I = Weff(t.Neff_I[j]-1.0); // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // case MM: // GD->MM // Weff_D = Weff(t.Neff_D[j]-1.0); // Weff_I = Weff(t.Neff_I[j]-1.0); // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2M]*wk[j]*Weff_M; // Q.tr_lin[i][I2M]+= t.tr_lin[j][M2D]*wk[j]*Weff_M; // Q.tr_lin[i][I2I]+= t.tr_lin[j][M2I]*wk[j]*Weff_M; // Q.tr_lin[i][D2D]+= t.tr_lin[j][D2D]*wk[j]*Weff_D; // Q.tr_lin[i][I2M]+= t.tr_lin[j][I2M]*wk[j]*Weff_I; // Q.tr_lin[i][I2I]+= t.tr_lin[j][I2I]*wk[j]*Weff_I; // break; // } // break; // } // } // i = this->i[step]; // j = this->j[step]; // Weff_M = Weff(t.Neff_M[j]-1.0); // for (a=0; a<20; a++) Q.f[i][a] += t.f[j][a]*wk[j]*Weff_M; // } #ifdef CLUSTALO /* @* Hit::ClobberGlobal (eg, hit) * */ void Hit::ClobberGlobal(void){ if (i){ //delete[] i; i = NULL; } if (j){ //delete[] j; j = NULL; } if (states){ //delete[] states; states = NULL; } if (S){ //delete[] S; S = NULL; } if (S_ss){ //delete[] S_ss; S_ss = NULL; } if (P_posterior){ //delete[] P_posterior; P_posterior = NULL; } if (Xcons){ //delete[] Xcons; Xcons = NULL; } // delete[] l; l = NULL; i = j = NULL; states = NULL; S = S_ss = P_posterior = NULL; Xcons = NULL; if (irep==1) // if irep>1 then longname etc point to the same memory locations as the first repeat. { // but these have already been deleted. // printf("Delete name = %s\n",name);////////////////////////// //delete[] longname; longname = NULL; //delete[] name; name = NULL; //delete[] file; file = NULL; //delete[] dbfile; dbfile = NULL; /*for (int k=0; k<n_display; k++) { delete[] sname[k]; sname[k] = NULL; delete[] seq[k]; seq[k] = NULL; }*/ //delete[] sname; sname = NULL; //delete[] seq; seq = NULL; } score = score_sort = score_aass = 0.0; Pval = Pvalt = Eval = Probab = 0; Pforward = sum_of_probs = 0.00; L = irep = nrep = n_display = nsteps = 0; i1 = i2 = j1 = j2 = matched_cols = min_overlap = 0; } #endif /* * EOF hhhit-C.h */