clustal-omega-1.2.4/src/clustal-omega.c
/**
 * @brief Defines the alignment order by calculating a guide tree. In
 * a first-step pairwise distances will be calculated (or read from a
 * file). In a second step those distances will be clustered and a
 * guide-tree created. Steps 1 and 2 will be skipped if a guide-tree
 * file was given, in which case the guide-tree will be just read from
 * the file.
 *
 * @param[out] piOrderLR_p
 * order in which nodes/profiles are to be merged/aligned
 * @param[out] pdSeqWeights_p
 * Sequence weights
 * @param[out] pdSeqWeights_p
 * Sequence weights
 * @param[in] prMSeq
 * The sequences from which the alignment order is to be calculated
 * @param[in] iPairDistType
 * Method of pairwise distance comparison
 * @param[in] pcDistmatInfile
 * If not NULL distances will be read from this file instead of being
 * calculated
 * @param[in] pcDistmatOutfile
 * If not NULL computed pairwise distances will be written to this file
 * @param[in] iClusteringType
 * Clustering method to be used to cluster the pairwise distances
 * @param[in] pcGuidetreeInfile
 * If not NULL guidetree will be read from this file. Skips pairwise
 * distance and guidetree computation
 * @param[in] pcGuidetreeOutfile
 * If not NULL computed guidetree will be written to this file
 * @param[in] bUseMbed
 * If TRUE, fast mBed guidetree computation will be employed
 *
 * @return Non-zero on error
 *
 */
@brief Defines the alignment order by calculating a guide tree. In
a first-step pairwise distances will be calculated (or read from a
file). In a second step those distances will be clustered and a
guide-tree created. Steps 1 and 2 will be skipped if a guide-tree
file was given, in which case the guide-tree will be just read from
the file.



@return Non-zero on error
@brief Defines the alignment order by calculating a guide tree . In a first - step pairwise distances will be calculated ( or read from a file ) . In a second step those distances will be clustered and a guide - tree created . Steps 1 and 2 will be skipped if a guide - tree file was given in which case the guide - tree will be just read from the file . @return Non - zero on error
int
AlignmentOrder(int **piOrderLR_p, double **pdSeqWeights_p, mseq_t *prMSeq,
               int iPairDistType, char *pcDistmatInfile, char *pcDistmatOutfile,
               int iClusteringType, int iClustersizes, 
               char *pcGuidetreeInfile, char *pcGuidetreeOutfile, char *pcClusterFile, 
               bool bUseMbed, bool bPercID)
{
    symmatrix_t *distmat = NULL;
    tree_t *prTree = NULL;
    int i = 0;
    if (2==prMSeq->nseqs) {
        Log(&rLog, LOG_VERBOSE,
            "Have only two sequences: No need to compute pairwise score and compute a tree.");
        if (NULL != pcDistmatOutfile){
            Log(&rLog, LOG_WARN, "Have only two sequences: Will not calculate/print distance matrix.");
        }
        (*piOrderLR_p) = (int*) CKMALLOC(DIFF_NODE * 3 * sizeof(int));
        (*piOrderLR_p)[DIFF_NODE*0+LEFT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*0+RGHT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*0+PRNT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*1+LEFT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*1+RGHT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*1+PRNT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*2+LEFT_NODE] = 0;
        (*piOrderLR_p)[DIFF_NODE*2+RGHT_NODE] = 1;
        (*piOrderLR_p)[DIFF_NODE*2+PRNT_NODE] = 2;
#if USE_WEIGHTS
         (*pdWeights_p) = (double *) CKMALLOC(uNodeCount * sizeof(double));
         (*pdWeights_p)[0] = 0.5;
         (*pdWeights_p)[1] = 0.5;
#endif
        return OK;
    }
    if (NULL != pcGuidetreeInfile) {
        Log(&rLog, LOG_INFO, "Reading guide-tree from %s", pcGuidetreeInfile);
        if (GuideTreeFromFile(&prTree, prMSeq, pcGuidetreeInfile)) {
            Log(&rLog, LOG_ERROR, "Reading of guide tree %s failed.", pcGuidetreeInfile);
            return FAILURE;
        }
    } else {
        if (bUseMbed) {
            if (NULL != pcDistmatInfile) {
                Log(&rLog, LOG_ERROR, "Can't input distance matrix when in mbed mode.");
                return FAILURE;                
            }
            if (Mbed(&prTree, prMSeq, iPairDistType, pcGuidetreeOutfile, iClustersizes, pcClusterFile)) {
                Log(&rLog, LOG_ERROR, "mbed execution failed.");
                return FAILURE;
            }
            Log(&rLog, LOG_INFO, "Guide-tree computation (mBed) done.");
            if (NULL != pcDistmatOutfile) {
                Log(&rLog, LOG_INFO,
                    "Ignoring request to write distance matrix (am in mBed mode)");
            }
        } else {
            if (PairDistances(&distmat, prMSeq, iPairDistType, bPercID, 
                              0, prMSeq->nseqs, 0, prMSeq->nseqs,
                              pcDistmatInfile, pcDistmatOutfile)) {
                Log(&rLog, LOG_ERROR, "Couldn't compute pair distances");
                return FAILURE;
            }
            if (CLUSTERING_UPGMA == iClusteringType) {
                char **labels;
                labels = (char**) CKMALLOC(prMSeq->nseqs * sizeof(char*));
                for (i=0; i<prMSeq->nseqs; i++) {
                    labels[i] = prMSeq->sqinfo[i].name;
                }
                GuideTreeUpgma(&prTree, labels, distmat, pcGuidetreeOutfile);
                Log(&rLog, LOG_INFO, "Guide-tree computation done.");
                CKFREE(labels);
            } else {
                Log(&rLog, LOG_FATAL, "INTERNAL ERROR %s",
                      "clustering method should have been checked before");
            }
        } 
    } 
#if USE_WEIGHTS
    Log(&rLog, LOG_INFO, "Calculating sequence weights");
    CalcClustalWeights(pdSeqWeights_p, prTree);
    for (i = 0; i < GetLeafCount(prTree); i++) {
        Log(&rLog, LOG_VERBOSE,
            "Weight for seq no %d: %s = %f",
             i, prMSeq->sqinfo[i].name, (*pdSeqWeights_p)[i]);
    }
#else
    Log(&rLog, LOG_DEBUG, "Not using weights");
#endif
    TraverseTree(piOrderLR_p, prTree, prMSeq);
    if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
        uint uNodeIndex;
        FILE *fp = LogGetFP(&rLog, LOG_INFO);
        Log(&rLog, LOG_DEBUG, "left/right order after tree traversal");
        for (uNodeIndex = 0; uNodeIndex < GetNodeCount(prTree); uNodeIndex++) {
            fprintf(fp, "%3d:\t%2d/%2d -> %d\n", i,
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+LEFT_NODE],
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+RGHT_NODE],
                   (*piOrderLR_p)[DIFF_NODE*uNodeIndex+PRNT_NODE]);
        }
    }
    FreeMuscleTree(prTree);
    FreeSymMatrix(&distmat);
#if 0
    Log(&rLog, LOG_FATAL, "DEBUG EXIT before leaving %s", __FUNCTION__);
#endif
    return OK;
}
int AlignmentOrder ( int * * piOrderLR_p , double * * pdSeqWeights_p , mseq_t * prMSeq , int iPairDistType , char * pcDistmatInfile , char * pcDistmatOutfile , int iClusteringType , int iClustersizes , char * pcGuidetreeInfile , char * pcGuidetreeOutfile , char * pcClusterFile , bool bUseMbed , bool bPercID ) { symmatrix_t * distmat = NULL ; tree_t * prTree = NULL ; int i = 0 ; if ( 2 == prMSeq -> nseqs ) { Log ( & rLog , LOG_VERBOSE , " " ) ; if ( NULL != pcDistmatOutfile ) { Log ( & rLog , LOG_WARN , " " ) ; } ( * piOrderLR_p ) = ( int * ) CKMALLOC ( DIFF_NODE * 3 * sizeof ( int ) ) ; ( * piOrderLR_p ) [ DIFF_NODE * 0 + LEFT_NODE ] = 0 ; ( * piOrderLR_p ) [ DIFF_NODE * 0 + RGHT_NODE ] = 0 ; ( * piOrderLR_p ) [ DIFF_NODE * 0 + PRNT_NODE ] = 0 ; ( * piOrderLR_p ) [ DIFF_NODE * 1 + LEFT_NODE ] = 1 ; ( * piOrderLR_p ) [ DIFF_NODE * 1 + RGHT_NODE ] = 1 ; ( * piOrderLR_p ) [ DIFF_NODE * 1 + PRNT_NODE ] = 1 ; ( * piOrderLR_p ) [ DIFF_NODE * 2 + LEFT_NODE ] = 0 ; ( * piOrderLR_p ) [ DIFF_NODE * 2 + RGHT_NODE ] = 1 ; ( * piOrderLR_p ) [ DIFF_NODE * 2 + PRNT_NODE ] = 2 ; #if USE_WEIGHTS ( * pdWeights_p ) = ( double * ) CKMALLOC ( uNodeCount * sizeof ( double ) ) ; ( * pdWeights_p ) [ 0 ] = 0.5 ; ( * pdWeights_p ) [ 1 ] = 0.5 ; #endif return OK ; } if ( NULL != pcGuidetreeInfile ) { Log ( & rLog , LOG_INFO , " " , pcGuidetreeInfile ) ; if ( GuideTreeFromFile ( & prTree , prMSeq , pcGuidetreeInfile ) ) { Log ( & rLog , LOG_ERROR , " " , pcGuidetreeInfile ) ; return FAILURE ; } } else { if ( bUseMbed ) { if ( NULL != pcDistmatInfile ) { Log ( & rLog , LOG_ERROR , " " ) ; return FAILURE ; } if ( Mbed ( & prTree , prMSeq , iPairDistType , pcGuidetreeOutfile , iClustersizes , pcClusterFile ) ) { Log ( & rLog , LOG_ERROR , " " ) ; return FAILURE ; } Log ( & rLog , LOG_INFO , " " ) ; if ( NULL != pcDistmatOutfile ) { Log ( & rLog , LOG_INFO , " " ) ; } } else { if ( PairDistances ( & distmat , prMSeq , iPairDistType , bPercID , 0 , prMSeq -> nseqs , 0 , prMSeq -> nseqs , pcDistmatInfile , pcDistmatOutfile ) ) { Log ( & rLog , LOG_ERROR , " " ) ; return FAILURE ; } if ( CLUSTERING_UPGMA == iClusteringType ) { char * * labels ; labels = ( char * * ) CKMALLOC ( prMSeq -> nseqs * sizeof ( char * ) ) ; for ( i = 0 ; i < prMSeq -> nseqs ; i ++ ) { labels [ i ] = prMSeq -> sqinfo [ i ] . name ; } GuideTreeUpgma ( & prTree , labels , distmat , pcGuidetreeOutfile ) ; Log ( & rLog , LOG_INFO , " " ) ; CKFREE ( labels ) ; } else { Log ( & rLog , LOG_FATAL , " " , " " ) ; } } } #if USE_WEIGHTS Log ( & rLog , LOG_INFO , " " ) ; CalcClustalWeights ( pdSeqWeights_p , prTree ) ; for ( i = 0 ; i < GetLeafCount ( prTree ) ; i ++ ) { Log ( & rLog , LOG_VERBOSE , " " , i , prMSeq -> sqinfo [ i ] . name , ( * pdSeqWeights_p ) [ i ] ) ; } #else Log ( & rLog , LOG_DEBUG , " " ) ; #endif TraverseTree ( piOrderLR_p , prTree , prMSeq ) ; if ( rLog . iLogLevelEnabled <= LOG_DEBUG ) { uint uNodeIndex ; FILE * fp = LogGetFP ( & rLog , LOG_INFO ) ; Log ( & rLog , LOG_DEBUG , " " ) ; for ( uNodeIndex = 0 ; uNodeIndex < GetNodeCount ( prTree ) ; uNodeIndex ++ ) { fprintf ( fp , " \t \n " , i , ( * piOrderLR_p ) [ DIFF_NODE * uNodeIndex + LEFT_NODE ] , ( * piOrderLR_p ) [ DIFF_NODE * uNodeIndex + RGHT_NODE ] , ( * piOrderLR_p ) [ DIFF_NODE * uNodeIndex + PRNT_NODE ] ) ; } } FreeMuscleTree ( prTree ) ; FreeSymMatrix ( & distmat ) ; #if 0 Log ( & rLog , LOG_FATAL , " " , __FUNCTION__ ) ; #endif return OK ; }
@brief Defines the alignment order by calculating a guide tree.
@brief Defines the alignment order by calculating a guide tree .
/* pairwise distance matrix (tmat in 1.83) */ /* guide tree */ /* Shortcut for only two sequences: Do not compute k-tuple
     * distances. Use the same logic as in TraverseTree() to setup
     * piOrderLR_p. Changes there will have to be reflected here as
     * well. */ /* root */ /* Same logic as CalcClustalWeights(). Changes there will
           have to be reflected here as well. */ /* compute distance & guide tree, alternatively read distances or
     * guide tree from file
     *
     */ /* clustering of distances to get guide tree
             */ /* did not use mBed */ /* had to calculate tree (did not read from file) */ /* derive sequence weights from tree
     *
     */ /* define traversing order of tree
     *
     */ /* FIXME: debug only, FS */