block_multiple_alignment.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: block_multiple_alignment.cpp,v $
 * PRODUCTION Revision 1000.3  2004/06/01 18:27:49  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.57
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: block_multiple_alignment.cpp,v 1000.3 2004/06/01 18:27:49 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:
*      Classes to hold alignment data
*
* ===========================================================================
*/

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

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

#include <objects/seqalign/Dense_diag.hpp>

#include "block_multiple_alignment.hpp"
#include "sequence_set.hpp"
#include "molecule.hpp"
#include "conservation_colorer.hpp"
#include "style_manager.hpp"
#include "structure_set.hpp"
#include "messenger.hpp"
#include "cn3d_colors.hpp"
#include "alignment_manager.hpp"
#include "cn3d_tools.hpp"
#include "molecule_identifier.hpp"
#include "cn3d_threader.hpp"
#include "cn3d_blast.hpp"

// hack so I can catch memory leaks specific to this module, at the line where allocation occurs
#ifdef _DEBUG
#ifdef MemNew
#undef MemNew
#endif
#define MemNew(sz) memset(malloc(sz), 0, (sz))
#endif

USING_NCBI_SCOPE;
USING_SCOPE(objects);


BEGIN_SCOPE(Cn3D)

BlockMultipleAlignment::BlockMultipleAlignment(SequenceList *sequenceList, AlignmentManager *alnMgr) :
    sequences(sequenceList), conservationColorer(NULL), alignmentManager(alnMgr),
    alignMasterFrom(-1), alignMasterTo(-1), alignSlaveFrom(-1), alignSlaveTo(-1),
    pssm(NULL), showGeometryViolations(false)
{
    InitCache();
    rowDoubles.resize(sequenceList->size(), 0.0);
    rowStrings.resize(sequenceList->size());

    // create conservation colorer
    conservationColorer = new ConservationColorer(this);
}

void BlockMultipleAlignment::InitCache(void)
{
    cachePrevRow = -1;
    cachePrevBlock = NULL;
    cacheBlockIterator = blocks.begin();
}

BlockMultipleAlignment::~BlockMultipleAlignment(void)
{
    BlockList::iterator i, ie = blocks.end();
    for (i=blocks.begin(); i!=ie; ++i) delete *i;
    delete sequences;
    delete conservationColorer;
    RemovePSSM();
}

BlockMultipleAlignment * BlockMultipleAlignment::Clone(void) const
{
    // must actually copy the list
    BlockMultipleAlignment *copy = new BlockMultipleAlignment(new SequenceList(*sequences), alignmentManager);
    BlockList::const_iterator b, be = blocks.end();
    for (b=blocks.begin(); b!=be; ++b) {
        Block *newBlock = (*b)->Clone(copy);
        copy->blocks.push_back(newBlock);
        if ((*b)->IsAligned()) {
            MarkBlockMap::const_iterator r = markBlocks.find(*b);
            if (r != markBlocks.end())
                copy->markBlocks[newBlock] = true;
        }
    }
    copy->UpdateBlockMapAndColors();
    copy->rowDoubles = rowDoubles;
    copy->rowStrings = rowStrings;
    copy->geometryViolations = geometryViolations;
    copy->showGeometryViolations = showGeometryViolations;
    copy->updateOrigin = updateOrigin;
    copy->alignMasterFrom = alignMasterFrom;
    copy->alignMasterTo = alignMasterTo;
    copy->alignSlaveFrom = alignSlaveFrom;
    copy->alignSlaveTo = alignSlaveTo;
    return copy;
}

static Int4 Round(double d)
{
    if (d >= 0.0)
        return (Int4) (d + 0.5);
    else
        return (Int4) (d - 0.5);
}

const BLAST_Matrix * BlockMultipleAlignment::GetPSSM(void) const
{
    if (pssm) return pssm;

    // for now, use threader's SeqMtf
    BLAST_KarlinBlkPtr karlinBlock = BlastKarlinBlkCreate();
    Seq_Mtf *seqMtf = Threader::CreateSeqMtf(this, 1.0, karlinBlock);

    pssm = (BLAST_Matrix *) MemNew(sizeof(BLAST_Matrix));
    pssm->is_prot = TRUE;
    pssm->name = StringSave("BLOSUM62");
    pssm->karlinK = karlinBlock->K;
    pssm->rows = seqMtf->n + 1;
    pssm->columns = 26;

#ifdef PRINT_PSSM
    FILE *f = fopen("blast_matrix.txt", "w");
#endif

    int i, j;
    pssm->matrix = (Int4 **) MemNew(pssm->rows * sizeof(Int4 *));
    for (i=0; i<pssm->rows; ++i) {
        pssm->matrix[i] = (Int4 *) MemNew(pssm->columns * sizeof(Int4));
#ifdef PRINT_PSSM
        fprintf(f, "matrix %i : ", i);
#endif

        // set scores from threader matrix
        if (i < seqMtf->n) {
            // initialize all rows with custom score, or BLAST_SCORE_MIN; to match what Aron's function creates
            for (j=0; j<pssm->columns; ++j)
                pssm->matrix[i][j] = (j == 21 ? -1 : (j == 25 ? -4 : BLAST_SCORE_MIN));

            for (j=0; j<seqMtf->AlphabetSize; ++j) {
                pssm->matrix[i][LookupBLASTResidueNumberFromThreaderResidueNumber(j)] =
                    Round(((double) seqMtf->ww[i][j]) / Threader::SCALING_FACTOR);
            }
        } else {
            // initialize last row with BLAST_SCORE_MIN
            for (j=0; j<pssm->columns; ++j)
                pssm->matrix[i][j] = BLAST_SCORE_MIN;
        }
#ifdef PRINT_PSSM
        for (j=0; j<pssm->columns; ++j)
            fprintf(f, "%i ", pssm->matrix[i][j]);
        fprintf(f, "\n");
#endif
    }

#ifdef PRINT_PSSM
    // for diffing with scoremat stored in ascii CD
    fprintf(f, "{\n");
    for (i=0; i<seqMtf->n; ++i) {
        for (j=0; j<pssm->columns; ++j) {
            fprintf(f, "      %i,\n", pssm->matrix[i][j]);
        }
    }
    fprintf(f, "}\n");
    // for diffing with .mtx file
    for (i=0; i<seqMtf->n; ++i) {
        for (j=0; j<pssm->columns; ++j) {
            fprintf(f, "%i  ", pssm->matrix[i][j]);
        }
        fprintf(f, "\n");
    }
#endif

#ifdef _DEBUG
    pssm->posFreqs = (Nlm_FloatHi **) MemNew(pssm->rows * sizeof(Nlm_FloatHi *));
    for (i=0; i<pssm->rows; ++i) {
        pssm->posFreqs[i] = (Nlm_FloatHi *) MemNew(pssm->columns * sizeof(Nlm_FloatHi));
#ifdef PRINT_PSSM
        fprintf(f, "freqs %i : ", i);
#endif
        if (i < seqMtf->n) {
            for (j=0; j<seqMtf->AlphabetSize; ++j) {
                pssm->posFreqs[i][LookupBLASTResidueNumberFromThreaderResidueNumber(j)] =
                    seqMtf->freqs[i][j] / Threader::SCALING_FACTOR;
            }
        }
#ifdef PRINT_PSSM
        for (j=0; j<pssm->columns; ++j)
            fprintf(f, "%g ", pssm->posFreqs[i][j]);
        fprintf(f, "\n");
#endif
    }

    // we're not actually using the frequency table; just printing it out for debugging/testing
    for (i=0; i<pssm->rows; ++i) MemFree(pssm->posFreqs[i]);
    MemFree(pssm->posFreqs);
#endif // _DEBUG

    pssm->posFreqs = NULL;

#ifdef PRINT_PSSM
    fclose(f);
#endif

    FreeSeqMtf(seqMtf);
    BlastKarlinBlkDestruct(karlinBlock);
    return pssm;
}

void BlockMultipleAlignment::RemovePSSM(void) const
{
    if (pssm) {
        BLAST_MatrixDestruct(pssm);
        pssm = NULL;
    }
}

void BlockMultipleAlignment::FreeColors(void)
{
    conservationColorer->FreeColors();
    RemovePSSM();
}

bool BlockMultipleAlignment::CheckAlignedBlock(const Block *block) const
{
    if (!block || !block->IsAligned()) {
        ERRORMSG("MultipleAlignment::CheckAlignedBlock() - checks aligned blocks only");
        return false;
    }
    if (block->NSequences() != sequences->size()) {
        ERRORMSG("MultipleAlignment::CheckAlignedBlock() - block size mismatch");
        return false;
    }

    // make sure ranges are reasonable for each sequence
    int row;
    const Block
        *prevBlock = GetBlockBefore(block),
        *nextBlock = GetBlockAfter(block);
    const Block::Range *range, *prevRange = NULL, *nextRange = NULL;
    SequenceList::const_iterator sequence = sequences->begin();
    for (row=0; row<block->NSequences(); ++row, ++sequence) {
        range = block->GetRangeOfRow(row);
        if (prevBlock) prevRange = prevBlock->GetRangeOfRow(row);
        if (nextBlock) nextRange = nextBlock->GetRangeOfRow(row);
        if (range->to - range->from + 1 != block->width ||       // check block width
            (prevRange && range->from <= prevRange->to) ||          // check for range overlap
            (nextRange && range->to >= nextRange->from) ||          // check for range overlap
            range->from > range->to ||                              // check range values
            range->to >= (*sequence)->Length()) {                   // check bounds of end
            ERRORMSG("MultipleAlignment::CheckAlignedBlock() - range error");
            return false;
        }
    }

    return true;
}

bool BlockMultipleAlignment::AddAlignedBlockAtEnd(UngappedAlignedBlock *newBlock)
{
    blocks.push_back(newBlock);
    return CheckAlignedBlock(newBlock);
}

UnalignedBlock * BlockMultipleAlignment::
    CreateNewUnalignedBlockBetween(const Block *leftBlock, const Block *rightBlock)
{
    if ((leftBlock && !leftBlock->IsAligned()) ||
        (rightBlock && !rightBlock->IsAligned())) {
        ERRORMSG("CreateNewUnalignedBlockBetween() - passed an unaligned block");
        return NULL;
    }

    int row, from, to, length;
    SequenceList::const_iterator s, se = sequences->end();

    UnalignedBlock *newBlock = new UnalignedBlock(this);
    newBlock->width = 0;

    for (row=0, s=sequences->begin(); s!=se; ++row, ++s) {

        if (leftBlock)
            from = leftBlock->GetRangeOfRow(row)->to + 1;
        else
            from = 0;

        if (rightBlock)
            to = rightBlock->GetRangeOfRow(row)->from - 1;
        else
            to = (*s)->Length() - 1;

        newBlock->SetRangeOfRow(row, from, to);

        length = to - from + 1;
        if (length < 0) { // just to make sure...
            ERRORMSG("CreateNewUnalignedBlockBetween() - unaligned length < 0");
            return NULL;
        }
        if (length > newBlock->width) newBlock->width = length;
    }

    if (newBlock->width == 0) {
        delete newBlock;
        return NULL;
    } else
        return newBlock;
}

bool BlockMultipleAlignment::AddUnalignedBlocks(void)
{
    BlockList::iterator a, ae = blocks.end();
    const Block *alignedBlock = NULL, *prevAlignedBlock = NULL;
    Block *newUnalignedBlock;

    // unaligned blocks to the left of each aligned block
    for (a=blocks.begin(); a!=ae; ++a) {
        alignedBlock = *a;
        newUnalignedBlock = CreateNewUnalignedBlockBetween(prevAlignedBlock, alignedBlock);