+++ /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: pair_dist.c 242 2011-05-27 14:04:21Z andreas $
- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <stdlib.h>
-#include <ctype.h>
-#include <assert.h>
-#include <time.h>
-
-/* only neededfor iNumberOfThreads */
-#include "clustal-omega.h"
-
-#include "ktuple_pair.h"
-#include "pair_dist.h"
-#include "progress.h"
-#include "util.h"
-
-
-/* Up to rev 173 we had a USE_SYM_KTUPLE switch implemented here. When active
- * ktuple distances were computed twice for each pair and averaged. Idea was
- * to avoid assymmetries in the pairwise scores (score(a, b) is often not the
- * same as score(b, a)). Results on BAliBASE indicate that this is overkill:
- *
- * r92_default core columns: avg-sp=0.800656 avg-tc=0.47711 (of total 218)
- * r93-mod--norm-ktuple/ core columns: avg-sp=0.800656 avg-tc=0.47711 (of total 218)
- * r93-mod--sym-ktuple/ core columns: avg-sp=0.801083 avg-tc=0.476544 (of total 217)
- * r93-mod--rand-ktuple-1 core columns: avg-sp=0.799289 avg-tc=0.468028 (of total 218)
- * r93-mod--rand-ktuple-2 core columns: avg-sp=0.801654 avg-tc=0.47659 (of total 217)
- * r93-mod--rand-ktuple-3 core columns: avg-sp=0.800234 avg-tc=0.474908 (of total 218)
- * r93-mod--rand-ktuple-4 core columns: avg-sp=0.800573 avg-tc=0.476514 (of total 218)
- * r93-mod--rand-ktuple-5 core columns: avg-sp=0.799679 avg-tc=0.468716 (of total 218)
- *
- */
-
-static double
-KimuraCorrection(double frac_id);
-
-static int
-SquidIdPairDist(symmatrix_t *tmat, mseq_t *mseq,
- int istart, int iend,
- int jstart, int jend,
- bool use_KimuraCorrection, progress_t *prProgress,
- unsigned long int *ulStepNo, unsigned long int ulTotalStepNo);
-
-/* Taken from Muscle's msadistkimura.cpp */
-static int DAYHOFF_PAMS[]={
- 195, /* 75.0% observed d; 195 PAMs estimated = 195% estimated d */
- 196, /* 75.1% observed d; 196 PAMs estimated */
- 197, 198, 199, 200, 200, 201, 202, 203,
- 204, 205, 206, 207, 208, 209, 209, 210, 211, 212,
- 213, 214, 215, 216, 217, 218, 219, 220, 221, 222,
- 223, 224, 226, 227, 228, 229, 230, 231, 232, 233,
- 234, 236, 237, 238, 239, 240, 241, 243, 244, 245,
- 246, 248, 249, 250, /* 250 PAMs = 80.3% observed d */
- 252, 253, 254, 255, 257, 258,
- 260, 261, 262, 264, 265, 267, 268, 270, 271, 273,
- 274, 276, 277, 279, 281, 282, 284, 285, 287, 289,
- 291, 292, 294, 296, 298, 299, 301, 303, 305, 307,
- 309, 311, 313, 315, 317, 319, 321, 323, 325, 328,
- 330, 332, 335, 337, 339, 342, 344, 347, 349, 352,
- 354, 357, 360, 362, 365, 368, 371, 374, 377, 380,
- 383, 386, 389, 393, 396, 399, 403, 407, 410, 414,
- 418, 422, 426, 430, 434, 438, 442, 447, 451, 456,
- 461, 466, 471, 476, 482, 487, 493, 498, 504, 511,
- 517, 524, 531, 538, 545, 553, 560, 569, 577, 586,
- 595, 605, 615, 626, 637, 649, 661, 675, 688, 703,
- 719, 736, 754, 775, 796, 819, 845, 874, 907, 945,
- /* 92.9% observed; 945 PAMs */
- 988 /* 93.0% observed; 988 PAMs */
-};
-static int DAYHOFF_TABLE_ENTRIES = sizeof(DAYHOFF_PAMS)/sizeof(DAYHOFF_PAMS[0]);
-
-
-
-/**
- *
- * @brief Compute Kimura corrected distance.
- *
- * Original Muscle documentation following:
- * """
- * This is defined to be:
- * log_e(1 - p - p*p/5)
- * where p is the fraction of residues that differ, i.e.:
- * p = (1 - fractional_conservation)
- * This measure is infinite for p = 0.8541 and is considered
- * unreliable for p >= 0.75 (according to the ClustalW docs).
- * ClustalW uses a table lookup for values > 0.75. The following table
- * was copied from the ClustalW file dayhoff.h.
- * """
- *
- * @note copied from Muscle's msadistkimura.cpp:KimuraDist()
- *
- * @warning For protein only (uses Dayhoff substitution parameters)
- *
- * @param[in] p
- * distance, e.g. 1.0 - fractional/relative identity
- *
- * @return The Kimura corrected distance
- *
- */
-double
-KimuraCorrection(double p)
-{
- int table_index;
-
- /* Typical case: use Kimura's empirical formula */
- if (p < 0.75)
- return -log(1 - p - (p*p)/5);
-
- /* Per ClustalW, return 10.0 for anything over 93% */
- if (p > 0.93)
- return 10.0;
-
- /* If 0.75 >= p <= 0.93, use table lookup */
- table_index = (int) ((p - 0.75)*1000 + 0.5);
- if (table_index < 0 || table_index >= DAYHOFF_TABLE_ENTRIES)
- Log(&rLog, LOG_FATAL, "Internal error in %s:%s", __FILE__, __FUNCTION__);
-
- return DAYHOFF_PAMS[table_index] / 100.0;
-}
-/*** end: KimuraCorrection() ***/
-
-
-
-
-/**
- * @brief Compute distances between all aligned sequence pairs using
- * squid's PairwiseIdentity, which is: idents / MIN(len1, len2)
- *
- * @param[out] tmat
- * Where to store the computed distances
- * @param[in] mseq
- * The aligned sequences
- * @param[in] istart
- * For distances [i][j] i>=istart, i<j
- * @param[in] iend
- * For distances [i][j] i<iend, i<j
- * @param[in] jstart
- * For distances [i][j] j>=jstart, i<j
- * @param[in] jend
- * For distances [i][j] i<j<jend, i<j
- * @param[in] use_kimura
- * Use Kimura corrected values (Proteins only)
- *
- * @return Non-zero on error
- *
- */
-int
-SquidIdPairDist(symmatrix_t *tmat, mseq_t *mseq,
- int istart, int iend,
- int jstart, int jend,
- bool use_kimura, progress_t *prProgress,
- unsigned long int *ulStepNo, unsigned long int ulTotalStepNo)
-{
- int i, j; /* aux */
- /* progress_t *prProgress; */
- bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE;
- /* unsigned long int ulStepNo;
- unsigned long ulTotalStepNo; */
-
- assert(NULL != tmat);
- assert(NULL != mseq);
-
- if (TRUE != mseq->aligned) {
- Log(&rLog, LOG_ERROR, "Sequences need to be aligned (%s)", __FUNCTION__);
- return -1;
- }
- if (SEQTYPE_PROTEIN != mseq->seqtype && TRUE == use_kimura) {
- Log(&rLog, LOG_WARN, "Using Kimura distance corretion which includes Dayhoff substitution table lookup for non-protein sequences");
- }
-
- NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
- "Pairwise distance calculation progress", bPrintCR);
- /* estimation of total number of steps (if istart and jstart are
- * both 0)
- */
- /* ulTotalStepNo = iend*jend - iend*iend/2 + iend/2;
- ulStepNo = 0; */
- /*LOG_DEBUG("istart=%d iend=%d jstart=%d jend=%d", istart, iend, jstart, jend);*/
- for (i=istart; i<iend; ++i) {
- /* by definition a sequence compared to itself should give a
- score of 0 */
- SymMatrixSetValue(tmat, i, i, 0.0);
-#ifdef HAVE_OPENMP
- #pragma omp critical(squidid)
-#endif
- {
- ProgressLog(prProgress, *ulStepNo, ulTotalStepNo, FALSE);
- }
-
- for (j=MAX(i+1, jstart); j<jend; ++j) {
- float dist;
- dist = 1.0 - PairwiseIdentity(mseq->seq[i], mseq->seq[j]);
-
-
-#ifdef HAVE_OPENMP
- #pragma omp atomic
-#endif
- (*ulStepNo)++;
- /*LOG_DEBUG("%d:%d raw dist = %f", i, j, dist);*/
- if (use_kimura) {
- dist = KimuraCorrection(dist);
- /*LOG_DEBUG("cor dist = %f", dist);*/
- }
- SymMatrixSetValue(tmat, i, j, dist);
-#ifdef HAVE_OPENMP
- #pragma omp critical(squidid)
-#endif
- {
- Log(&rLog, LOG_DEBUG, "Aligned distance for sequence pair %d:%d= %lg",
- i+1, j+1, dist);
- }
- }
- }
-
- return 0;
-}
-/*** end: SquidIdPairDist() ***/
-
-
-
-
-/**
- * @brief compute or read precomputed distances for given sequences
- *
- * @param[out] distmat
- * Distances will be written to this matrix. will be allocated here as
- * well. Caller must free with FreeSymMatrix()
- * @param[in] mseq
- * Distances will be computed for these sequences
- * @param[in] pairdist_type
- * Type of pairwise distance comparison
- * @param[in] fdist_in
- * If not NULL, sequences will be written from this file instead of
- * computing them
- * @param[in] istart
- * Compute distances for sequences i:j, i>=istart, i<j.
- * Usually 0.
- * @param[in] iend
- * Compute distances for sequences i:j, i<iend, i<j
- * Usually mseq->nseqs.
- * @param[in] jstart
- * Compute distances for sequences i:j, j>=jstart, i<j
- * Usually 0.
- * @param[in] jend
- * Compute distances for sequences i:j, j<iend, i<j
- * Usually mseq->nseqs.
- * @param[in] fdist_out
- * If not NULL, distances will be written to this files
- *
- *
- */
-int
-PairDistances(symmatrix_t **distmat, mseq_t *mseq, int pairdist_type,
- int istart, int iend,
- int jstart, int jend,
- char *fdist_in, char *fdist_out)
-{
- int uSeqIndex;
- unsigned long int ulStepNo = 0, ulTotalStepNo; /* DD: moved from SquidIdPairDist so progress bar works multithreaded */
- int iChunk, iChunkStart, iChunkEnd;
- int iChunkStarts[iNumberOfThreads];
- int iChunkEnds[iNumberOfThreads];
- progress_t *prProgress = NULL;
- int iSquidSuccess = 0;
- bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE;
-
- assert(NULL!=distmat);
- assert(NULL!=mseq);
- assert(istart<iend);
- assert(jstart<jend);
-
- if (NewSymMatrix(distmat, iend, jend)!=0) {
- Log(&rLog, LOG_FATAL, "%s", "Memory allocation for distance matrix failed");
- }
-
-
- /* compute pairwise distances or read from file
- *
- */
-#if 0
-#include "random-dist.h"
-#else
- if (NULL != fdist_in) {
- Log(&rLog, LOG_FATAL, "FIXME: reading of distance matrix from file not implemented");
-
- } else {
-
- /* break into chunks, one for each thread
- matrix is a triangle, not a square
- hence making even chunk sizes is slightly fiddlier
- */
- ulTotalStepNo = iend*jend - iend*iend/2 + iend/2;
-
- /* FIXME: can get rid of iChunkStart, iChunkEnd now that we're using the arrays */
- iChunkStart = iend;
- for(iChunk = 0; iChunk <= iNumberOfThreads; iChunk++)
- {
- iChunkEnd = iChunkStart;
- if(iChunk == iNumberOfThreads - 1)
- iChunkStart = 0;
- else
- iChunkStart = iend - ((double)(iend - istart) * sqrt(((double)iChunk + 1.0)/(double)iNumberOfThreads));
- iChunkStarts[iChunk] = iChunkStart;
- iChunkEnds[iChunk] = iChunkEnd;
- }
-
- if (PAIRDIST_KTUPLE == pairdist_type) {
-
- Log(&rLog, LOG_INFO, "Calculating pairwise ktuple-distances...");
-
- NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
- "Ktuple-distance calculation progress", bPrintCR);
-#ifdef HAVE_OPENMP
- #pragma omp parallel for private(iChunk) schedule(dynamic)
-#endif
- for(iChunk = 0; iChunk < iNumberOfThreads; iChunk++)
- {
- KTuplePairDist((*distmat), mseq, iChunkStarts[iChunk],
- iChunkEnds[iChunk], jstart, jend, NULL, prProgress,
- &ulStepNo, ulTotalStepNo);
- }
-
-#if 0
- printf("total ops %d\n", ulStepNo);
-#endif
- /* old format:
- KTuplePairDist((*distmat), mseq,
- istart, iend,
- jstart, jend, NULL); */
-
- } else if (PAIRDIST_SQUIDID == pairdist_type) {
- Log(&rLog, LOG_INFO, "Calculating pairwise aligned identity distances...");
-
- NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
- "Pairwise identity calculation progress", bPrintCR);
-#ifdef HAVE_OPENMP
- #pragma omp parallel for private(iChunk) schedule(dynamic)
-#endif
- for(iChunk = 0; iChunk < iNumberOfThreads; iChunk++)
- {
- iSquidSuccess = SquidIdPairDist((*distmat), mseq,
- iChunkStarts[iChunk], iChunkEnds[iChunk],
- jstart, jend, FALSE, prProgress,
- &ulStepNo, ulTotalStepNo);
- }
- if(iSquidSuccess != 0)
- return -1;
-
- } else if (PAIRDIST_SQUIDID_KIMURA == pairdist_type) {
- Log(&rLog, LOG_INFO, "Calculating Kimura-corrected pairwise aligned identity distances...");
- NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
- "Pairwise identity calculation progress", bPrintCR);
-#ifdef HAVE_OPENMP
- #pragma omp parallel for private(iChunk) schedule(dynamic)
-#endif
- for(iChunk = 0; iChunk < iNumberOfThreads; iChunk++)
- {
- iSquidSuccess = SquidIdPairDist((*distmat), mseq,
- iChunkStarts[iChunk], iChunkEnds[iChunk],
- jstart, jend, TRUE, prProgress,
- &ulStepNo, ulTotalStepNo);
- }
- if(iSquidSuccess != 0)
- return -1;
- } else {
- Log(&rLog, LOG_FATAL, "INTERNAL ERROR: don't know about pairdist_type %d",
- pairdist_type);
- }
- }
-#endif /* random/proper distance calculation */
-
-
- /* optional printing of matrix to file
- */
- if (NULL != fdist_out) {
- /* need a copy of sequence names for printing */
- char **names;
- names = (char **)CKMALLOC(mseq->nseqs * sizeof(char*));
- for (uSeqIndex=0; uSeqIndex<mseq->nseqs; uSeqIndex++) {
- names[uSeqIndex] = mseq->sqinfo[uSeqIndex].name;
- }
-
- SymMatrixPrint((*distmat), names, fdist_out);
-
- Log(&rLog, LOG_INFO, "Pairwise distance matrix written to %s",
- fdist_out);
- CKFREE(names);
- }
-
-#if 0
-#include "distance-distrib.h"
-#endif
-
- if (NULL != prProgress) {
- ProgressDone(prProgress);
- FreeProgress(&prProgress);
- }
-
- return 0;
-}
-/*** end: PairDistances() ***/
-
-
-
git://source.jalview.org / jabaws.git / blobdiff