********************************************************************/
/*
- * 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:
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
* @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
*
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[]) {
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
// 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();
tL = strlen(ppcSecndProf[0]);
}
else {
- tL = prHMM->L;
+ tL = prHMM2->L;
}
// 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."<<endl;
par.forward=0;
}
-
// Read input file (HMM, HHM, or alignment format), and add pseudocounts etc.
q.cQT = 'q';
if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
and add pseudocounts etc. */
t.cQT = 't';
if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
- ppcSecndProf, iSecndCnt, prHMM,
+ ppcSecndProf, iSecndCnt, prHMM2,
pcPrealigned2, pcRepresent2, pdWeightsR,
infile, t)) {
sprintf(zcError, "%s:%s:%d: Problem Reading/Preparing t-profile (len=%d)\n",
// Factor Null model into HMM t
t.IncludeNullModelInHMM(q,t);
+ /* alignment will fail if one profile does not contain useful characters, FS, r259 -> 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;
/* 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");
}
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;
}
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;
hitlist.PrintHitList(q, par.outfile);
hitlist.PrintAlignments(
#ifdef CLUSTALO
- ppcFirstProf, ppcSecndProf,
+ ppcFirstProf, ppcSecndProf, zcAux, zcError,
#endif
q, par.outfile);
if (v>=2) {
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
hitlist.Reset();
while (!hitlist.End())
hitlist.ReadNext().Delete(); // Delete content of hit object
-
+
if (strucfile && par.wstruc>0)
{
for (int i=0; i<q.L+2; i++){
t.ClobberGlobal();
}
hitlist.ClobberGlobal();
-
+
return iRetVal;
} /* this is the end of hhalign() //end main */
git://source.jalview.org / jabaws.git / blobdiff