+++ /dev/null
-/* -*- 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: hhfunc-C.h 256 2011-06-23 13:57:28Z fabian $
- */
-
-/*
- * Changelog: Michael Remmert made changes to hhalign stand-alone code
- * FS implemented some of the changes on 2010-11-11 -> MR1
- *
- * did not incorporate SS prediction PSIpred (yet); functions are:
- * CalculateSS(3),
- */
-
-
-// hhfunc.C
-
-
-/**
- * AddBackgroundFrequencies()
- *
- * @brief add background frequencies (derived from HMM) to
- * sequence/profile
- *
- * @param[in,out] ppfFreq,
- * [in] residue frequencies of sequence/profile,
- * [out] overlayed with HMM background frequencies
- * @param[in,out] ppfPseudoF,
- * [in] residue frequencies+pseudocounts of sequence/profile,
- * [out] overlayed with HMM background frequencies+pseudocounts
- * @param[in] iSeqLen,
- * length of sequence/profile (not aligned to HMM)
- * @pram[in] prHMM,
- * background HMM
- * @param[in] ppcSeq,
- * sequences/profile to be 'softened'
- * @param[in] pcPrealigned,
- * sequence aligned to HMM, this is not quite a consensus,
- * it is identical to 1st sequence but over-writes gap information,
- * if other sequences in profile have (non-gap) residues
- * @param[in] iPreCnt
- * number of sequences pre-aligned (pcPrealigned is 'consensus' of these sequences)
- * @param[in] pcRepresent
- * sequence representative of HMM, aligned to pcSeq0
- */
-void
-AddBackgroundFrequencies(float **ppfFreq, float **ppfPseudoF, float **ppfPseudoTr,
- int iSeqLen, hmm_light *prHMM,
- char **ppcSeq, char *pcPrealigned, int iPreCnt, char *pcRepresent)
-{
-
- char *pcS = pcPrealigned; /* current residue in pre-aligned sequence */
- char *pcH = pcRepresent; /* current residue in pre-aligned HMM */
- int iS = 0; /* position in un-aligned sequence (corresponds to pcS) */
- int iH = 0; /* position in un-aligned HMM (corresponds to pcH) */
- int iA; /* residue iterator */
- int iT; /* transition state iterator */
- float fFWeight = 0.50 / sqrt((float)(iPreCnt)); /* weight of backgroud frequencies */ /* FIXME: tune value, 0.50 default */
- //float fFWeight = 0.75;
- float fFThgiew = UNITY - fFWeight; /* weight of 'true' frequencies */
- float fGWeight = 0.50 / sqrt((float)(iPreCnt)); /* weight of backgroud transitions */ /* FIXME: tune value, 0.50 default */
- //float fGWeight = 0.50 /*/ (float)(iPreCnt)*/; /* weight of backgroud transitions */ /* FIXME: tune value, 0.50 default */
- float fGThgiew = UNITY - fGWeight; /* weight of 'true' transitions */
- float fAux;
- /* zf1SeqTrans[] are the (phenomenological) transition probabilities (+ pseudo-counts)
- for a single sequence. It is almost certain (0.99) to stay in a match state, and very
- unlikely (0.05) to go from match to insertion/deletion */
- float zf1SeqTrans[] = {0.98, 0.01, 0.01, 0.25, 0.75, 0.25, 0.75};
- float zf1SeqInit[] = {0.98, 0.01, 0.01, 0.99, 0.01, 0.99, 0.01};
- float zf1SeqDel[] = {0.24, 0.01, 0.75, 0.01, 0.01, 0.01, 0.01};
- float zf1SeqRevrt[] = {0.98, 0.01, 0.01, 0.01, 0.01, 0.99, 0.01};
-
- if ( (NULL == pcPrealigned) || (NULL == pcRepresent) ){
- /*printf("%s/%s:%d: WARNING HMM=NULL -- didn't think I would get here (carry on, no danger)\n",
- __FUNCTION__, __FILE__, __LINE__);*/
- return;
- }
-
- if (NULL == prHMM->p){
- printf("%s:%s:%d: WARNING -- Background Pseudocounts point to NULL\n"
- "\tthis is not intended - don't add background but CONTINUE\n",
- __FUNCTION__, __FILE__, __LINE__);
- return;
-
- }
-
- /* FIXME: should be 0 (FS thinks) but -1 gives better results */
- iH = iS = 0/*-1*//*+1*/;
- while ( ('\0' != *pcS) && ('\0' != *pcH) ){
-
- if ( ('-' != *pcH) && ('-' != *pcS) ){
- iH++;
- iS++;
-
-#if 0
- /* match state
- * - HMM had a gap in previous position (now closed)
- * FIXME: this does not really work */
- if ((iH > 0) && ('-' == *(pcH-1))){
-
- for (iT = 0; iT < STATE_TRANSITIONS; iT++){
- ppfPseudoTr[iS-1][iT] = log2f(zf1SeqRevrt[iT]);
- }
- }
-#endif
-
-#if 1
- /* do frequencies -- this is not really useful;
- frequencies are derived from HMM
- and contain already pseudocounts (PCs),
- adding frequencies and then PCs will add PCs _twice_
- results are better than not to add them,
- but not as good as PCs */
- for (iA = 0; iA < AMINOACIDS; iA++){
- fAux = fFThgiew * ppfFreq[iS][iA] + fFWeight * prHMM->f[iH][iA];
- ppfFreq[iS][iA] = fAux;
- }
-#endif
- /* do pseudo-counts */
- for (iA = 0; iA < AMINOACIDS; iA++){
- fAux =
- fFThgiew * ppfPseudoF[iS][iA] + fFWeight * prHMM->p[iH][iA];
- ppfPseudoF[iS][iA] = fAux;
- }
-#if 0 /* do state transitions */
- for (iT = 0; iT < STATE_TRANSITIONS; iT++){
-#if 1
- /* this is a very crude method,
- which averages the logarithms of the transitions,
- effectively performing a geometric mean -
- this presumably violates normalisation.
- however, results are surprisingly good */
- fAux =
- fGThgiew * ppfPseudoTr[iS][iT] +
- fGWeight * prHMM->tr[iH][iT];
- ppfPseudoTr[iS][iT] = fAux;
-#else /* crude averaging */
- /* There are 2 things to consider:
- (1) one should really blend the probabilities of
- the transitions, however, by default we have
- the logarithms thereof.
- So must exponentiate, blend, and then take log again.
- This is expensive, and does not seem to lead to better
- results than blending the logarithms
- (and violating normalisation)
- (2) transition probabilities for a single sequence
- are really easy, there are no insert/delete transitions.
- However, there is a begin state that is different
- from the main body.
- But again, this does not seem to make a blind bit
- of difference
- */
- if (iS > 0){
- fAux =
- fGThgiew * zf1SeqTrans[iT] +
- fGWeight * prHMM->linTr[iH][iT];
- ppfPseudoTr[iS][iT] = log2f(fAux);
- }
- else {
- fAux =
- fGThgiew * zf1SeqInit[iT] +
- fGWeight * prHMM->linTr[iH][iT];
- ppfPseudoTr[iS][iT] = log2f(fAux);
- }
-#endif /* mixing of linear/log */
- } /* 0 <= iT < STATE_TRANSITIONS */
-#endif /* did state transitions */
- } /* was Match -- neither HMM nor sequence have gap */
-
- else if ('-' == *pcH){
- /* sequence opened up gap in HMM */
-#if 0
- if ((iH > 0) && ('-' != *(pcH-1)) && (iS > 0)){
- /* this is the first gap in HMM
- * FIXME: this does not really work */
- for (iT = 0; iT < STATE_TRANSITIONS; iT++){
- ppfPseudoTr[iS-1][iT] = log2f(zf1SeqDel[iT]);
- }
- }
- else {
- /* do nothing, keep single sequence values exactly as they are*/
- }
-#endif
- iS++;
- }
- else if ('-' == *pcS){
- /* here the single sequence has a gap,
- and the HMM (as a whole) does not. There may be individual gaps
- in the HMM at this stage. By ignoring this we say that the HMM
- dominates the overall behaviour - as in the M2M state as well
- */
- iH++;
- }
-
- pcH++;
- pcS++;
-
- } /* !EO[seq/hmm] */
-
- return;
-
-} /* this is the end of AddBackgroundFrequencies() */
-
-
-
-/**
- * ReadAndPrepare()
- *
- * @brief Read HMM input file or transfer alignment
- * and add pseudocounts etc.
- *
- * @param[in] iRnPtype
- * type of read/prepare
- * enum {INTERN_ALN_2_HMM = 0, READ_ALN_2_HMM, READ_HMM_2_HMM};
- * @param[in] ppcProf
- * alignment
- * @param[in] iCnt
- * number of seqs in alignment
- * @param[in,out] prHMM,
- * [in] if sequences read/prepared, [out] if HMM from file
- * @param[in] pcPrealigned,
- * (single) sequence aligned to background HMM
- * @param[in] pcRepresent,
- * sequence representing HMM aligned to individual sequence
- * param[in] pdExWeights
- * (external) sequence weights, derived from tree
- * @param[in] infile
- * name of file with HMM info (formerly also alignment)
- * @param[out] q
- * HMM structure with transition probabilities, residue frequencies etc
- * @param[???] qali
- * FIXME: what is qali?
- *
- * @return FAILURE on error, OK otherwise
- */
-int
-ReadAndPrepare(int iRnPtype,
- char **ppcProf, int iCnt, hmm_light *prHMM,
- char *pcPrealigned, char *pcRepresent, double *pdExWeights,
- char* infile, HMM& q, Alignment* qali=NULL) {
-
- //#ifndef CLUSTALO_NOFILE
- char path[NAMELEN];
-
- /* NOTE: there are different scenarios:
-
- (i) ("" != infile) - read HMM from file,
- transfer frequency/transition/pseudo-count (f/t/p) info into prHMM
-
- (ii) ('\0' != ppcProf[0]) - transfer sequence/alignment into q/qali,
- don't save f/t/p into prHMM,
- on the contrary, if prior f/t/p available then add it to q/qali,
- this is only done if (1==iCnt)
-
- (iii) ('\0' == ppcProf[0]) - re-cycle old HMM information
- recreate a HMM from previous data
- */
-
- /********************************/
- /*** (o) Recycle internal HMM ***/
- /********************************/
- if ( (INTERN_ALN_2_HMM == iRnPtype) && (iCnt <= 0) ){
-
- /* NOTE: here we are writing into a HMM structure/class;
- memory has been allocated for this in hhalign.cpp;
- however, as iCnt<=0, there may not be memory for
- prHMM->n_display sequences/names.
- But then, there doesn't have to be.
- At this stage we are just copying one HMM into another HMM,
- sequences are irrelevant. The only sequence of (marginal)
- interest is the consensus sequence */
- /* FIXME: check that prHMM is valid -- how? */
-
- const int ciCons = 0;
- const int ciNoof = ciCons+1;
- const int ciInvalid = -1;
- q.n_display = ciNoof; /* only consensus */
- q.sname = NULL;
- q.ncons = ciCons;
- q.nfirst = ciCons;//prHMM->nfirst;
- q.nss_dssp = ciInvalid;//prHMM->nss_dssp;
- q.nsa_dssp = ciInvalid;//prHMM->nsa_dssp;
- q.nss_pred = ciInvalid;//prHMM->nss_pred;
- q.nss_conf = ciInvalid;//prHMM->nss_conf;
- q.L = prHMM->L;
- q.N_in = prHMM->N_in;
- q.N_filtered = prHMM->N_filtered;
- /* NOTE: I (FS) think that only ever 1 sequence will be transferred here,
- that is, the consensus sequence. However, we might want to allow
- (in the future) to transfer more sequences,
- hence the awkward for() loop */
-#if 0
- for (int i = prHMM->ncons+0; i < prHMM->ncons+q.n_display; i++){
- /* NOTE: In the original hhalign code the first position
- is kept open ('\0'). This makes it difficult to use
- string functions like strlen/strdup etc.
- Insert a temporary gap (.) to facilitate string operations */
- char cHead = prHMM->seq[i][0];
- prHMM->seq[i][0] = '.';
- q.seq[i] = strdup(prHMM->seq[i]);
- prHMM->seq[i][0] = q.seq[i][0] = cHead;
- }
-#else
- {
- char cHead = prHMM->seq[prHMM->ncons][0];
- prHMM->seq[prHMM->ncons][0] = '.';
- q.seq[q.ncons] = strdup(prHMM->seq[prHMM->ncons]);
- prHMM->seq[prHMM->ncons][0] = q.seq[q.ncons][0] = cHead;
- }
-#endif
- const int NEFF_SCORE = 10; /* FIXME: Magic Number */
- for (int i = 0; i < prHMM->L+1; i++){
- q.Neff_M[i] = prHMM->Neff_M[i];
- q.Neff_I[i] = prHMM->Neff_I[i];
- q.Neff_D[i] = prHMM->Neff_D[i];
- }
- q.Neff_HMM = prHMM->Neff_HMM;
- /* skip longname,name,file,fam,sfam,fold,cl */
- q.lamda = prHMM->lamda;
- q.mu = prHMM->mu;
-
- HMMshadow rShadow = {0}; /* friend of HMM to access private members */
- rShadow.copyShadowToHMM(q, *prHMM);
-
- /* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
- /* pav already done in copyShadowToHMM */
- /* skip pnul */
-
- return OK;
-
-
- } /* INTERN_ALN_2_HMM && iCnt<=0 */
-
- /******************************/
- /*** (i) Read HMM from file ***/
- /******************************/
- char line[LINELEN] = {0}; // input line
- FILE *inf = NULL;
- //if ( (0 != strcmp(infile,"")) /*&& (iCnt > 0)*/ )
- if ( (READ_HMM_2_HMM == iRnPtype) || (READ_ALN_2_HMM == iRnPtype)) {
-
- if (0 == strcmp(infile,"")){
- printf("%s:%s:%d:\n"
- "\texpected to re %s from file but no file specified\n"
- ""
- , __FUNCTION__, __FILE__, __LINE__
- , (READ_HMM_2_HMM==iRnPtype?"HMM":"alignment"));
- return FAILURE;
- }
- inf = fopen(infile, "r");
- if (!inf) OpenFileError(infile);
- Pathname(path,infile);
-
- //}
- //else {
- //inf = stdin;
- //if (v>=2) printf("Reading HMM / multiple alignment from standard input ...\n(To get a help list instead, quit and type %s -h.)\n",program_name);
- //*path='\0';
- //}
-
- fgetline(line,LINELEN-1,inf);
- }
-
- //if ( (0 != strcmp(infile,"")) && (iCnt > 0) )
- if ( (READ_HMM_2_HMM == iRnPtype) ){
-
- if (0 == strcmp(infile,"")){
- printf("%s:%s:%d: expected to read HMM from file but no file-name\n",
- __FUNCTION__, __FILE__, __LINE__);
- }
-
- // Is infile a HMMER3 file? /* MR1 */
- if (!strncmp(line,"HMMER3",6))
- {
- if (v>=2) {
- cout<<"Query file is in HMMER3 format\n";
- cout<<"WARNING! Use of HMMER3 format as input results in dramatically loss of sensitivity!\n";
- }
-
- // Read 'query HMMER file
- rewind(inf);
- q.ReadHMMer3(inf,path);
-
- // Don't add transition pseudocounts to query!!
-
- // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
- q.PreparePseudocounts();
-
- // DON'T ADD amino acid pseudocounts to query: pcm=0! q.p[i][a] = f[i][a]
- q.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
-
- /* further down there is a q.Log2LinTransitionProbs
- but only if (iCnt>0), however, we still need it it here (i think),
- there is no danger of doing this twice, as trans_lin is checked
- FIXME (FS, 2011-01-12) */
- /* further down there is a q.Log2LinTransitionProbs
- but only if (iCnt>0), however, we still need it it here (i think),
- there is no danger of doing this twice, as trans_lin is checked
- FIXME (FS, 2011-01-12) */
- //if (par.forward >= 1)
- {
- q.Log2LinTransitionProbs(1.0);
- }
-
- }
- // ... or Is infile an old HMMER file?
- else if (!strncmp(line,"HMMER",5)) {
- if (v>=2) {
- cout<<"Query file is in HMMER format\n";
- cout<<"WARNING! Use of HMMER format as input results in dramatically loss of sensitivity!\n";
- }
-
- // Read 'query HMMER file
- q.ReadHMMer(inf,path);
- if (v>=2 && q.Neff_HMM>11.0)
- fprintf(stderr,"WARNING: HMM %s looks too diverse (Neff=%.1f>11). Better check the underlying alignment... \n",q.name,q.Neff_HMM);
-
- // Don't add transition pseudocounts to query!!
-
- // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
- q.PreparePseudocounts();
-
- // DON'T ADD amino acid pseudocounts to query: pcm=0! q.p[i][a] = f[i][a]
- q.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
-
- /* further down there is a q.Log2LinTransitionProbs
- but only if (iCnt>0), however, we still need it it here (i think),
- there is no danger of doing this twice, as trans_lin is checked
- FIXME (FS, 2011-01-12) */
- //if (par.forward >= 1)
- {
- q.Log2LinTransitionProbs(1.0);
- }
-
- } /* it was a HMMer file */
-
- // ... or is it an hhm file?
- else if (!strncmp(line,"NAME",4) || !strncmp(line,"HH",2)) {
-
- if (v>=2) cout<<"Query file is in HHM format\n";
-
- // Rewind to beginning of line and read query hhm file
- rewind(inf);
- q.Read(inf,path);
- if (v>=2 && q.Neff_HMM>11.0)
- fprintf(stderr,"WARNING: HMM %s looks too diverse (Neff=%.1f>11). Better check the underlying alignment... \n",q.name,q.Neff_HMM);
-
- // Add transition pseudocounts to query -> q.p[i][a]
- q.AddTransitionPseudocounts();
-
- // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
- q.PreparePseudocounts();
-
- // Add amino acid pseudocounts to query: q.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
- q.AddAminoAcidPseudocounts();
-
- } /* it was a HHM file */
- else {
- fprintf(stderr, "%s:%s:%d: Unknown HMM format in infile\n"
- "infile=%s, #seq=%d\n"
- , __FUNCTION__, __FILE__, __LINE__
- , infile, iCnt);
- return FAILURE;
- }
-
- /*fclose(inf);*/
-
- /*** transfer class info to struct */
- prHMM->n_display = q.n_display;
- /* ignore sname*/
- prHMM->seq = (char **)calloc((q.n_display+1), sizeof(char *));
- /* FIXME valgrind says bytes get lost in the above calloc during
- * hmm-iteration
- */
- for (int i = 0; i < q.n_display; i++){
- /* NOTE: In the original hhalign code the first position
- is kept open ('\0'). This makes it difficult to use
- string functions like strlen/strdup etc.
- Insert a temporary gap (.) to facilitate string operations */
- char cHead = q.seq[i][0];
- q.seq[i][0] = '.';
- prHMM->seq[i] = strdup(q.seq[i]);
- q.seq[i][0] = prHMM->seq[i][0] = cHead;
- }
- prHMM->ncons = q.ncons;
- prHMM->nfirst = q.nfirst;
- prHMM->nss_dssp = q.nss_dssp;
- prHMM->nsa_dssp = q.nsa_dssp;
- prHMM->nss_pred = q.nss_pred;
- prHMM->nss_conf = q.nss_conf;
- prHMM->L = q.L;
- prHMM->N_in = q.N_in;
- prHMM->N_filtered = q.N_filtered;
- const int NEFF_SCORE = 10; /* FIXME: Magic Number */
- prHMM->Neff_M = (float *)calloc(prHMM->L+1, sizeof(float));
- prHMM->Neff_I = (float *)calloc(prHMM->L+1, sizeof(float));
- prHMM->Neff_D = (float *)calloc(prHMM->L+1, sizeof(float));
- /*for (int i = 0; i < prHMM->L+1; i++){
- prHMM->Neff_M[i] = prHMM->Neff_I[i] = prHMM->Neff_D[i] = NEFF_SCORE;
- }*/
- prHMM->Neff_HMM = q.Neff_HMM;
- /* skip longname,name,file,fam,sfam,fold,cl */
- prHMM->lamda = q.lamda;
- prHMM->mu = q.mu;
-
- HMMshadow rShadow = {0}; /* friend of HMM to access private members */
- rShadow.copyHMMtoShadow(q);
-
- prHMM->f = (float **)calloc(prHMM->L+1, sizeof(float *));
- prHMM->g = (float **)calloc(prHMM->L+1, sizeof(float *));
- prHMM->p = (float **)calloc(prHMM->L+1, sizeof(float *));
- prHMM->tr = (float **)calloc(prHMM->L+1, sizeof(float *));
- prHMM->linTr = (float **)calloc(prHMM->L+1, sizeof(float *));
- for (int i = 0; i < prHMM->L+1; i++){
- prHMM->f[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- prHMM->g[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- prHMM->p[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- for (int j = 0; j < AMINOACIDS; j++){
- prHMM->f[i][j] = (float)rShadow.f[i][j];
- prHMM->g[i][j] = (float)rShadow.g[i][j];
- prHMM->p[i][j] = (float)rShadow.p[i][j];
- }
- prHMM->tr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
- prHMM->linTr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
- for (int j = 0; j< STATE_TRANSITIONS; j++){
- prHMM->tr[i][j] = (float)rShadow.tr[i][j];
- prHMM->linTr[i][j] = fpow2(rShadow.tr[i][j]);
- }
- } /*0 <= i < prHMM->L+1 */
- /* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
- for (int j = 0; j < AMINOACIDS; j++){
- prHMM->pav[j] = (float)rShadow.pav[j];
- }
- /* skip pnul */
-
- } /* have read HHM from file */
- /*else if ( ((NULL != ppcProf) && (iCnt > 0) && ('\0' != ppcProf[0][0])) ||
- ( (0 != strcmp(infile,"") && (0 == iCnt) )) )*/
- else if ( (INTERN_ALN_2_HMM == iRnPtype) || (READ_ALN_2_HMM == iRnPtype) ) {
-
- if ( (INTERN_ALN_2_HMM == iRnPtype) &&
- ( (NULL == ppcProf) || (iCnt <= 0) || ('\0' == ppcProf[0][0]) ) ){
- printf("%s:%s:%d: was expecting internal alignment but\n"
- "\tppcProf=%p, #seq=%d, ppcProf[0][0]=%c\n"
- , __FUNCTION__, __FILE__, __LINE__
- , ppcProf, iCnt, ppcProf[0][0]);
- return FAILURE;
- }
- else if ( (READ_ALN_2_HMM == iRnPtype) &&
- (0 == strcmp(infile,"")) ){
- printf("%s:%s:%d: was expecting to read alignment from file but no filename\n"
- , __FUNCTION__, __FILE__, __LINE__);
- return FAILURE;
- }
-
- /*******************************/
- /*** (ii) it is an alignment ***/
- /*******************************/
- /* transfer alignment information from clustal character array
- into pali/q/t classes */
- /*
- * NOTE that emissions in HMMer file format contain pseudo-counts.
- * HHM file format does not contain emission pseudo-counts.
- * the structure that stores background HMM information does contain pseudo-counts
- */
- Alignment* pali;
- if (qali==NULL){
- pali=new Alignment(iCnt); /* FIXME: pass in iCnt to get rid of MAXSEQ */
- }
- else{
- pali=qali;
- }
-
- if (par.calibrate) {
- printf("\nError in %s: only HHM files can be calibrated.\n",program_name);
- printf("Build an HHM file from your alignment with 'hhmake -i %s' and rerun hhsearch with the hhm file\n\n",infile);
- exit(1);
- }
-
- if (v>=2 && strcmp(infile,"stdin")) cout<<infile<<" is in A2M, A3M or FASTA format\n";
-
- /* Read alignment from infile into matrix X[k][l] as ASCII
- (and supply first line as extra argument) */
- //if (iCnt > 0)
- if (INTERN_ALN_2_HMM == iRnPtype){
- pali->Transfer(ppcProf, iCnt);
- }
- //else if (0 == iCnt)
- else if (READ_ALN_2_HMM == iRnPtype){
- pali->Read(inf, infile, line);
- }
- else {
- printf("%s:%s:%d: FATAL problem\n"
- "infile = (%s), #seq = %d\n"
- , __FUNCTION__, __FILE__, __LINE__
- , infile, iCnt);
- return FAILURE;
- }
-
- /* Convert ASCII to int (0-20), throw out all insert states,
- record their number in I[k][i]
- and store marked sequences in name[k] and seq[k] */
- pali->Compress(infile);
-
- /* Sort out the nseqdis most dissimilar sequences for display
- in the output alignments */
- pali->FilterForDisplay(par.max_seqid, par.coverage, par.qid,
- par.qsc,par.nseqdis);
-
- // Filter alignment for min score per column with core query profile, defined by coverage_core and qsc_core
- //if (par.coresc>-10) pali->HomologyFilter(par.coverage_core, par.qsc_core, par.coresc);
-
- /* Remove sequences with seq. identity larger than seqid percent
- (remove the shorter of two) */
- pali->N_filtered = pali->Filter(par.max_seqid, par.coverage,
- par.qid, par.qsc, par.Ndiff);
-
- /* Calculate pos-specific weights,
- AA frequencies and transitions -> f[i][a], tr[i][a] */
- pali->FrequenciesAndTransitions(q);
- if (v>=2 && q.Neff_HMM>11.0)
- fprintf(stderr,"WARNING: alignment %s looks too diverse (Neff=%.1f>11). Better check it with an alignment viewer... \n",q.name,q.Neff_HMM);
-
- /*printf("%d %d %f %d (N,Nf,Neff,L) %s:%s:%d\n"
- , q.N_in, q.N_filtered, q.Neff_HMM, q.L, __FUNCTION__, __FILE__, __LINE__);*/
-
- // Add transition pseudocounts to query -> p[i][a]
- q.AddTransitionPseudocounts();
-
- /* Generate an amino acid frequency matrix from f[i][a]
- with full pseudocount admixture (tau=1) -> g[i][a] */
- q.PreparePseudocounts();
-
- /* Add amino acid pseudocounts to query:
- p[i][a] = (1-tau)*f[i][a] + tau*g[i][a] */
- q.AddAminoAcidPseudocounts();
-
- /* ****** add aligned background pseudocounts ***** */
- HMMshadow rShadowQ = {0};
- rShadowQ.copyHMMtoShadow(q);
-
- AddBackgroundFrequencies(rShadowQ.f, rShadowQ.p, rShadowQ.tr,
- q.L, prHMM,
- q.seq, pcPrealigned, iCnt, pcRepresent);
-
-
- if (qali==NULL){
- delete(pali); pali = NULL;
- }
-
- } /* else if (NULL != ppcProf) // not hmmer/hhalign but alignment */
-
- //else if ((prHMM->L > 0) && ('\0' == ppcProf[0][0]))
- else if (INTERN_HMM_2_HMM == iRnPtype){
-
- /******************************************/
- /*** (iii) re-cycle old HMM information ***/
- /******************************************/
-
- if (prHMM->L <= 0){
- printf("%s:%s:%d: was expecting to copy HMM structure but L=%d\n"
- , __FUNCTION__, __FILE__, __LINE__, prHMM->L);
- }
-
- printf("%s:%s:%d: RE-CYCLING HMM\n", __FUNCTION__, __FILE__, __LINE__);
-
-#if 0
- q.n_display = prHMM->n_display;
- /* ignore sname*/
- for (int i = 0; i < q.n_display; i++){
- /* NOTE: In the original hhalign code the first position
- is kept open ('\0'). This makes it difficult to use
- string functions like strlen/strdup etc.
- Insert a temporary gap (.) to facilitate string operations */
- char cHead = prHMM->seq[i][0];
- prHMM->seq[i][0] = '.';
- q.seq[i] = strdup(prHMM->seq[i]);
- q.seq[i][0] = prHMM->seq[i][0] = cHead;
- }
- q.nfirst = prHMM->nfirst;
-#else
- q.n_display = 1;
- q.nfirst = 0;
- char cHead = prHMM->seq[prHMM->nfirst][0];
- prHMM->seq[prHMM->nfirst][0] = '.';
- q.seq[0] = strdup(prHMM->seq[prHMM->nfirst]);
- q.seq[q.n_display-1][0] = prHMM->seq[prHMM->nfirst][0] = cHead;
-#endif
- q.ncons = prHMM->ncons;
- q.nss_dssp = prHMM->nss_dssp;
- q.nsa_dssp = prHMM->nsa_dssp;
- q.nss_pred = prHMM->nss_pred;
- q.nss_conf = prHMM->nss_conf;
- q.L = prHMM->L;
- q.N_in = prHMM->N_in;
- q.N_filtered = prHMM->N_filtered;
-#define NEFF_SCORE 10 /* FIXME: Magic Number */
- /*for (int i; i < prHMM->L+1; i++){
- q.Neff_M[i] = q.Neff_I[i] = q.Neff_D[i] = NEFF_SCORE;
- }*/
- q.Neff_HMM = prHMM->Neff_HMM;
- /* skip longname,name,file,fam,sfam,fold,cl */
- q.lamda = prHMM->lamda;
- q.mu = prHMM->mu;
-
- HMMshadow rShadow = {0}; /* friend of HMM to access private members */
- rShadow.copyShadowToHMM(q, *prHMM);
-
- }
-
- if (iCnt > 0){
- if (par.forward>=1) q.Log2LinTransitionProbs(1.0);
- }
-
- if (NULL != inf){
- fclose(inf);
- }
-
- return OK;
-
-} /*** end: ReadAndPrepare() ***/
-
-/**
- * FreeHMMstruct()
- *
- * @brief FIXME
- *
- * @param[in,out]
- */
-extern "C" void
-FreeHMMstruct(hmm_light *prHMM)
-{
- int i;
-
- if (NULL != prHMM->f){
- for (i = 0; i < prHMM->L+1; i++){
- if (NULL != prHMM->f[i]){
- free(prHMM->f[i]); prHMM->f[i] = NULL;
- }
- } /* 0 <= i < prHMM->L+1 */
- free(prHMM->f); prHMM->f = NULL;
- }
- if (NULL != prHMM->g){
- for (i = 0; i < prHMM->L+1; i++){
- if (NULL != prHMM->g[i]){
- free(prHMM->g[i]); prHMM->g[i] = NULL;
- }
- } /* 0 <= i < prHMM->L+1 */
- free(prHMM->g); prHMM->g = NULL;
- }
- if (NULL != prHMM->p){
- for (i = 0; i < prHMM->L+1; i++){
- if (NULL != prHMM->p[i]){
- free(prHMM->p[i]); prHMM->p[i] = NULL;
- }
- } /* 0 <= i < prHMM->L+1 */
- free(prHMM->p); prHMM->p = NULL;
- }
- if (NULL != prHMM->tr){
- for (i = 0; i < prHMM->L+1; i++){
- if (NULL != prHMM->tr[i]){
- free(prHMM->tr[i]); prHMM->tr[i] = NULL;
- }
- } /* 0 <= i < prHMM->L+1 */
- free(prHMM->tr); prHMM->tr = NULL;
- }
- if (NULL != prHMM->linTr){
- for (i = 0; i < prHMM->L+1; i++){
- if (NULL != prHMM->linTr[i]){
- free(prHMM->linTr[i]); prHMM->linTr[i] = NULL;
- }
- } /* 0 <= i < prHMM->L+1 */
- free(prHMM->linTr); prHMM->linTr = NULL;
- }
-
- if (NULL != prHMM->Neff_M){
- free(prHMM->Neff_M); prHMM->Neff_M = NULL;
- }
- if (NULL != prHMM->Neff_I){
- free(prHMM->Neff_I); prHMM->Neff_I = NULL;
- }
- if (NULL != prHMM->Neff_D){
- free(prHMM->Neff_D); prHMM->Neff_D = NULL;
- }
-
- if (NULL != prHMM->seq){
- for (i = 0; i < prHMM->n_display; i++){
- if (NULL != prHMM->seq[i]){
- free(prHMM->seq[i]); prHMM->seq[i] = NULL;
- }
- }
- free(prHMM->seq); prHMM->seq = NULL;
- }
-
- memset(prHMM, 0, sizeof(hmm_light));
-
-} /*** end: FreeHMMstruct() ***/
-
-
-/**
- * @brief comparisin function for ascending qsort() of doubles
- */
-int
-CompDblAsc(const void *pv1, const void *pv2){
-
- double d1 = *(double *)pv1;
- double d2 = *(double *)pv2;
-
- if (d1 > d2) { return +1; }
- else if (d1 < d2) { return -1; }
- else { return 0; }
-
-} /*** end: CompDblAsc() ***/
-
-
-/**
- * AlnToHMM2()
- *
- * @brief convert alignment into HMM (hhmake)
- *
- * @param[out] prHMM
- * structure with pseudocounts etc
- * @param[in] pcHMM_input
- * name of file with HMM
- *
- */
-extern "C" int
-AlnToHMM2(hmm_light *rHMM_p,
- char **ppcSeq, int iN){
-
- if (0 == par.M){
- SetDefaults();
- SetSubstitutionMatrix();
- par.cons = 1;
- par.M = 2;
- par.forward=2;
- par.Mgaps=100;
- const int ciGoodMeasureSeq = 10;
- int iAuxLen = strlen(ppcSeq[0]);
- par.maxResLen = iAuxLen;
- par.maxResLen += ciGoodMeasureSeq;
- par.maxColCnt = iAuxLen + ciGoodMeasureSeq;
- par.max_seqid=DEFAULT_FILTER;
- par.loc=0; par.mact=0;
- /* some minor parameters, affecting alignment (i think) */
- par.p = 0.0; /* minimum threshold for inclusion in hit list */
- par.Z = 100; /* minimum threshold for inclusion in hit list and alignment listing */
- par.z = 1; /* min number of lines in hit list */
- par.B = 100; /* max number of alignments */
- par.b = 1; /* min number of alignments */
- par.E = 1e6; /* maximum threshold for inclusion in hit list and alignment listing */
- par.altali=1;par.hitrank=0;par.showcons=1; par.showdssp=1;par.showpred=1;par.nseqdis=iN;par.cons=1;
- }
-
- const int ciGoodMeasure = 10;
- int iLen = strlen(ppcSeq[0]) + ciGoodMeasure;
- if (iLen > par.maxResLen){
- par.maxResLen = par.maxColCnt = iLen;
- }
- HMM rTemp(iN+2, iLen); /* temporary template */
- Alignment rTempAli(iN+2, iLen); /* temporary alignment */
- int iParCons = par.cons;
-
- /*printf(">>>>>>>>>>> %s:%s:%d: there are %d sequences\n", __FUNCTION__, __FILE__, __LINE__, iN);*/
-
- par.cons = 1;
- if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
- ppcSeq, iN, rHMM_p,
- NULL/*prealigned*/, NULL/*representative*/, NULL/*weights*/, //YES/*transfer*/,
- (char*)("")/*in-file*/, rTemp, &rTempAli)) {
- return FAILURE;
- }
- par.cons = iParCons;
-
- /******/
- /*** transfer class info to struct */
- rHMM_p->n_display = rTemp.n_display;
- rHMM_p->sname = NULL;
- rHMM_p->seq = (char **)calloc((rTemp.n_display+1), sizeof(char *));
-
- for (int i = 0; i < rTemp.n_display; i++){
- /* NOTE: In the original hhalign code the first position
- is kept open ('\0'). This makes it difficult to use
- string functions like strlen/strdup etc.
- Insert a temporary gap (.) to facilitate string operations */
- char cHead = rTemp.seq[i][0];
- rTemp.seq[i][0] = '.';
- rHMM_p->seq[i] = strdup(rTemp.seq[i]);
- rTemp.seq[i][0] = rHMM_p->seq[i][0] = cHead;
- }
- rHMM_p->ncons = rTemp.ncons;
- rHMM_p->nfirst = rTemp.nfirst;
- if (-1 == rHMM_p->ncons){
- /* ncons needed later for alignment of
- representative sequence and copy of profile.
- ncons not always set (-1 default),
- this will cause segmentation fault.
- nfirst (probably) right index -
- no problem if not */
- rHMM_p->ncons = rTemp.nfirst;
- }
- rHMM_p->nss_dssp = rTemp.nss_dssp;
- rHMM_p->nsa_dssp = rTemp.nsa_dssp;
- rHMM_p->nss_pred = rTemp.nss_pred;
- rHMM_p->nss_conf = rTemp.nss_conf;
- rHMM_p->L = rTemp.L;
- rHMM_p->N_in = rTemp.N_in;
- rHMM_p->N_filtered = rTemp.N_filtered;
-#define NEFF_SCORE 10 /* FIXME: Magic Number */ //// get rid of that, FS, 2010-12-22
- rHMM_p->Neff_M = (float *)calloc(rHMM_p->L+1, sizeof(float));
- rHMM_p->Neff_I = (float *)calloc(rHMM_p->L+1, sizeof(float));
- rHMM_p->Neff_D = (float *)calloc(rHMM_p->L+1, sizeof(float));
- for (int i = 0; i < rHMM_p->L+1; i++){
- rHMM_p->Neff_M[i] = rHMM_p->Neff_I[i] = rHMM_p->Neff_D[i] = NEFF_SCORE; //// get rid of that, FS, 2010-12-22
- }
- rHMM_p->Neff_HMM = rTemp.Neff_HMM;
- /* skip longname,name,file,fam,sfam,fold,cl */
- rHMM_p->lamda = rTemp.lamda;
- rHMM_p->mu = rTemp.mu;
-
- HMMshadow rShadow = {0}; /* friend of HMM to access private members */
- rShadow.copyHMMtoShadow(rTemp);
-
- rHMM_p->Neff_M = (float *)calloc(rHMM_p->L+1, sizeof(float));
- rHMM_p->Neff_I = (float *)calloc(rHMM_p->L+1, sizeof(float));
- rHMM_p->Neff_D = (float *)calloc(rHMM_p->L+1, sizeof(float));
- rHMM_p->f = (float **)calloc(rHMM_p->L+1, sizeof(float *));
- rHMM_p->g = (float **)calloc(rHMM_p->L+1, sizeof(float *));
- rHMM_p->p = (float **)calloc(rHMM_p->L+1, sizeof(float *));
- rHMM_p->tr = (float **)calloc(rHMM_p->L+1, sizeof(float *));
- rHMM_p->linTr = (float **)calloc(rHMM_p->L+1, sizeof(float *));
-
- for (int i = 0; i < rHMM_p->L+1; i++){
- rHMM_p->Neff_M[i] = (float)rShadow.Neff_M[i];
- rHMM_p->Neff_I[i] = (float)rShadow.Neff_I[i];
- rHMM_p->Neff_D[i] = (float)rShadow.Neff_D[i];
- rHMM_p->f[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- rHMM_p->g[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- rHMM_p->p[i] = (float *)calloc(AMINOACIDS, sizeof(float));
- for (int j = 0; j < AMINOACIDS; j++){
- rHMM_p->f[i][j] = (float)rShadow.f[i][j];
- rHMM_p->g[i][j] = (float)rShadow.g[i][j];
- rHMM_p->p[i][j] = (float)rShadow.p[i][j];
- }
- rHMM_p->tr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
- rHMM_p->linTr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
- for (int j = 0; j< STATE_TRANSITIONS; j++){
- rHMM_p->tr[i][j] = (float)rShadow.tr[i][j];
- rHMM_p->linTr[i][j] = fpow2(rShadow.tr[i][j]);
- }
- } /*0 <= i < prHMM->L+1 */
- /* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
- for (int j = 0; j < AMINOACIDS; j++){
- rHMM_p->pav[j] = (float)rShadow.pav[j];
- }
- /* skip pnul */
-
-
- /******/
-
-
- rTemp.ClobberGlobal();
- rTempAli.ClobberGlobal();
-
- return OK;
-
-} /*** end of AlnToHMM2() ***/
-
-
-/**
- * HHMake_Wrapper()
- *
- * @brief turn alignment (from file) into HHM (HMM) on file
- *
- * @param[out] hmm_out
- * name of file with HMM info corresponding to tmp_aln
- * @param[in] tmp_aln
- * name of file with alignment, to be turned into HMM (HHM)
- *
- */
-extern "C" int
-HHMake_Wrapper(char *tmp_aln, char *hmm_out)
-{
-
- HMM rTemp; /* temporary template */
- Alignment rTempAli; /* temporary alignment */
- hmm_light *rHMM_p = NULL;
-
- /** Note:
- this is a wrapper for a stand-alone program hhmake.
- hhmake uses a different set of parameters than hhalign.
- However, all parameters are GLOBAL.
- So at this point we register the hhalign settings,
- reset them to hhmake settings and set them back
- at the end of the function
- */
-
- /* save settings from hhalign */
- int iParShowcons=par.showcons;
- int iParAppend=par.append;
- int iParNseqdis=par.nseqdis;
- int iParMark=par.mark;
- int iParMaxSeqid=par.max_seqid;
- int iParQid=par.qid;
- double dParQsc=par.qsc;
- int iParCoverage=par.coverage;
- int iParNdiff=par.Ndiff;
- int iParCoverageCore=par.coverage_core;
- double dParQscCore=par.qsc_core;
- double dParCoresc=par.coresc;
- int iParM=par.M;
- int iParMgaps=par.Mgaps;
- int iParMatrix=par.matrix;
- int iParPcm=par.pcm;
- double dParPcw=par.pcw;
- double dParGapb=par.gapb;
- int iParWg=par.wg;
-
- /* these are settings suitable for hhmake */
- par.showcons=1; // write consensus sequence into hhm file
- par.append=0; // overwrite output file
- par.nseqdis=10; // maximum number of query or template sequences to be recoreded in HMM and diplayed in output alignments
- par.mark=0; // 1: only marked sequences (or first) get displayed; 0: most divergent ones get displayed
- par.max_seqid=90; // default for maximum sequence identity threshold
- par.qid=0; // default for maximum sequence identity threshold
- par.qsc=-20.0f; // default for minimum score per column with query
- par.coverage=0; // default for minimum coverage threshold
- par.Ndiff=100; // pick Ndiff most different sequences from alignment
- par.coverage_core=30; // Minimum coverage for sequences in core alignment
- par.qsc_core=0.3f; // Minimum score per column with core alignment (HMM)
- par.coresc=-20.0f; // Minimum score per column with core alignment (HMM)
- par.M=1; // match state assignment is by A2M/A3M
- par.Mgaps=50; // above this percentage of gaps, columns are assigned to insert states
- par.matrix=0; // Subst.matrix 0: Gonnet, 1: HSDM, 2: BLOSUM50 3: BLOSUM62
- par.pcm=0; // no amino acid and transition pseudocounts added
- par.pcw=0; // wc>0 weighs columns according to their intra-clomun similarity
- par.gapb=0.0; // default values for transition pseudocounts; 0.0: add no transition pseudocounts!
- par.wg=0; // 0: use local sequence weights 1: use local ones
-
-
- if (OK != ReadAndPrepare(READ_ALN_2_HMM,
- NULL, 0, rHMM_p, NULL, NULL, NULL,
- tmp_aln, rTemp, &rTempAli)) {
- return FAILURE;
- }
-
- rTemp.WriteToFile(hmm_out);
-
-
-
- /* restore settings from hhalign */
- par.showcons=iParShowcons;
- par.append=iParAppend;
- par.nseqdis=iParNseqdis;
- par.mark=iParMark;
- par.max_seqid=iParMaxSeqid;
- par.qid=iParQid;
- par.qsc=dParQsc;
- par.coverage=iParCoverage;
- par.Ndiff=iParNdiff;
- par.coverage_core=iParCoverageCore;
- par.qsc_core=dParQscCore;
- par.coresc=dParCoresc;
- par.M=iParM;
- par.Mgaps=iParMgaps;
- par.matrix=iParMatrix;
- par.pcm=iParPcm;
- par.pcw=dParPcw;
- par.gapb=dParGapb;
- par.wg=iParWg;
-
-
- /* prepare exit */
- rTemp.ClobberGlobal();
- rTempAli.ClobberGlobal();
-
- return OK;
-}
-
-/**
- * readHMMWrapper()
- *
- * @brief read HMM from file, copy (relevant) info into struct
- *
- * @param[out] prHMM
- * structure with pseudocounts etc, probably uninitialised on entry
- * @param[in] pcHMM_input
- * name of file with HMM
- *
- */
-extern "C" int
-readHMMWrapper(hmm_light *rHMM_p,
- char *pcHMM_input)
-{
- par.maxResLen = 15002;
- HMM rTemp(1000,par.maxResLen); /* temporary template */
- Alignment rTempAli; /* temporary alignment */
- /* AW changed init from {0} to 0 because it failed to compile with
- * gcc 4.3.3 with the following error:
- * error: braces around initializer for non-aggregate type
- */
- /* FS taken out initialiser alltogether */
-
- /* 0th arg of RnP is the type of RnP, ie,
- enum {INTERN_ALN_2_HMM = 0, READ_ALN_2_HMM, READ_HMM_2_HMM};*/
- /* 1st arg of ReadAndPrepare() is ppcSeqs, 2nd is #seq */
- /* FIXME: at this stage the 3rd arg, rHMM_p, only acts as a dummy.
- this is rather silly, as it is the actual struct that will
- carry the HMM info at the end.
- If we write it already in ReadAndPrepare() then we could
- dispense with all this friend-class nonsense */
- if (OK != ReadAndPrepare(READ_HMM_2_HMM,
- NULL, 0, rHMM_p, NULL, NULL, NULL,
- pcHMM_input, rTemp, &rTempAli)) {
- return FAILURE;
- }
-
- /* an important piece of information I want to get out of here
- is the consenssus sequence. there are however certain
- Pfam HMMs that don't trigger consensus calculation.
- at the moment I (FS) don't understand why this is
- (or rather why this is not). The proper place to do this
- should be inside ReadAndPrepare():ReadHMMer(), but
- I have not quite penetrated the logic there.
- Therefore I try to catch this condition at this point (here)
- and rectify it.
- */
- if (-1 == rHMM_p->ncons){
- int i, iA;
- rHMM_p->ncons = rHMM_p->nfirst;
-
- for (i = 0; i < rHMM_p->L; i++){
- double dPmax = 0.00;
- int iAmax = -1;
- for (iA = 0; iA < AMINOACIDS; iA++){
- if (rHMM_p->f[i][iA] > dPmax){
- iAmax = iA;
- dPmax = rHMM_p->f[i][iA];
- }
- } /* (0 <= iA < AMINOACIDS) */
- rHMM_p->seq[rHMM_p->ncons][i] = i2aa(iAmax);
- } /* (0 <= i < rHMM_p->L) */
-
- } /* ncons not set */
-
-
- rTemp.ClobberGlobal();
- rTempAli.ClobberGlobal();
-
- return OK;
-
-} /*** end: readHMMWrapper() ***/
-
-
-
-
-
-
-/////////////////////////////////////////////////////////////////////////////
-/**
- * @brief Do precalculations for q and t to prepare comparison
- */
-void
-PrepareTemplate(HMM& q, HMM& t, int format)
-{
- if (format==0) // HHM format
- {
- // Add transition pseudocounts to template
- t.AddTransitionPseudocounts();
-
- // Generate an amino acid frequency matrix from t.f[i][a] with max pseudocounts (tau=1) -> t.g[i][a]
- t.PreparePseudocounts();
-
- // Add amino acid pseudocounts to template: t.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
- t.AddAminoAcidPseudocounts();
- }
- else // HHMER format
- {
- // Don't add transition pseudocounts to template
- // t.AddTransitionPseudocounts(par.gapd, par.gape, par.gapf, par.gapg, par.gaph, par.gapi, 0.0);
-
- // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
- t.PreparePseudocounts();
-
- // DON'T ADD amino acid pseudocounts to temlate: pcm=0! t.p[i][a] = t.f[i][a]
- t.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
- }
-
- // Modify transition probabilities to include SS-dependent penalties
- if (par.ssgap) t.UseSecStrucDependentGapPenalties();
-
- if (par.forward>=1) t.Log2LinTransitionProbs(1.0);
-
- // Factor Null model into HMM t
- // ATTENTION! t.p[i][a] is divided by pnul[a] (for reasons of efficiency) => do not reuse t.p
- t.IncludeNullModelInHMM(q,t); // Can go BEFORE the loop if not dependent on template
-
- return;
-}
-
-
-/*** end of hhfunc-C.h ***/
git://source.jalview.org / jabaws.git / blobdiff