开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
block_align.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：38k
源码类别：

生物技术

开发平台：

C/C++

block_align.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: block_align.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:11:09  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: block_align.cpp,v 1000.2 2004/06/01 18:11:09 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors:  Paul Thiessen
  36. *
  37. * File Description:
  38. *      Dynamic programming-based alignment algorithms: block aligner
  39. *
  40. * ===========================================================================
  41. */
  43. #include <ncbi_pch.hpp>
  44. #include <corelib/ncbistd.hpp>
  45. #include <corelib/ncbidiag.hpp>
  46. #include <corelib/ncbi_limits.h>
  48. #include <vector>
  49. #include <list>
  50. #include <algorithm>
  52. #include <algo/structure/struct_dp/struct_dp.h>
  54. USING_NCBI_SCOPE;
  57. const int DP_NEGATIVE_INFINITY = kMin_Int;
  58. const unsigned int DP_POSITIVE_INFINITY = kMax_UInt;
  59. const unsigned int DP_UNFROZEN_BLOCK = kMax_UInt;
  61. BEGIN_SCOPE(struct_dp)
  63. #define ERROR_MESSAGE(s) ERR_POST(Error << "block_align: " << s << '!')
  64. #define WARNING_MESSAGE(s) ERR_POST(Warning << "block_align: " << s)
  65. #define INFO_MESSAGE(s) ERR_POST(Info << "block_align: " << s)
  67. #define NO_TRACEBACK kMax_UInt
  69. class Cell
  70. {
  71. public:
  72.     int score;
  73.     unsigned int tracebackResidue;
  74.     Cell(void) : score(DP_NEGATIVE_INFINITY), tracebackResidue(NO_TRACEBACK) { }
  75. };
  77. class Matrix
  78. {
  79. public:
  80.     typedef vector < Cell > ResidueRow;
  81.     typedef vector < ResidueRow > Grid;
  82.     Grid grid;
  83.     Matrix(unsigned int nBlocks, unsigned int nResidues) : grid(nBlocks + 1)
  84.         { for (unsigned int i=0; i<nBlocks; ++i) grid[i].resize(nResidues + 1); }
  85.     ResidueRow& operator [] (unsigned int block) { return grid[block]; }
  86.     const ResidueRow& operator [] (unsigned int block) const { return grid[block]; }
  87. };
  89. // checks to make sure frozen block positions are legal (borrowing my own code from blockalign.c)
  90. int ValidateFrozenBlockPositions(const DP_BlockInfo *blocks,
  91.     unsigned int queryFrom, unsigned int queryTo, bool checkGapSum)
  92. {
  93.     static const unsigned int NONE = kMax_UInt;
  94.     unsigned int
  95.         block,
  96.         unfrozenBlocksLength = 0,
  97.         prevFrozenBlock = NONE,
  98.         maxGapsLength = 0;
  100.     for (block=0; block<blocks->nBlocks; ++block) {
  102.         /* keep track of max gap space between frozen blocks */
  103.         if (block > 0)
  104.             maxGapsLength += blocks->maxLoops[block - 1];
  106.         /* to allow room for unfrozen blocks between frozen ones */
  107.         if (blocks->freezeBlocks[block] == DP_UNFROZEN_BLOCK) {
  108.             unfrozenBlocksLength += blocks->blockSizes[block];
  109.             continue;
  110.         }
  112.         /* check absolute block end positions */
  113.         if (blocks->freezeBlocks[block] < queryFrom) {
  114.             ERROR_MESSAGE("Frozen block " << (block+1) << " can't start < " << (queryFrom+1));
  115.             return STRUCT_DP_PARAMETER_ERROR;
  116.         }
  117.         if (blocks->freezeBlocks[block] + blocks->blockSizes[block] - 1 > queryTo) {
  118.             ERROR_MESSAGE("Frozen block " << (block+1) << " can't end > " << (queryTo+1));
  119.             return STRUCT_DP_PARAMETER_ERROR;
  120.         }
  122.         /* checks for legal distances between frozen blocks */
  123.         if (prevFrozenBlock != NONE) {
  124.             unsigned int prevFrozenBlockEnd =
  125.                 blocks->freezeBlocks[prevFrozenBlock] + blocks->blockSizes[prevFrozenBlock] - 1;
  127.             /* check for adequate room for unfrozen blocks between frozen blocks */
  128.             if (blocks->freezeBlocks[block] <= (prevFrozenBlockEnd + unfrozenBlocksLength)) {
  129.                 ERROR_MESSAGE("Frozen block " << (block+1) << " starts before end of prior frozen block, "
  130.                     "or doesn't leave room for intervening unfrozen block(s)");
  131.                 return STRUCT_DP_PARAMETER_ERROR;
  132.             }
  134.             /* check for too much gap space since last frozen block,
  135.                but if frozen blocks are adjacent, issue warning only */
  136.             if (checkGapSum &&
  137.                 blocks->freezeBlocks[block] > prevFrozenBlockEnd + 1 + unfrozenBlocksLength + maxGapsLength)
  138.             {
  139.                 if (prevFrozenBlock == block - 1) {
  140.                     WARNING_MESSAGE("Frozen block " << (block+1) << " is further than allowed loop length"
  141.                         " from adjacent frozen block " << (prevFrozenBlock+1));
  142.                 } else {
  143.                     ERROR_MESSAGE("Frozen block " << (block+1) << " is too far away from prior frozen block"
  144.                         " given allowed intervening loop length(s) " << maxGapsLength
  145.                         << " plus unfrozen block length(s) " << unfrozenBlocksLength);
  146.                     return STRUCT_DP_PARAMETER_ERROR;
  147.                 }
  148.             }
  149.         }
  151.         /* reset counters after each frozen block */
  152.         prevFrozenBlock = block;
  153.         unfrozenBlocksLength = maxGapsLength = 0;
  154.     }
  156.     return STRUCT_DP_OKAY;
  157. }
  159. // fill matrix values; return true on success. Matrix must have only default values when passed in.
  160. int CalculateGlobalMatrix(Matrix& matrix,
  161.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
  162.     unsigned int queryFrom, unsigned int queryTo)
  163. {
  164.     unsigned int block, residue, prevResidue, lastBlock = blocks->nBlocks - 1;
  165.     int score = 0, sum;
  167.     // find possible block positions, based purely on block lengths
  168.     vector < unsigned int > firstPos(blocks->nBlocks), lastPos(blocks->nBlocks);
  169.     for (block=0; block<=lastBlock; ++block) {
  170.         if (block == 0) {
  171.             firstPos[0] = queryFrom;
  172.             lastPos[lastBlock] = queryTo - blocks->blockSizes[lastBlock] + 1;
  173.         } else {
  174.             firstPos[block] = firstPos[block - 1] + blocks->blockSizes[block - 1];
  175.             lastPos[lastBlock - block] =
  176.                 lastPos[lastBlock - block + 1] - blocks->blockSizes[lastBlock - block];
  177.         }
  178.     }
  180.     // further restrict the search if blocks are frozen
  181.     for (block=0; block<=lastBlock; ++block) {
  182.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  183.             if (blocks->freezeBlocks[block] < firstPos[block] ||
  184.                 blocks->freezeBlocks[block] > lastPos[block])
  185.             {
  186.                 ERROR_MESSAGE("CalculateGlobalMatrix() - frozen block "
  187.                     << (block+1) << " does not leave room for unfrozen blocks");
  188.                 return STRUCT_DP_PARAMETER_ERROR;
  189.             }
  190.             firstPos[block] = lastPos[block] = blocks->freezeBlocks[block];
  191.         }
  192.     }
  194.     // fill in first row with scores of first block at all possible positions
  195.     for (residue=firstPos[0]; residue<=lastPos[0]; ++residue)
  196.         matrix[0][residue - queryFrom].score = BlockScore(0, residue);
  198.     // for each successive block, find the best allowed pairing of the block with the previous block
  199.     bool blockScoreCalculated;
  200.     for (block=1; block<=lastBlock; ++block) {
  201.         for (residue=firstPos[block]; residue<=lastPos[block]; ++residue) {
  202.             blockScoreCalculated = false;
  204.             for (prevResidue=firstPos[block - 1]; prevResidue<=lastPos[block - 1]; ++prevResidue) {
  206.                 // current block must come after the previous block
  207.                 if (residue < prevResidue + blocks->blockSizes[block - 1])
  208.                     break;
  210.                 // cut off at max loop length from previous block, but not if both blocks are frozen
  211.                 if (residue > prevResidue + blocks->blockSizes[block - 1] + blocks->maxLoops[block - 1] &&
  212.                         (blocks->freezeBlocks[block] == DP_UNFROZEN_BLOCK ||
  213.                          blocks->freezeBlocks[block - 1] == DP_UNFROZEN_BLOCK))
  214.                     continue;
  216.                 // make sure previous block is at an allowed position
  217.                 if (matrix[block - 1][prevResidue - queryFrom].score == DP_NEGATIVE_INFINITY)
  218.                     continue;
  220.                 // get score at this position
  221.                 if (!blockScoreCalculated) {
  222.                     score = BlockScore(block, residue);
  223.                     if (score == DP_NEGATIVE_INFINITY)
  224.                         break;
  225.                     blockScoreCalculated = true;
  226.                 }
  228.                 // find highest sum of scores for allowed pairing of this block with any previous
  229.                 sum = score + matrix[block - 1][prevResidue - queryFrom].score;
  230.                 if (sum > matrix[block][residue - queryFrom].score) {
  231.                     matrix[block][residue - queryFrom].score = sum;
  232.                     matrix[block][residue - queryFrom].tracebackResidue = prevResidue;
  233.                 }
  234.             }
  235.         }
  236.     }
  238.     return STRUCT_DP_OKAY;
  239. }
  241. // fill matrix values; return true on success. Matrix must have only default values when passed in.
  242. int CalculateGlobalMatrixGeneric(Matrix& matrix,
  243.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
  244.     unsigned int queryFrom, unsigned int queryTo)
  245. {
  246.     unsigned int block, residue, prevResidue, lastBlock = blocks->nBlocks - 1;
  247.     int blockScore = 0, sum;
  248.     unsigned int loopPenalty;
  250.     // find possible block positions, based purely on block lengths
  251.     vector < unsigned int > firstPos(blocks->nBlocks), lastPos(blocks->nBlocks);
  252.     for (block=0; block<=lastBlock; ++block) {
  253.         if (block == 0) {
  254.             firstPos[0] = queryFrom;
  255.             lastPos[lastBlock] = queryTo - blocks->blockSizes[lastBlock] + 1;
  256.         } else {
  257.             firstPos[block] = firstPos[block - 1] + blocks->blockSizes[block - 1];
  258.             lastPos[lastBlock - block] =
  259.                 lastPos[lastBlock - block + 1] - blocks->blockSizes[lastBlock - block];
  260.         }
  261.     }
  263.     // further restrict the search if blocks are frozen
  264.     for (block=0; block<=lastBlock; ++block) {
  265.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  266.             if (blocks->freezeBlocks[block] < firstPos[block] ||
  267.                 blocks->freezeBlocks[block] > lastPos[block])
  268.             {
  269.                 ERROR_MESSAGE("CalculateGlobalMatrix() - frozen block "
  270.                     << (block+1) << " does not leave room for unfrozen blocks");
  271.                 return STRUCT_DP_PARAMETER_ERROR;
  272.             }
  273.             firstPos[block] = lastPos[block] = blocks->freezeBlocks[block];
  274.         }
  275.     }
  277.     // fill in first row with scores of first block at all possible positions
  278.     for (residue=firstPos[0]; residue<=lastPos[0]; ++residue)
  279.         matrix[0][residue - queryFrom].score = BlockScore(0, residue);
  281.     // for each successive block, find the best allowed pairing of the block with the previous block
  282.     bool blockScoreCalculated;
  283.     for (block=1; block<=lastBlock; ++block) {
  284.         for (residue=firstPos[block]; residue<=lastPos[block]; ++residue) {
  285.             blockScoreCalculated = false;
  287.             for (prevResidue=firstPos[block - 1]; prevResidue<=lastPos[block - 1]; ++prevResidue) {
  289.                 // current block must come after the previous block
  290.                 if (residue < prevResidue + blocks->blockSizes[block - 1])
  291.                     break;
  293.                 // make sure previous block is at an allowed position
  294.                 if (matrix[block - 1][prevResidue - queryFrom].score == DP_NEGATIVE_INFINITY)
  295.                     continue;
  297.                 // get loop score at this position; assume loop score zero if both frozen
  298.                 if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK &&
  299.                     blocks->freezeBlocks[block - 1] != DP_UNFROZEN_BLOCK)
  300.                 {
  301.                     loopPenalty = 0;
  302.                 } else {
  303.                     loopPenalty = LoopScore(block - 1, residue - prevResidue - blocks->blockSizes[block - 1]);
  304.                     if (loopPenalty == DP_POSITIVE_INFINITY)
  305.                         continue;
  306.                 }
  308.                 // get score at this position
  309.                 if (!blockScoreCalculated) {
  310.                     blockScore = BlockScore(block, residue);
  311.                     if (blockScore == DP_NEGATIVE_INFINITY)
  312.                         break;
  313.                     blockScoreCalculated = true;
  314.                 }
  316.                 // find highest sum of scores + loop score for allowed pairing of this block with previous
  317.                 sum = blockScore + matrix[block - 1][prevResidue - queryFrom].score - loopPenalty;
  318.                 if (sum > matrix[block][residue - queryFrom].score) {
  319.                     matrix[block][residue - queryFrom].score = sum;
  320.                     matrix[block][residue - queryFrom].tracebackResidue = prevResidue;
  321.                 }
  322.             }
  323.         }
  324.     }
  326.     return STRUCT_DP_OKAY;
  327. }
  329. // fill matrix values; return true on success. Matrix must have only default values when passed in.
  330. int CalculateLocalMatrix(Matrix& matrix,
  331.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
  332.     unsigned int queryFrom, unsigned int queryTo)
  333. {
  334.     unsigned int block, residue, prevResidue, lastBlock = blocks->nBlocks - 1, tracebackResidue = 0;
  335.     int score, sum, bestPrevScore;
  337.     // find last possible block positions, based purely on block lengths
  338.     vector < unsigned int > lastPos(blocks->nBlocks);
  339.     for (block=0; block<=lastBlock; ++block) {
  340.         if (blocks->blockSizes[block] > queryTo - queryFrom + 1) {
  341.             ERROR_MESSAGE("Block " << (block+1) << " too large for this query range");
  342.             return STRUCT_DP_PARAMETER_ERROR;
  343.         }
  344.         lastPos[block] = queryTo - blocks->blockSizes[block] + 1;
  345.     }
  347.     // first row: positive scores of first block at all possible positions
  348.     for (residue=queryFrom; residue<=lastPos[0]; ++residue) {
  349.         score = BlockScore(0, residue);
  350.         matrix[0][residue - queryFrom].score = (score > 0) ? score : 0;
  351.     }
  353.     // first column: positive scores of all blocks at first positions
  354.     for (block=1; block<=lastBlock; ++block) {
  355.         score = BlockScore(block, queryFrom);
  356.         matrix[block][0].score = (score > 0) ? score : 0;
  357.     }
  359.     // for each successive block, find the best positive scoring with a previous block, if any
  360.     for (block=1; block<=lastBlock; ++block) {
  361.         for (residue=queryFrom+1; residue<=lastPos[block]; ++residue) {
  363.             // get score at this position
  364.             score = BlockScore(block, residue);
  365.             if (score == DP_NEGATIVE_INFINITY)
  366.                 continue;
  368.             // find max score of any allowed previous block
  369.             bestPrevScore = DP_NEGATIVE_INFINITY;
  370.             for (prevResidue=queryFrom; prevResidue<=lastPos[block - 1]; ++prevResidue) {
  372.                 // current block must come after the previous block
  373.                 if (residue < prevResidue + blocks->blockSizes[block - 1])
  374.                     break;
  376.                 // cut off at max loop length from previous block
  377.                 if (residue > prevResidue + blocks->blockSizes[block - 1] + blocks->maxLoops[block - 1])
  378.                     continue;
  380.                 // keep maximum score
  381.                 if (matrix[block - 1][prevResidue - queryFrom].score > bestPrevScore) {
  382.                     bestPrevScore = matrix[block - 1][prevResidue - queryFrom].score;
  383.                     tracebackResidue = prevResidue;
  384.                 }
  385.             }
  387.             // extend alignment if the sum of this block's + previous block's score is positive
  388.             if (bestPrevScore > 0 && (sum=bestPrevScore+score) > 0) {
  389.                 matrix[block][residue - queryFrom].score = sum;
  390.                 matrix[block][residue - queryFrom].tracebackResidue = tracebackResidue;
  391.             }
  393.             // otherwise, start new alignment if score is positive
  394.             else if (score > 0)
  395.                 matrix[block][residue - queryFrom].score = score;
  396.         }
  397.     }
  399.     return STRUCT_DP_OKAY;
  400. }
  402. // fill matrix values; return true on success. Matrix must have only default values when passed in.
  403. int CalculateLocalMatrixGeneric(Matrix& matrix,
  404.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
  405.     unsigned int queryFrom, unsigned int queryTo)
  406. {
  407.     unsigned int block, residue, prevResidue, loopPenalty,
  408.         lastBlock = blocks->nBlocks - 1, tracebackResidue = 0;
  409.     int score, sum, bestPrevScore;
  411.     // find last possible block positions, based purely on block lengths
  412.     vector < unsigned int > lastPos(blocks->nBlocks);
  413.     for (block=0; block<=lastBlock; ++block) {
  414.         if (blocks->blockSizes[block] > queryTo - queryFrom + 1) {
  415.             ERROR_MESSAGE("Block " << (block+1) << " too large for this query range");
  416.             return STRUCT_DP_PARAMETER_ERROR;
  417.         }
  418.         lastPos[block] = queryTo - blocks->blockSizes[block] + 1;
  419.     }
  421.     // first row: positive scores of first block at all possible positions
  422.     for (residue=queryFrom; residue<=lastPos[0]; ++residue) {
  423.         score = BlockScore(0, residue);
  424.         matrix[0][residue - queryFrom].score = (score > 0) ? score : 0;
  425.     }
  427.     // first column: positive scores of all blocks at first positions
  428.     for (block=1; block<=lastBlock; ++block) {
  429.         score = BlockScore(block, queryFrom);
  430.         matrix[block][0].score = (score > 0) ? score : 0;
  431.     }
  433.     // for each successive block, find the best positive scoring with a previous block, if any
  434.     for (block=1; block<=lastBlock; ++block) {
  435.         for (residue=queryFrom+1; residue<=lastPos[block]; ++residue) {
  437.             // get score at this position
  438.             score = BlockScore(block, residue);
  439.             if (score == DP_NEGATIVE_INFINITY)
  440.                 continue;
  442.             // find max score of any allowed previous block
  443.             bestPrevScore = DP_NEGATIVE_INFINITY;
  444.             for (prevResidue=queryFrom; prevResidue<=lastPos[block - 1]; ++prevResidue) {
  446.                 // current block must come after the previous block
  447.                 if (residue < prevResidue + blocks->blockSizes[block - 1])
  448.                     break;
  450.                 // make sure previous block is at an allowed position
  451.                 if (matrix[block - 1][prevResidue - queryFrom].score == DP_NEGATIVE_INFINITY)
  452.                     continue;
  454.                 // get loop score
  455.                 loopPenalty = LoopScore(block - 1, residue - prevResidue - blocks->blockSizes[block - 1]);
  456.                 if (loopPenalty == DP_POSITIVE_INFINITY)
  457.                     continue;
  459.                 // keep maximum score
  460.                 sum = matrix[block - 1][prevResidue - queryFrom].score - loopPenalty;
  461.                 if (sum > bestPrevScore) {
  462.                     bestPrevScore = sum;
  463.                     tracebackResidue = prevResidue;
  464.                 }
  465.             }
  467.             // extend alignment if the sum of this block's + previous block's score is positive
  468.             if (bestPrevScore > 0 && (sum=bestPrevScore+score) > 0) {
  469.                 matrix[block][residue - queryFrom].score = sum;
  470.                 matrix[block][residue - queryFrom].tracebackResidue = tracebackResidue;
  471.             }
  473.             // otherwise, start new alignment if score is positive
  474.             else if (score > 0)
  475.                 matrix[block][residue - queryFrom].score = score;
  476.         }
  477.     }
  479.     return STRUCT_DP_OKAY;
  480. }
  482. int TracebackAlignment(const Matrix& matrix, unsigned int lastBlock, unsigned int lastBlockPos,
  483.     unsigned int queryFrom, DP_AlignmentResult *alignment)
  484. {
  485.     // trace backwards from last block/pos
  486.     vector < unsigned int > blockPositions;
  487.     unsigned int block = lastBlock, pos = lastBlockPos;
  488.     do {
  489.         blockPositions.push_back(pos);  // list is backwards after this...
  490.         pos = matrix[block][pos - queryFrom].tracebackResidue;
  491.         --block;
  492.     } while (pos != NO_TRACEBACK);
  493.     unsigned int firstBlock = block + 1; // last block traced to == first block of the alignment
  495.     // allocate and fill in alignment result structure
  496.     alignment->score = matrix[lastBlock][lastBlockPos - queryFrom].score;
  497.     alignment->firstBlock = firstBlock;
  498.     alignment->nBlocks = blockPositions.size();
  499.     alignment->blockPositions = new unsigned int[blockPositions.size()];
  500.     for (block=0; block<blockPositions.size(); ++block)
  501.         alignment->blockPositions[block] = blockPositions[lastBlock - firstBlock - block];
  503.     return STRUCT_DP_FOUND_ALIGNMENT;
  504. }
  506. int TracebackGlobalAlignment(const Matrix& matrix,
  507.     const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
  508.     DP_AlignmentResult **alignment)
  509. {
  510.     if (!alignment) {
  511.         ERROR_MESSAGE("TracebackGlobalAlignment() - NULL alignment handle");
  512.         return STRUCT_DP_PARAMETER_ERROR;
  513.     }
  514.     *alignment = NULL;
  516.     // find max score (e.g., best-scoring position of last block)
  517.     int score = DP_NEGATIVE_INFINITY;
  518.     unsigned int residue, lastBlockPos = 0;
  519.     for (residue=queryFrom; residue<=queryTo; ++residue) {
  520.         if (matrix[blocks->nBlocks - 1][residue - queryFrom].score > score) {
  521.             score = matrix[blocks->nBlocks - 1][residue - queryFrom].score;
  522.             lastBlockPos = residue;
  523.         }
  524.     }
  526.     if (score == DP_NEGATIVE_INFINITY) {
  527.         ERROR_MESSAGE("TracebackGlobalAlignment() - somehow failed to find any allowed global alignment");
  528.         return STRUCT_DP_ALGORITHM_ERROR;
  529.     }
  531. //    INFO_MESSAGE("Score of best global alignment: " << score);
  532.     *alignment = new DP_AlignmentResult;
  533.     return TracebackAlignment(matrix, blocks->nBlocks - 1, lastBlockPos, queryFrom, *alignment);
  534. }
  536. int TracebackLocalAlignment(const Matrix& matrix,
  537.     const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
  538.     DP_AlignmentResult **alignment)
  539. {
  540.     if (!alignment) {
  541.         ERROR_MESSAGE("TracebackLocalAlignment() - NULL alignment handle");
  542.         return STRUCT_DP_PARAMETER_ERROR;
  543.     }
  544.     *alignment = NULL;
  546.     // find max score (e.g., best-scoring position of any block)
  547.     int score = DP_NEGATIVE_INFINITY;
  548.     unsigned int block, residue, lastBlock = 0, lastBlockPos = 0;
  549.     for (block=0; block<blocks->nBlocks; ++block) {
  550.         for (residue=queryFrom; residue<=queryTo; ++residue) {
  551.             if (matrix[block][residue - queryFrom].score > score) {
  552.                 score = matrix[block][residue - queryFrom].score;
  553.                 lastBlock = block;
  554.                 lastBlockPos = residue;
  555.             }
  556.         }
  557.     }
  559.     if (score <= 0) {
  560. //        INFO_MESSAGE("No positive-scoring local alignment found.");
  561.         return STRUCT_DP_NO_ALIGNMENT;
  562.     }
  564. //    INFO_MESSAGE("Score of best local alignment: " << score);
  565.     *alignment = new DP_AlignmentResult;
  566.     return TracebackAlignment(matrix, lastBlock, lastBlockPos, queryFrom, *alignment);
  567. }
  569. typedef struct {
  570.     unsigned int block;
  571.     unsigned int residue;
  572. } Traceback;
  574. int TracebackMultipleLocalAlignments(const Matrix& matrix,
  575.     const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
  576.     DP_MultipleAlignmentResults **alignments, unsigned int maxAlignments)
  577. {
  578.     if (!alignments) {
  579.         ERROR_MESSAGE("TracebackMultipleLocalAlignments() - NULL alignments handle");
  580.         return STRUCT_DP_PARAMETER_ERROR;
  581.     }
  582.     *alignments = NULL;
  584.     unsigned int block, residue;
  585.     vector < vector < bool > > usedCells(blocks->nBlocks);
  586.     for (block=0; block<blocks->nBlocks; ++block)
  587.         usedCells[block].resize(queryTo - queryFrom + 1, false);
  589.     // find N max scores
  590.     list < Traceback > tracebacks;
  591.     while (maxAlignments == 0 || tracebacks.size() < maxAlignments) {
  593.         // find next best scoring cell that's not part of an alignment already reported
  594.         int score = DP_NEGATIVE_INFINITY;
  595.         unsigned int lastBlock = 0, lastBlockPos = 0;
  596.         for (block=0; block<blocks->nBlocks; ++block) {
  597.             for (residue=queryFrom; residue<=queryTo; ++residue) {
  598.                 if (!usedCells[block][residue - queryFrom] &&
  599.                     matrix[block][residue - queryFrom].score > score)
  600.                 {
  601.                     score = matrix[block][residue - queryFrom].score;
  602.                     lastBlock = block;
  603.                     lastBlockPos = residue;
  604.                 }
  605.             }
  606.         }
  607.         if (score <= 0)
  608.             break;
  610.         // mark cells of this alignment as used, and see if any part of this alignment
  611.         // has been reported before
  612.         block = lastBlock;
  613.         residue = lastBlockPos;
  614.         bool repeated = false;
  615.         do {
  616.             if (usedCells[block][residue - queryFrom]) {
  617.                 repeated = true;
  618.                 break;
  619.             } else {
  620.                 usedCells[block][residue - queryFrom] = true;
  621.             }
  622.             residue = matrix[block][residue - queryFrom].tracebackResidue;
  623.             --block;
  624.         } while (residue != NO_TRACEBACK);
  625.         if (repeated)
  626.             continue;
  628.         // add traceback start to list
  629.         tracebacks.resize(tracebacks.size() + 1);
  630.         tracebacks.back().block = lastBlock;
  631.         tracebacks.back().residue = lastBlockPos;
  632.     }
  634.     if (tracebacks.size() == 0) {
  635. //        INFO_MESSAGE("No positive-scoring local alignment found.");
  636.         return STRUCT_DP_NO_ALIGNMENT;
  637.     }
  639.     // allocate result structure
  640.     *alignments = new DP_MultipleAlignmentResults;
  641.     (*alignments)->nAlignments = 0;
  642.     (*alignments)->alignments = new DP_AlignmentResult[tracebacks.size()];
  643.     unsigned int a;
  644.     for (a=0; a<tracebacks.size(); ++a)
  645.         (*alignments)->alignments[a].blockPositions = NULL;
  647.     // fill in results from tracebacks
  648.     list < Traceback >::const_iterator t, te = tracebacks.end();
  649.     for (t=tracebacks.begin(), a=0; t!=te; ++t, ++a) {
  650.         if (TracebackAlignment(matrix, t->block, t->residue, queryFrom, &((*alignments)->alignments[a]))
  651.                 == STRUCT_DP_FOUND_ALIGNMENT)
  652.         {
  653.             ++((*alignments)->nAlignments);
  654. //            if (a == 0)
  655. //                INFO_MESSAGE("Score of best local alignment: " << (*alignments)->alignments[a].score);
  656. //            else
  657. //                INFO_MESSAGE("Score of next local alignment: " << (*alignments)->alignments[a].score);
  658.         } else
  659.             return STRUCT_DP_ALGORITHM_ERROR;
  660.     }
  662.     return STRUCT_DP_FOUND_ALIGNMENT;
  663. }
  665. END_SCOPE(struct_dp)
  668. // leave the C-called functions outside any namespace, just in case that might
  669. // cause any problems when linking to C code...
  670. USING_SCOPE(struct_dp);
  673. int DP_GlobalBlockAlign(
  674.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
  675.     unsigned int queryFrom, unsigned int queryTo,
  676.     DP_AlignmentResult **alignment)
  677. {
  678.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || !BlockScore || queryTo < queryFrom) {
  679.         ERROR_MESSAGE("DP_GlobalBlockAlign() - invalid parameters");
  680.         return STRUCT_DP_PARAMETER_ERROR;
  681.     }
  683.     unsigned int i, sumBlockLen = 0;
  684.     for (i=0; i<blocks->nBlocks; ++i)
  685.         sumBlockLen += blocks->blockSizes[i];
  686.     if (sumBlockLen > queryTo - queryFrom + 1) {
  687.         ERROR_MESSAGE("DP_GlobalBlockAlign() - sum of block lengths longer than query region");
  688.         return STRUCT_DP_PARAMETER_ERROR;
  689.     }
  691.     int status = ValidateFrozenBlockPositions(blocks, queryFrom, queryTo, true);
  692.     if (status != STRUCT_DP_OKAY) {
  693.         ERROR_MESSAGE("DP_GlobalBlockAlign() - ValidateFrozenBlockPositions() returned error");
  694.         return status;
  695.     }
  697.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  699.     status = CalculateGlobalMatrix(matrix, blocks, BlockScore, queryFrom, queryTo);
  700.     if (status != STRUCT_DP_OKAY) {
  701.         ERROR_MESSAGE("DP_GlobalBlockAlign() - CalculateGlobalMatrix() failed");
  702.         return status;
  703.     }
  705.     return TracebackGlobalAlignment(matrix, blocks, queryFrom, queryTo, alignment);
  706. }
  708. int DP_GlobalBlockAlignGeneric(
  709.     const DP_BlockInfo *blocks,
  710.     DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
  711.     unsigned int queryFrom, unsigned int queryTo,
  712.     DP_AlignmentResult **alignment)
  713. {
  714.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || queryTo < queryFrom ||
  715.             !BlockScore || !LoopScore) {
  716.         ERROR_MESSAGE("DP_GlobalBlockAlignGeneric() - invalid parameters");
  717.         return STRUCT_DP_PARAMETER_ERROR;
  718.     }
  720.     unsigned int i, sumBlockLen = 0;
  721.     for (i=0; i<blocks->nBlocks; ++i)
  722.         sumBlockLen += blocks->blockSizes[i];
  723.     if (sumBlockLen > queryTo - queryFrom + 1) {
  724.         ERROR_MESSAGE("DP_GlobalBlockAlignGeneric() - sum of block lengths longer than query region");
  725.         return STRUCT_DP_PARAMETER_ERROR;
  726.     }
  728.     int status = ValidateFrozenBlockPositions(blocks, queryFrom, queryTo, false);
  729.     if (status != STRUCT_DP_OKAY) {
  730.         ERROR_MESSAGE("DP_GlobalBlockAlignGeneric() - ValidateFrozenBlockPositions() returned error");
  731.         return status;
  732.     }
  734.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  736.     status = CalculateGlobalMatrixGeneric(matrix, blocks, BlockScore, LoopScore, queryFrom, queryTo);
  737.     if (status != STRUCT_DP_OKAY) {
  738.         ERROR_MESSAGE("DP_GlobalBlockAlignGeneric() - CalculateGlobalMatrixGeneric() failed");
  739.         return status;
  740.     }
  742.     return TracebackGlobalAlignment(matrix, blocks, queryFrom, queryTo, alignment);
  743. }
  745. int DP_LocalBlockAlign(
  746.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
  747.     unsigned int queryFrom, unsigned int queryTo,
  748.     DP_AlignmentResult **alignment)
  749. {
  750.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || !BlockScore || queryTo < queryFrom) {
  751.         ERROR_MESSAGE("DP_LocalBlockAlign() - invalid parameters");
  752.         return STRUCT_DP_PARAMETER_ERROR;
  753.     }
  754.     for (unsigned int block=0; block<blocks->nBlocks; ++block) {
  755.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  756.             WARNING_MESSAGE("DP_LocalBlockAlign() - frozen block specifications are ignored...");
  757.             break;
  758.         }
  759.     }
  761.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  763.     int status = CalculateLocalMatrix(matrix, blocks, BlockScore, queryFrom, queryTo);
  764.     if (status != STRUCT_DP_OKAY) {
  765.         ERROR_MESSAGE("DP_LocalBlockAlign() - CalculateLocalMatrix() failed");
  766.         return status;
  767.     }
  769.     return TracebackLocalAlignment(matrix, blocks, queryFrom, queryTo, alignment);
  770. }
  772. int DP_LocalBlockAlignGeneric(
  773.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
  774.     unsigned int queryFrom, unsigned int queryTo,
  775.     DP_AlignmentResult **alignment)
  776. {
  777.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || !BlockScore || queryTo < queryFrom) {
  778.         ERROR_MESSAGE("DP_LocalBlockAlignGeneric() - invalid parameters");
  779.         return STRUCT_DP_PARAMETER_ERROR;
  780.     }
  781.     for (unsigned int block=0; block<blocks->nBlocks; ++block) {
  782.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  783.             WARNING_MESSAGE("DP_LocalBlockAlignGeneric() - frozen block specifications are ignored...");
  784.             break;
  785.         }
  786.     }
  788.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  790.     int status = CalculateLocalMatrixGeneric(matrix, blocks, BlockScore, LoopScore, queryFrom, queryTo);
  791.     if (status != STRUCT_DP_OKAY) {
  792.         ERROR_MESSAGE("DP_LocalBlockAlignGeneric() - CalculateLocalMatrixGeneric() failed");
  793.         return status;
  794.     }
  796.     return TracebackLocalAlignment(matrix, blocks, queryFrom, queryTo, alignment);
  797. }
  799. int DP_MultipleLocalBlockAlign(
  800.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
  801.     unsigned int queryFrom, unsigned int queryTo,
  802.     DP_MultipleAlignmentResults **alignments, unsigned int maxAlignments)
  803. {
  804.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || !BlockScore || queryTo < queryFrom) {
  805.         ERROR_MESSAGE("DP_MultipleLocalBlockAlign() - invalid parameters");
  806.         return STRUCT_DP_PARAMETER_ERROR;
  807.     }
  808.     for (unsigned int block=0; block<blocks->nBlocks; ++block) {
  809.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  810.             WARNING_MESSAGE("DP_MultipleLocalBlockAlign() - frozen block specifications are ignored...");
  811.             break;
  812.         }
  813.     }
  815.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  817.     int status = CalculateLocalMatrix(matrix, blocks, BlockScore, queryFrom, queryTo);
  818.     if (status != STRUCT_DP_OKAY) {
  819.         ERROR_MESSAGE("DP_MultipleLocalBlockAlign() - CalculateLocalMatrix() failed");
  820.         return status;
  821.     }
  823.     return TracebackMultipleLocalAlignments(matrix, blocks, queryFrom, queryTo, alignments, maxAlignments);
  824. }
  826. int DP_MultipleLocalBlockAlignGeneric(
  827.     const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
  828.     unsigned int queryFrom, unsigned int queryTo,
  829.     DP_MultipleAlignmentResults **alignments, unsigned int maxAlignments)
  830. {
  831.     if (!blocks || blocks->nBlocks < 1 || !blocks->blockSizes || !BlockScore || queryTo < queryFrom) {
  832.         ERROR_MESSAGE("DP_MultipleLocalBlockAlignGeneric() - invalid parameters");
  833.         return STRUCT_DP_PARAMETER_ERROR;
  834.     }
  835.     for (unsigned int block=0; block<blocks->nBlocks; ++block) {
  836.         if (blocks->freezeBlocks[block] != DP_UNFROZEN_BLOCK) {
  837.             WARNING_MESSAGE("DP_MultipleLocalBlockAlignGeneric() - frozen block specifications are ignored...");
  838.             break;
  839.         }
  840.     }
  842.     Matrix matrix(blocks->nBlocks, queryTo - queryFrom + 1);
  844.     int status = CalculateLocalMatrixGeneric(matrix, blocks, BlockScore, LoopScore, queryFrom, queryTo);
  845.     if (status != STRUCT_DP_OKAY) {
  846.         ERROR_MESSAGE("DP_MultipleLocalBlockAlignGeneric() - CalculateLocalMatrixGeneric() failed");
  847.         return status;
  848.     }
  850.     return TracebackMultipleLocalAlignments(matrix, blocks, queryFrom, queryTo, alignments, maxAlignments);
  851. }
  853. DP_BlockInfo * DP_CreateBlockInfo(unsigned int nBlocks)
  854. {
  855.     DP_BlockInfo *blocks = new DP_BlockInfo;
  856.     blocks->nBlocks = nBlocks;
  857.     blocks->blockPositions = new unsigned int[nBlocks];
  858.     blocks->blockSizes = new unsigned int[nBlocks];
  859.     blocks->maxLoops = new unsigned int[nBlocks - 1];
  860.     blocks->freezeBlocks = new unsigned int[nBlocks];
  861.     for (unsigned int i=0; i<nBlocks; ++i)
  862.         blocks->freezeBlocks[i] = DP_UNFROZEN_BLOCK;
  863.     return blocks;
  864. }
  866. void DP_DestroyBlockInfo(DP_BlockInfo *blocks)
  867. {
  868.     if (!blocks)
  869.         return;
  870.     delete[] blocks->blockPositions;
  871.     delete[] blocks->blockSizes;
  872.     delete[] blocks->maxLoops;
  873.     delete[] blocks->freezeBlocks;
  874.     delete blocks;
  875. }
  877. void DP_DestroyAlignmentResult(DP_AlignmentResult *alignment)
  878. {
  879.     if (!alignment)
  880.         return;
  881.     delete[] alignment->blockPositions;
  882.     delete alignment;
  883. }
  885. void DP_DestroyMultipleAlignmentResults(DP_MultipleAlignmentResults *alignments)
  886. {
  887.     if (!alignments) return;
  888.     for (unsigned int i=0; i<alignments->nAlignments; ++i)
  889.         if (alignments->alignments[i].blockPositions)
  890.             delete[] alignments->alignments[i].blockPositions;
  891.     delete[] alignments->alignments;
  892.     delete alignments;
  893. }
  895. unsigned int DP_CalculateMaxLoopLength(
  896.         unsigned int nLoops, const unsigned int *loops,
  897.         double percentile, unsigned int extension, unsigned int cutoff)
  898. {
  899.     vector < unsigned int > loopLengths(nLoops);
  900.     unsigned int index, max;
  901.     for (index=0; index<nLoops; ++index)
  902.         loopLengths[index] = loops[index];
  904.     stable_sort(loopLengths.begin(), loopLengths.end());
  906.     if (percentile < 1.0) {
  907.         index = (unsigned int)(percentile * (nLoops - 1) + 0.5);
  908.         max = loopLengths[index] + extension;
  909.     } else {  // percentile >= 1.0
  910.         max = ((unsigned int)(percentile * loopLengths[nLoops - 1] + 0.5)) + extension;
  911.     }
  913.     if (cutoff > 0 && max > cutoff)
  914.         max = cutoff;
  916.     return max;
  917. }
  919. /*
  920. * ---------------------------------------------------------------------------
  921. * $Log: block_align.cpp,v $
  922. * Revision 1000.2  2004/06/01 18:11:09  gouriano
  923. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  924. *
  925. * Revision 1.20  2004/05/21 21:41:04  gorelenk
  926. * Added PCH ncbi_pch.hpp
  927. *
  928. * Revision 1.19  2004/03/15 18:54:57  thiessen
  929. * prefer prefix vs. postfix ++/-- operators
  930. *
  931. * Revision 1.18  2004/02/19 02:21:19  thiessen
  932. * fix struct_dp paths
  933. *
  934. * Revision 1.17  2003/12/19 14:37:50  thiessen
  935. * add local generic loop function alignment routines
  936. *
  937. * Revision 1.16  2003/12/08 16:33:58  thiessen
  938. * fix signed/unsigned mix
  939. *
  940. * Revision 1.15  2003/12/08 16:21:35  thiessen
  941. * add generic loop scoring function interface
  942. *
  943. * Revision 1.14  2003/09/07 01:11:25  thiessen
  944. * fix small memory leak
  945. *
  946. * Revision 1.13  2003/09/07 00:06:19  thiessen
  947. * add multiple tracebacks for local alignments
  948. *
  949. * Revision 1.12  2003/08/23 22:10:05  thiessen
  950. * fix local alignment block edge bug
  951. *
  952. * Revision 1.11  2003/08/22 14:28:49  thiessen
  953. * add standard loop calculating function
  954. *
  955. * Revision 1.10  2003/08/19 19:36:48  thiessen
  956. * avoid annoying 'might be used uninitialized' warnings
  957. *
  958. * Revision 1.9  2003/08/19 19:26:21  thiessen
  959. * error if block size > query range length
  960. *
  961. * Revision 1.8  2003/07/11 15:27:48  thiessen
  962. * add DP_ prefix to globals
  963. *
  964. * Revision 1.7  2003/06/22 12:18:16  thiessen
  965. * fixes for unsigned/signed comparison
  966. *
  967. * Revision 1.6  2003/06/22 12:11:37  thiessen
  968. * add local alignment
  969. *
  970. * Revision 1.5  2003/06/19 13:48:23  thiessen
  971. * cosmetic/typo fixes
  972. *
  973. * Revision 1.4  2003/06/18 21:55:15  thiessen
  974. * remove unused params
  975. *
  976. * Revision 1.3  2003/06/18 21:46:09  thiessen
  977. * add traceback, alignment return structure
  978. *
  979. * Revision 1.2  2003/06/18 19:10:17  thiessen
  980. * fix lf issues
  981. *
  982. * Revision 1.1  2003/06/18 16:54:12  thiessen
  983. * initial checkin; working global block aligner and demo app
  984. *
  985. */