X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=binaries%2Fsrc%2Fclustalo%2Fsrc%2Fhhalign%2Fhhalign.cpp;h=6fdcc964dac5f10f449a35852597dae155fc69fd;hb=1eab6cadc31107fabe7ef060de1ae91eeb4a025b;hp=c2f73292fa14389a4c8e8bcb9aab2999015b4e28;hpb=b222b987978b6b2d4b6fb1b4a3a2a0ea87a19c3d;p=jabaws.git diff --git a/binaries/src/clustalo/src/hhalign/hhalign.cpp b/binaries/src/clustalo/src/hhalign/hhalign.cpp index c2f7329..6fdcc96 100644 --- a/binaries/src/clustalo/src/hhalign/hhalign.cpp +++ b/binaries/src/clustalo/src/hhalign/hhalign.cpp @@ -15,7 +15,7 @@ ********************************************************************/ /* - * RCS $Id: hhalign.cpp 250 2011-06-17 13:06:20Z fabian $ + * RCS $Id: hhalign.cpp 309 2016-06-13 14:10:11Z fabian $ */ /* hhalign.C: @@ -113,10 +113,7 @@ float** Pstruc=NULL; /* structure matrix which can be multiplied to prob ratios from aa column comparisons in Forward() */ float** Sstruc=NULL; /* structure matrix which can be added to log odds from aa column comparisons in Forward() */ - - - - +bool nucleomode = false; //////////////////////////////////////////////////////////////////////////// // Help functions @@ -670,6 +667,8 @@ ProcessArguments(int argc, char** argv) * @param[in] rHhalignPara, * various parameters passed down from commandline * e.g., iMaxRamMB + * @param[out] rHHscores + * qualify goodness of alignment * * iFlag,zcAux,zcError are debugging arguments * @@ -678,11 +677,11 @@ ProcessArguments(int argc, char** argv) extern "C" int hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, char **ppcSecndProf, int iSecndCnt, double *pdWeightsR, - double *dScore_p, hmm_light *prHMM, + double *dScore_p, hmm_light *prHMM, hmm_light *prHMM2, char *pcPrealigned1, char *pcRepresent1, char *pcPrealigned2, char *pcRepresent2, - hhalign_para rHhalignPara, - int iFlag, + hhalign_para rHhalignPara, hhalign_scores *rHHscores, + int iFlag, int iVerbosity, char zcAux[], char zcError[]) { @@ -715,7 +714,7 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, int Nali; /* number of normally backtraced alignments in dot plot */ - SetDefaults(); + SetDefaults(rHhalignPara); strcpy(par.tfile,""); /* FS, Initialise Argument */ strcpy(par.alnfile,""); par.p=0.0 ; /* minimum threshold for inclusion in hit list @@ -845,8 +844,21 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, // Set (global variable) substitution matrix and derived matrices - SetSubstitutionMatrix(); - + // DD: new experimental matrices and params for nucleotides + if(rHhalignPara.bIsDna) + { + nucleomode = true; + SetDnaDefaults(rHhalignPara); + SetDnaSubstitutionMatrix(); + } + else if(rHhalignPara.bIsRna) + { + nucleomode = true; + SetRnaDefaults(rHhalignPara); + SetRnaSubstitutionMatrix(); + } + else + SetSubstitutionMatrix(); // Set secondary structure substitution matrix SetSecStrucSubstitutionMatrix(); @@ -867,7 +879,7 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, tL = strlen(ppcSecndProf[0]); } else { - tL = prHMM->L; + tL = prHMM2->L; } @@ -877,6 +889,8 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, // longest allowable length of database HMM int Lmaxmem=(int)((float)MEMSPACE_DYNPROG/qL/6/8); if (par.forward==2 && tL+2>=Lmaxmem) { + sprintf(zcError, "%s:%s:%d: switch to Viterbi (qL=%d, tL=%d, MAC-RAM=%d\n)", + __FUNCTION__, __FILE__, __LINE__, qL, tL, rHhalignPara.iMacRamMB); if (v>=1) cerr<<"WARNING: Not sufficient memory to realign with MAC algorithm. Using Viterbi algorithm."< r260 */ + if ( (q.L <= 0) || (t.L <= 0) ){ + sprintf(zcError, "%s:%s:%d: Problem Reading/Preparing profiles (len(q)=%d/len(t)=%d)\n", + __FUNCTION__, __FILE__, __LINE__, q.L, t.L); + iRetVal = RETURN_FROM_RNP; + goto this_is_the_end; + } + /* switch at this stage is ineffectual as log/lin already decided in RnP(). FS, r228 -> r229 */ /*const float MEMSPACE_DYNPROG = 512*1024*1024; @@ -1007,10 +1028,11 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, /* Do (self-)comparison, store results if score>SMIN, and try next best alignment */ /* FIXME very ambigous and possibly faulty if-else */ - if (v>=2) { - if (par.forward==2) { - printf("Using maximum accuracy (MAC) alignment algorithm ...\n"); - } + if (v>=2) + { + if (par.forward==2) { + printf("Using maximum accuracy (MAC) alignment algorithm ...\n"); + } else if (par.forward==0) { printf("Using Viterbi algorithm ...\n"); } @@ -1064,7 +1086,7 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, hit.BacktraceMAC(q,t); if ((isnan(hit.score)) || (isnan(hit.Pforward))){ printf("nan after backtrace\n"); - } + } } /* use MAC algorithm */ // printf ("%-12.12s %-12.12s irep=%-2i score=%6.2f hit.Pvalt=%.2g\n",hit.name,hit.fam,hit.irep,hit.score,hit.Pvalt); *dScore_p = hit.score; @@ -1188,12 +1210,13 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, } int iPrAliRtn = hitlist.PrintAlignments( #ifdef CLUSTALO - ppcFirstProf, ppcSecndProf, + ppcFirstProf, ppcSecndProf, zcAux, zcError, #endif q, par.pairwisealisfile, par.outformat); - if (OK != iPrAliRtn){ - fprintf(stderr, "%s:%s:%d: Could not print alignments\n", + if ( (OK != iPrAliRtn) ){ + sprintf(zcAux, "%s:%s:%d: Could not print alignments\n", __FUNCTION__, __FILE__, __LINE__); + strcat(zcError, zcAux); iRetVal = RETURN_FROM_PRINT_ALI; /* this is where mis-alignment was originally spotted, hope to trap it now earlier. FS, r241 -> r243 */ goto this_is_the_end; @@ -1207,7 +1230,7 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, hitlist.PrintHitList(q, par.outfile); hitlist.PrintAlignments( #ifdef CLUSTALO - ppcFirstProf, ppcSecndProf, + ppcFirstProf, ppcSecndProf, zcAux, zcError, #endif q, par.outfile); if (v>=2) { @@ -1291,7 +1314,20 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, q.name,t.name,par.hitrank,hit.score,hit.Pval); } - + rHHscores->hhScore = hit.score; + rHHscores->forwardProb = hit.Pforward; + rHHscores->sumPP = hit.sum_of_probs; + + /* next few lines commented out as they caused segfaults in RNA mode */ + /* rHHscores->PP = (double *)malloc((hit.L+GOOD_MEASURE) * sizeof(double));*/ + rHHscores->L = hit.L; +/* for (int i = 0; i < hit.L; i++){ + cout << "rHHscores->PP[i]" << rHHscores->PP[i] << endl; + cout << "hit.P_posterior[i+1]" << hit.P_posterior[i+1] << endl; + rHHscores->PP[i] = hit.P_posterior[i+1]; + } */ + + this_is_the_end: // Delete memory for dynamic programming matrix @@ -1308,7 +1344,7 @@ hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL, hitlist.Reset(); while (!hitlist.End()) hitlist.ReadNext().Delete(); // Delete content of hit object - + if (strucfile && par.wstruc>0) { for (int i=0; i