cn3d_ba_interface.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: cn3d_ba_interface.cpp,v $
 * PRODUCTION Revision 1000.2  2004/06/01 18:28:07  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.33
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: cn3d_ba_interface.cpp,v 1000.2 2004/06/01 18:28:07 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================
*
* Authors:  Paul Thiessen
*
* File Description:
*      Interface to Alejandro's block alignment algorithm
*
* ===========================================================================
*/

#ifdef _MSC_VER
#pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
#endif

#include <ncbi_pch.hpp>
#include <corelib/ncbistd.hpp>
#include <corelib/ncbi_limits.h>

#include <objects/seqalign/Seq_align.hpp>
#include <objects/seqalign/Seq_align_set.hpp>
#include <objects/seqalign/Dense_seg.hpp>
#include <objects/seqalign/Score.hpp>
#include <objects/general/Object_id.hpp>

#ifdef __WXMSW__
#include <windows.h>
#include <wx/msw/winundef.h>
#endif
#include <wx/wx.h>

#include <algo/structure/struct_dp/struct_dp.h>

#include "cn3d_ba_interface.hpp"
#include "block_multiple_alignment.hpp"
#include "sequence_set.hpp"
#include "structure_set.hpp"
#include "asn_converter.hpp"
#include "molecule_identifier.hpp"
#include "wx_tools.hpp"
#include "cn3d_tools.hpp"
#include "cn3d_blast.hpp"

USING_NCBI_SCOPE;
USING_SCOPE(objects);


BEGIN_SCOPE(Cn3D)

class BlockAlignerOptionsDialog : public wxDialog
{
public:
    BlockAlignerOptionsDialog(wxWindow* parent, const BlockAligner::BlockAlignerOptions& init);
    ~BlockAlignerOptionsDialog(void);

    bool GetValues(BlockAligner::BlockAlignerOptions *options);

private:
    IntegerSpinCtrl *iExtension, *iCutoff;
    FloatingPointSpinCtrl *fpPercent;
    wxCheckBox *cGlobal, *cKeep, *cMerge;

    void OnCloseWindow(wxCloseEvent& event);
    void OnButton(wxCommandEvent& event);
    DECLARE_EVENT_TABLE()
};

BlockAligner::BlockAligner(void)
{
    // default options
    currentOptions.loopPercentile = 0.6;
    currentOptions.loopExtension = 10;
    currentOptions.loopCutoff = 0;
    currentOptions.globalAlignment = false;
    currentOptions.keepExistingBlocks = false;
    currentOptions.mergeAfterEachSequence = false;
}

static void FreezeBlocks(const BlockMultipleAlignment *multiple,
    const BlockMultipleAlignment *pairwise, DP_BlockInfo *blockInfo)
{
    BlockMultipleAlignment::UngappedAlignedBlockList multipleABlocks, pairwiseABlocks;
    multiple->GetUngappedAlignedBlocks(&multipleABlocks);
    pairwise->GetUngappedAlignedBlocks(&pairwiseABlocks);

    // if a block in the multiple is contained in the pairwise (looking at master coords),
    // then add a constraint to keep it there
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator
        m, me = multipleABlocks.end(), p, pe = pairwiseABlocks.end();
    const Block::Range *multipleRangeMaster, *pairwiseRangeMaster, *pairwiseRangeSlave;
    int i;
    for (i=0, m=multipleABlocks.begin(); m!=me; ++i, ++m) {
        multipleRangeMaster = (*m)->GetRangeOfRow(0);
        for (p=pairwiseABlocks.begin(); p!=pe; ++p) {
            pairwiseRangeMaster = (*p)->GetRangeOfRow(0);
            if (pairwiseRangeMaster->from <= multipleRangeMaster->from &&
                    pairwiseRangeMaster->to >= multipleRangeMaster->to) {
                pairwiseRangeSlave = (*p)->GetRangeOfRow(1);
                blockInfo->freezeBlocks[i] = pairwiseRangeSlave->from +
                    (multipleRangeMaster->from - pairwiseRangeMaster->from);
                break;
            }
        }
        if (p == pe)
            blockInfo->freezeBlocks[i] = -1;
//        if (blockInfo->freezeBlocks[i] >= 0)
//            TESTMSG("block " << (i+1) << " frozen at query pos " << (blockInfo->freezeBlocks[i]+1));
//        else
//            TESTMSG("block " << (i+1) << " unfrozen");
    }
}

// global stuff for DP block aligner score callback
DP_BlockInfo *dpBlocks = NULL;
const BLAST_Matrix *dpPSSM = NULL;
const Sequence *dpQuery = NULL;

// sum of scores for residue vs. PSSM
int dpScoreFunction(unsigned int block, unsigned int queryPos)
{
    if (!dpBlocks || !dpPSSM || !dpQuery || block >= dpBlocks->nBlocks ||
        queryPos > dpQuery->Length() - dpBlocks->blockSizes[block])
    {
        ERRORMSG("dpScoreFunction() - bad parameters");
        return DP_NEGATIVE_INFINITY;
    }

    int i, masterPos = dpBlocks->blockPositions[block], score = 0;
    for (i=0; i<dpBlocks->blockSizes[block]; ++i)
        score += dpPSSM->matrix[masterPos + i]
            [LookupBLASTResidueNumberFromCharacter(dpQuery->sequenceString[queryPos + i])];

    return score;
}

static BlockMultipleAlignment * UnpackDPResult(DP_BlockInfo *blocks, DP_AlignmentResult *result,
    const Sequence *master, const Sequence *query)
{
    // create new alignment structure
    BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
    (*seqs)[0] = master;
    (*seqs)[1] = query;
    auto_ptr<BlockMultipleAlignment>
        bma(new BlockMultipleAlignment(seqs, master->parentSet->alignmentManager));

    // unpack result blocks
    for (unsigned int b=0; b<result->nBlocks; ++b) {
        UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(bma.get());
        newBlock->width = blocks->blockSizes[b + result->firstBlock];
        newBlock->SetRangeOfRow(0,
            blocks->blockPositions[b + result->firstBlock],
            blocks->blockPositions[b + result->firstBlock] + newBlock->width - 1);
        newBlock->SetRangeOfRow(1,
            result->blockPositions[b],
            result->blockPositions[b] + newBlock->width - 1);
        bma->AddAlignedBlockAtEnd(newBlock);
        TRACEMSG("block " << (b+1)
            << " position: " << (result->blockPositions[b]+1)
            << " score: " << dpScoreFunction(b, result->blockPositions[b]));
    }

    // finalize the alignment
    if (!bma->AddUnalignedBlocks() || !bma->UpdateBlockMapAndColors(false)) {
        ERRORMSG("Error finalizing new alignment!");
        return NULL;
    }

    // get scores
    wxString score;
    score.Printf(" raw score: %i", result->score);
    bma->SetRowStatusLine(0, score.c_str());
    bma->SetRowStatusLine(1, score.c_str());

    // success
    return bma.release();
}

bool BlockAligner::CreateNewPairwiseAlignmentsByBlockAlignment(BlockMultipleAlignment *multiple,
    const AlignmentList& toRealign, AlignmentList *newAlignments, int *nRowsAddedToMultiple,
    bool canMerge)
{
    // show options dialog each time block aligner is run
    if (!SetOptions(NULL)) return false;
    if (currentOptions.mergeAfterEachSequence && !canMerge) {
        ERRORMSG("Can only merge when editing is enabled in the sequence window");
        return false;
    }

    newAlignments->clear();
    BlockMultipleAlignment::UngappedAlignedBlockList blocks;
    multiple->GetUngappedAlignedBlocks(&blocks);
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
    if (nRowsAddedToMultiple) *nRowsAddedToMultiple = 0;
    int i;

    // set up block info
    dpBlocks = DP_CreateBlockInfo(blocks.size());
    unsigned int *loopLengths = new unsigned int[multiple->NRows()];
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator n;
    const Block::Range *range, *nextRange;
    for (i=0, b=blocks.begin(); b!=be; ++i, ++b) {
        range = (*b)->GetRangeOfRow(0);
        dpBlocks->blockPositions[i] = range->from;
        dpBlocks->blockSizes[i] = range->to - range->from + 1;

        // calculate max loop lengths
        if (i < blocks.size() - 1) {
            n = b;
            ++n;
            for (int r=0; r<multiple->NRows(); ++r) {
                range = (*b)->GetRangeOfRow(r);
                nextRange = (*n)->GetRangeOfRow(r);
                loopLengths[r] = nextRange->from - range->to - 1;
            }
            dpBlocks->maxLoops[i] = DP_CalculateMaxLoopLength(multiple->NRows(), loopLengths,
                currentOptions.loopPercentile, currentOptions.loopExtension, currentOptions.loopCutoff);
            TRACEMSG("allowed gap after block " << (i+1) << ": " << dpBlocks->maxLoops[i]);
        }
    }
    delete loopLengths;

    // set up PSSM
    dpPSSM = multiple->GetPSSM();

    bool errorsEncountered = false;

    AlignmentList::const_iterator s, se = toRealign.end();
    for (s=toRealign.begin(); s!=se; ++s) {
        if (multiple && (*s)->GetMaster() != multiple->GetMaster())
            ERRORMSG("master sequence mismatch");

        // slave sequence info
        dpQuery = (*s)->GetSequenceOfRow(1);
        int startQueryPosition =
            ((*s)->alignSlaveFrom >= 0 && (*s)->alignSlaveFrom < dpQuery->Length()) ?
                (*s)->alignSlaveFrom : 0;
        int endQueryPosition =
            ((*s)->alignSlaveTo >= 0 && (*s)->alignSlaveTo < dpQuery->Length()) ?
                (*s)->alignSlaveTo : dpQuery->Length() - 1;
        TRACEMSG("query region: " << (startQueryPosition+1) << " to " << (endQueryPosition+1));

        // set frozen blocks
        if (!currentOptions.keepExistingBlocks) {
            for (i=0; i<blocks.size(); ++i)
                dpBlocks->freezeBlocks[i] = DP_UNFROZEN_BLOCK;
        } else {
            FreezeBlocks(multiple, *s, dpBlocks);
        }

        SetDialogSevereErrors(false);

        // actually do the block alignment
        DP_AlignmentResult *dpResult = NULL;
        int dpStatus;
        if (currentOptions.globalAlignment)
            dpStatus = DP_GlobalBlockAlign(dpBlocks, dpScoreFunction,
                startQueryPosition, endQueryPosition, &dpResult);
        else
            dpStatus = DP_LocalBlockAlign(dpBlocks, dpScoreFunction,
                startQueryPosition, endQueryPosition, &dpResult);

        SetDialogSevereErrors(true);

        // create new alignment from DP result
        BlockMultipleAlignment *dpAlignment = NULL;
        if (dpResult) {
            dpAlignment = UnpackDPResult(dpBlocks, dpResult, multiple->GetMaster(), dpQuery);
            if (!dpAlignment)
                ERRORMSG("Couldn't create BlockMultipleAlignment from DP result");
        }

        if (dpStatus == STRUCT_DP_FOUND_ALIGNMENT && dpAlignment) {

            // merge or add alignment to list
            if (currentOptions.mergeAfterEachSequence && multiple->MergeAlignment(dpAlignment)) {
                delete dpAlignment; // if merge is successful, we can delete this alignment;