block_align.cpp

ncbi源码

int CalculateLocalMatrix(Matrix& matrix,
    const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore,
    unsigned int queryFrom, unsigned int queryTo)
{
    unsigned int block, residue, prevResidue, lastBlock = blocks->nBlocks - 1, tracebackResidue = 0;
    int score, sum, bestPrevScore;

    // find last possible block positions, based purely on block lengths
    vector < unsigned int > lastPos(blocks->nBlocks);
    for (block=0; block<=lastBlock; ++block) {
        if (blocks->blockSizes[block] > queryTo - queryFrom + 1) {
            ERROR_MESSAGE("Block " << (block+1) << " too large for this query range");
            return STRUCT_DP_PARAMETER_ERROR;
        }
        lastPos[block] = queryTo - blocks->blockSizes[block] + 1;
    }

    // first row: positive scores of first block at all possible positions
    for (residue=queryFrom; residue<=lastPos[0]; ++residue) {
        score = BlockScore(0, residue);
        matrix[0][residue - queryFrom].score = (score > 0) ? score : 0;
    }

    // first column: positive scores of all blocks at first positions
    for (block=1; block<=lastBlock; ++block) {
        score = BlockScore(block, queryFrom);
        matrix[block][0].score = (score > 0) ? score : 0;
    }

    // for each successive block, find the best positive scoring with a previous block, if any
    for (block=1; block<=lastBlock; ++block) {
        for (residue=queryFrom+1; residue<=lastPos[block]; ++residue) {

            // get score at this position
            score = BlockScore(block, residue);
            if (score == DP_NEGATIVE_INFINITY)
                continue;

            // find max score of any allowed previous block
            bestPrevScore = DP_NEGATIVE_INFINITY;
            for (prevResidue=queryFrom; prevResidue<=lastPos[block - 1]; ++prevResidue) {

                // current block must come after the previous block
                if (residue < prevResidue + blocks->blockSizes[block - 1])
                    break;

                // cut off at max loop length from previous block
                if (residue > prevResidue + blocks->blockSizes[block - 1] + blocks->maxLoops[block - 1])
                    continue;

                // keep maximum score
                if (matrix[block - 1][prevResidue - queryFrom].score > bestPrevScore) {
                    bestPrevScore = matrix[block - 1][prevResidue - queryFrom].score;
                    tracebackResidue = prevResidue;
                }
            }

            // extend alignment if the sum of this block's + previous block's score is positive
            if (bestPrevScore > 0 && (sum=bestPrevScore+score) > 0) {
                matrix[block][residue - queryFrom].score = sum;
                matrix[block][residue - queryFrom].tracebackResidue = tracebackResidue;
            }

            // otherwise, start new alignment if score is positive
            else if (score > 0)
                matrix[block][residue - queryFrom].score = score;
        }
    }

    return STRUCT_DP_OKAY;
}

// fill matrix values; return true on success. Matrix must have only default values when passed in.
int CalculateLocalMatrixGeneric(Matrix& matrix,
    const DP_BlockInfo *blocks, DP_BlockScoreFunction BlockScore, DP_LoopPenaltyFunction LoopScore,
    unsigned int queryFrom, unsigned int queryTo)
{
    unsigned int block, residue, prevResidue, loopPenalty,
        lastBlock = blocks->nBlocks - 1, tracebackResidue = 0;
    int score, sum, bestPrevScore;

    // find last possible block positions, based purely on block lengths
    vector < unsigned int > lastPos(blocks->nBlocks);
    for (block=0; block<=lastBlock; ++block) {
        if (blocks->blockSizes[block] > queryTo - queryFrom + 1) {
            ERROR_MESSAGE("Block " << (block+1) << " too large for this query range");
            return STRUCT_DP_PARAMETER_ERROR;
        }
        lastPos[block] = queryTo - blocks->blockSizes[block] + 1;
    }

    // first row: positive scores of first block at all possible positions
    for (residue=queryFrom; residue<=lastPos[0]; ++residue) {
        score = BlockScore(0, residue);
        matrix[0][residue - queryFrom].score = (score > 0) ? score : 0;
    }

    // first column: positive scores of all blocks at first positions
    for (block=1; block<=lastBlock; ++block) {
        score = BlockScore(block, queryFrom);
        matrix[block][0].score = (score > 0) ? score : 0;
    }

    // for each successive block, find the best positive scoring with a previous block, if any
    for (block=1; block<=lastBlock; ++block) {
        for (residue=queryFrom+1; residue<=lastPos[block]; ++residue) {

            // get score at this position
            score = BlockScore(block, residue);
            if (score == DP_NEGATIVE_INFINITY)
                continue;

            // find max score of any allowed previous block
            bestPrevScore = DP_NEGATIVE_INFINITY;
            for (prevResidue=queryFrom; prevResidue<=lastPos[block - 1]; ++prevResidue) {

                // current block must come after the previous block
                if (residue < prevResidue + blocks->blockSizes[block - 1])
                    break;

                // make sure previous block is at an allowed position
                if (matrix[block - 1][prevResidue - queryFrom].score == DP_NEGATIVE_INFINITY)
                    continue;

                // get loop score
                loopPenalty = LoopScore(block - 1, residue - prevResidue - blocks->blockSizes[block - 1]);
                if (loopPenalty == DP_POSITIVE_INFINITY)
                    continue;

                // keep maximum score
                sum = matrix[block - 1][prevResidue - queryFrom].score - loopPenalty;
                if (sum > bestPrevScore) {
                    bestPrevScore = sum;
                    tracebackResidue = prevResidue;
                }
            }

            // extend alignment if the sum of this block's + previous block's score is positive
            if (bestPrevScore > 0 && (sum=bestPrevScore+score) > 0) {
                matrix[block][residue - queryFrom].score = sum;
                matrix[block][residue - queryFrom].tracebackResidue = tracebackResidue;
            }

            // otherwise, start new alignment if score is positive
            else if (score > 0)
                matrix[block][residue - queryFrom].score = score;
        }
    }

    return STRUCT_DP_OKAY;
}

int TracebackAlignment(const Matrix& matrix, unsigned int lastBlock, unsigned int lastBlockPos,
    unsigned int queryFrom, DP_AlignmentResult *alignment)
{
    // trace backwards from last block/pos
    vector < unsigned int > blockPositions;
    unsigned int block = lastBlock, pos = lastBlockPos;
    do {
        blockPositions.push_back(pos);  // list is backwards after this...
        pos = matrix[block][pos - queryFrom].tracebackResidue;
        --block;
    } while (pos != NO_TRACEBACK);
    unsigned int firstBlock = block + 1; // last block traced to == first block of the alignment

    // allocate and fill in alignment result structure
    alignment->score = matrix[lastBlock][lastBlockPos - queryFrom].score;
    alignment->firstBlock = firstBlock;
    alignment->nBlocks = blockPositions.size();
    alignment->blockPositions = new unsigned int[blockPositions.size()];
    for (block=0; block<blockPositions.size(); ++block)
        alignment->blockPositions[block] = blockPositions[lastBlock - firstBlock - block];

    return STRUCT_DP_FOUND_ALIGNMENT;
}

int TracebackGlobalAlignment(const Matrix& matrix,
    const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
    DP_AlignmentResult **alignment)
{
    if (!alignment) {
        ERROR_MESSAGE("TracebackGlobalAlignment() - NULL alignment handle");
        return STRUCT_DP_PARAMETER_ERROR;
    }
    *alignment = NULL;

    // find max score (e.g., best-scoring position of last block)
    int score = DP_NEGATIVE_INFINITY;
    unsigned int residue, lastBlockPos = 0;
    for (residue=queryFrom; residue<=queryTo; ++residue) {
        if (matrix[blocks->nBlocks - 1][residue - queryFrom].score > score) {
            score = matrix[blocks->nBlocks - 1][residue - queryFrom].score;
            lastBlockPos = residue;
        }
    }

    if (score == DP_NEGATIVE_INFINITY) {
        ERROR_MESSAGE("TracebackGlobalAlignment() - somehow failed to find any allowed global alignment");
        return STRUCT_DP_ALGORITHM_ERROR;
    }

//    INFO_MESSAGE("Score of best global alignment: " << score);
    *alignment = new DP_AlignmentResult;
    return TracebackAlignment(matrix, blocks->nBlocks - 1, lastBlockPos, queryFrom, *alignment);
}

int TracebackLocalAlignment(const Matrix& matrix,
    const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
    DP_AlignmentResult **alignment)
{
    if (!alignment) {
        ERROR_MESSAGE("TracebackLocalAlignment() - NULL alignment handle");
        return STRUCT_DP_PARAMETER_ERROR;
    }
    *alignment = NULL;

    // find max score (e.g., best-scoring position of any block)
    int score = DP_NEGATIVE_INFINITY;
    unsigned int block, residue, lastBlock = 0, lastBlockPos = 0;
    for (block=0; block<blocks->nBlocks; ++block) {
        for (residue=queryFrom; residue<=queryTo; ++residue) {
            if (matrix[block][residue - queryFrom].score > score) {
                score = matrix[block][residue - queryFrom].score;
                lastBlock = block;
                lastBlockPos = residue;
            }
        }
    }

    if (score <= 0) {
//        INFO_MESSAGE("No positive-scoring local alignment found.");
        return STRUCT_DP_NO_ALIGNMENT;
    }

//    INFO_MESSAGE("Score of best local alignment: " << score);
    *alignment = new DP_AlignmentResult;
    return TracebackAlignment(matrix, lastBlock, lastBlockPos, queryFrom, *alignment);
}

typedef struct {
    unsigned int block;
    unsigned int residue;
} Traceback;

int TracebackMultipleLocalAlignments(const Matrix& matrix,
    const DP_BlockInfo *blocks, unsigned int queryFrom, unsigned int queryTo,
    DP_MultipleAlignmentResults **alignments, unsigned int maxAlignments)
{
    if (!alignments) {
        ERROR_MESSAGE("TracebackMultipleLocalAlignments() - NULL alignments handle");
        return STRUCT_DP_PARAMETER_ERROR;
    }
    *alignments = NULL;

    unsigned int block, residue;
    vector < vector < bool > > usedCells(blocks->nBlocks);
    for (block=0; block<blocks->nBlocks; ++block)
        usedCells[block].resize(queryTo - queryFrom + 1, false);

    // find N max scores
    list < Traceback > tracebacks;
    while (maxAlignments == 0 || tracebacks.size() < maxAlignments) {

        // find next best scoring cell that's not part of an alignment already reported
        int score = DP_NEGATIVE_INFINITY;
        unsigned int lastBlock = 0, lastBlockPos = 0;
        for (block=0; block<blocks->nBlocks; ++block) {
            for (residue=queryFrom; residue<=queryTo; ++residue) {
                if (!usedCells[block][residue - queryFrom] &&
                    matrix[block][residue - queryFrom].score > score)
                {
                    score = matrix[block][residue - queryFrom].score;
                    lastBlock = block;
                    lastBlockPos = residue;
                }
            }
        }
        if (score <= 0)
            break;

        // mark cells of this alignment as used, and see if any part of this alignment
        // has been reported before
        block = lastBlock;
        residue = lastBlockPos;
        bool repeated = false;
        do {
            if (usedCells[block][residue - queryFrom]) {
                repeated = true;
                break;
            } else {
                usedCells[block][residue - queryFrom] = true;
            }
            residue = matrix[block][residue - queryFrom].tracebackResidue;
            --block;
        } while (residue != NO_TRACEBACK);
        if (repeated)
            continue;

        // add traceback start to list
        tracebacks.resize(tracebacks.size() + 1);
        tracebacks.back().block = lastBlock;
        tracebacks.back().residue = lastBlockPos;
    }

    if (tracebacks.size() == 0) {
//        INFO_MESSAGE("No positive-scoring local alignment found.");
        return STRUCT_DP_NO_ALIGNMENT;
    }

    // allocate result structure
    *alignments = new DP_MultipleAlignmentResults;
    (*alignments)->nAlignments = 0;
    (*alignments)->alignments = new DP_AlignmentResult[tracebacks.size()];
    unsigned int a;
    for (a=0; a<tracebacks.size(); ++a)
        (*alignments)->alignments[a].blockPositions = NULL;

    // fill in results from tracebacks
    list < Traceback >::const_iterator t, te = tracebacks.end();
    for (t=tracebacks.begin(), a=0; t!=te; ++t, ++a) {
        if (TracebackAlignment(matrix, t->block, t->residue, queryFrom, &((*alignments)->alignments[a]))
                == STRUCT_DP_FOUND_ALIGNMENT)
        {
            ++((*alignments)->nAlignments);
//            if (a == 0)
//                INFO_MESSAGE("Score of best local alignment: " << (*alignments)->alignments[a].score);
//            else
//                INFO_MESSAGE("Score of next local alignment: " << (*alignments)->alignments[a].score);
        } else