block_multiple_alignment.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:71k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: block_multiple_alignment.cpp,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:27:49  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.57
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: block_multiple_alignment.cpp,v 1000.3 2004/06/01 18:27:49 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE
  13. *               National Center for Biotechnology Information
  14. *
  15. *  This software/database is a "United States Government Work" under the
  16. *  terms of the United States Copyright Act.  It was written as part of
  17. *  the author's official duties as a United States Government employee and
  18. *  thus cannot be copyrighted.  This software/database is freely available
  19. *  to the public for use. The National Library of Medicine and the U.S.
  20. *  Government have not placed any restriction on its use or reproduction.
  21. *
  22. *  Although all reasonable efforts have been taken to ensure the accuracy
  23. *  and reliability of the software and data, the NLM and the U.S.
  24. *  Government do not and cannot warrant the performance or results that
  25. *  may be obtained by using this software or data. The NLM and the U.S.
  26. *  Government disclaim all warranties, express or implied, including
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.
  29. *
  30. *  Please cite the author in any work or product based on this material.
  31. *
  32. * ===========================================================================
  33. *
  34. * Authors:  Paul Thiessen
  35. *
  36. * File Description:
  37. *      Classes to hold alignment data
  38. *
  39. * ===========================================================================
  40. */
  41. #ifdef _MSC_VER
  42. #pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
  43. #endif
  44. #include <ncbi_pch.hpp>
  45. #include <corelib/ncbistd.hpp>
  46. #include <objects/seqalign/Dense_diag.hpp>
  47. #include "block_multiple_alignment.hpp"
  48. #include "sequence_set.hpp"
  49. #include "molecule.hpp"
  50. #include "conservation_colorer.hpp"
  51. #include "style_manager.hpp"
  52. #include "structure_set.hpp"
  53. #include "messenger.hpp"
  54. #include "cn3d_colors.hpp"
  55. #include "alignment_manager.hpp"
  56. #include "cn3d_tools.hpp"
  57. #include "molecule_identifier.hpp"
  58. #include "cn3d_threader.hpp"
  59. #include "cn3d_blast.hpp"
  60. // hack so I can catch memory leaks specific to this module, at the line where allocation occurs
  61. #ifdef _DEBUG
  62. #ifdef MemNew
  63. #undef MemNew
  64. #endif
  65. #define MemNew(sz) memset(malloc(sz), 0, (sz))
  66. #endif
  67. USING_NCBI_SCOPE;
  68. USING_SCOPE(objects);
  69. BEGIN_SCOPE(Cn3D)
  70. BlockMultipleAlignment::BlockMultipleAlignment(SequenceList *sequenceList, AlignmentManager *alnMgr) :
  71.     sequences(sequenceList), conservationColorer(NULL), alignmentManager(alnMgr),
  72.     alignMasterFrom(-1), alignMasterTo(-1), alignSlaveFrom(-1), alignSlaveTo(-1),
  73.     pssm(NULL), showGeometryViolations(false)
  74. {
  75.     InitCache();
  76.     rowDoubles.resize(sequenceList->size(), 0.0);
  77.     rowStrings.resize(sequenceList->size());
  78.     // create conservation colorer
  79.     conservationColorer = new ConservationColorer(this);
  80. }
  81. void BlockMultipleAlignment::InitCache(void)
  82. {
  83.     cachePrevRow = -1;
  84.     cachePrevBlock = NULL;
  85.     cacheBlockIterator = blocks.begin();
  86. }
  87. BlockMultipleAlignment::~BlockMultipleAlignment(void)
  88. {
  89.     BlockList::iterator i, ie = blocks.end();
  90.     for (i=blocks.begin(); i!=ie; ++i) delete *i;
  91.     delete sequences;
  92.     delete conservationColorer;
  93.     RemovePSSM();
  94. }
  95. BlockMultipleAlignment * BlockMultipleAlignment::Clone(void) const
  96. {
  97.     // must actually copy the list
  98.     BlockMultipleAlignment *copy = new BlockMultipleAlignment(new SequenceList(*sequences), alignmentManager);
  99.     BlockList::const_iterator b, be = blocks.end();
  100.     for (b=blocks.begin(); b!=be; ++b) {
  101.         Block *newBlock = (*b)->Clone(copy);
  102.         copy->blocks.push_back(newBlock);
  103.         if ((*b)->IsAligned()) {
  104.             MarkBlockMap::const_iterator r = markBlocks.find(*b);
  105.             if (r != markBlocks.end())
  106.                 copy->markBlocks[newBlock] = true;
  107.         }
  108.     }
  109.     copy->UpdateBlockMapAndColors();
  110.     copy->rowDoubles = rowDoubles;
  111.     copy->rowStrings = rowStrings;
  112.     copy->geometryViolations = geometryViolations;
  113.     copy->showGeometryViolations = showGeometryViolations;
  114.     copy->updateOrigin = updateOrigin;
  115.     copy->alignMasterFrom = alignMasterFrom;
  116.     copy->alignMasterTo = alignMasterTo;
  117.     copy->alignSlaveFrom = alignSlaveFrom;
  118.     copy->alignSlaveTo = alignSlaveTo;
  119.     return copy;
  120. }
  121. static Int4 Round(double d)
  122. {
  123.     if (d >= 0.0)
  124.         return (Int4) (d + 0.5);
  125.     else
  126.         return (Int4) (d - 0.5);
  127. }
  128. const BLAST_Matrix * BlockMultipleAlignment::GetPSSM(void) const
  129. {
  130.     if (pssm) return pssm;
  131.     // for now, use threader's SeqMtf
  132.     BLAST_KarlinBlkPtr karlinBlock = BlastKarlinBlkCreate();
  133.     Seq_Mtf *seqMtf = Threader::CreateSeqMtf(this, 1.0, karlinBlock);
  134.     pssm = (BLAST_Matrix *) MemNew(sizeof(BLAST_Matrix));
  135.     pssm->is_prot = TRUE;
  136.     pssm->name = StringSave("BLOSUM62");
  137.     pssm->karlinK = karlinBlock->K;
  138.     pssm->rows = seqMtf->n + 1;
  139.     pssm->columns = 26;
  140. #ifdef PRINT_PSSM
  141.     FILE *f = fopen("blast_matrix.txt", "w");
  142. #endif
  143.     int i, j;
  144.     pssm->matrix = (Int4 **) MemNew(pssm->rows * sizeof(Int4 *));
  145.     for (i=0; i<pssm->rows; ++i) {
  146.         pssm->matrix[i] = (Int4 *) MemNew(pssm->columns * sizeof(Int4));
  147. #ifdef PRINT_PSSM
  148.         fprintf(f, "matrix %i : ", i);
  149. #endif
  150.         // set scores from threader matrix
  151.         if (i < seqMtf->n) {
  152.             // initialize all rows with custom score, or BLAST_SCORE_MIN; to match what Aron's function creates
  153.             for (j=0; j<pssm->columns; ++j)
  154.                 pssm->matrix[i][j] = (j == 21 ? -1 : (j == 25 ? -4 : BLAST_SCORE_MIN));
  155.             for (j=0; j<seqMtf->AlphabetSize; ++j) {
  156.                 pssm->matrix[i][LookupBLASTResidueNumberFromThreaderResidueNumber(j)] =
  157.                     Round(((double) seqMtf->ww[i][j]) / Threader::SCALING_FACTOR);
  158.             }
  159.         } else {
  160.             // initialize last row with BLAST_SCORE_MIN
  161.             for (j=0; j<pssm->columns; ++j)
  162.                 pssm->matrix[i][j] = BLAST_SCORE_MIN;
  163.         }
  164. #ifdef PRINT_PSSM
  165.         for (j=0; j<pssm->columns; ++j)
  166.             fprintf(f, "%i ", pssm->matrix[i][j]);
  167.         fprintf(f, "n");
  168. #endif
  169.     }
  170. #ifdef PRINT_PSSM
  171.     // for diffing with scoremat stored in ascii CD
  172.     fprintf(f, "{n");
  173.     for (i=0; i<seqMtf->n; ++i) {
  174.         for (j=0; j<pssm->columns; ++j) {
  175.             fprintf(f, "      %i,n", pssm->matrix[i][j]);
  176.         }
  177.     }
  178.     fprintf(f, "}n");
  179.     // for diffing with .mtx file
  180.     for (i=0; i<seqMtf->n; ++i) {
  181.         for (j=0; j<pssm->columns; ++j) {
  182.             fprintf(f, "%i  ", pssm->matrix[i][j]);
  183.         }
  184.         fprintf(f, "n");
  185.     }
  186. #endif
  187. #ifdef _DEBUG
  188.     pssm->posFreqs = (Nlm_FloatHi **) MemNew(pssm->rows * sizeof(Nlm_FloatHi *));
  189.     for (i=0; i<pssm->rows; ++i) {
  190.         pssm->posFreqs[i] = (Nlm_FloatHi *) MemNew(pssm->columns * sizeof(Nlm_FloatHi));
  191. #ifdef PRINT_PSSM
  192.         fprintf(f, "freqs %i : ", i);
  193. #endif
  194.         if (i < seqMtf->n) {
  195.             for (j=0; j<seqMtf->AlphabetSize; ++j) {
  196.                 pssm->posFreqs[i][LookupBLASTResidueNumberFromThreaderResidueNumber(j)] =
  197.                     seqMtf->freqs[i][j] / Threader::SCALING_FACTOR;
  198.             }
  199.         }
  200. #ifdef PRINT_PSSM
  201.         for (j=0; j<pssm->columns; ++j)
  202.             fprintf(f, "%g ", pssm->posFreqs[i][j]);
  203.         fprintf(f, "n");
  204. #endif
  205.     }
  206.     // we're not actually using the frequency table; just printing it out for debugging/testing
  207.     for (i=0; i<pssm->rows; ++i) MemFree(pssm->posFreqs[i]);
  208.     MemFree(pssm->posFreqs);
  209. #endif // _DEBUG
  210.     pssm->posFreqs = NULL;
  211. #ifdef PRINT_PSSM
  212.     fclose(f);
  213. #endif
  214.     FreeSeqMtf(seqMtf);
  215.     BlastKarlinBlkDestruct(karlinBlock);
  216.     return pssm;
  217. }
  218. void BlockMultipleAlignment::RemovePSSM(void) const
  219. {
  220.     if (pssm) {
  221.         BLAST_MatrixDestruct(pssm);
  222.         pssm = NULL;
  223.     }
  224. }
  225. void BlockMultipleAlignment::FreeColors(void)
  226. {
  227.     conservationColorer->FreeColors();
  228.     RemovePSSM();
  229. }
  230. bool BlockMultipleAlignment::CheckAlignedBlock(const Block *block) const
  231. {
  232.     if (!block || !block->IsAligned()) {
  233.         ERRORMSG("MultipleAlignment::CheckAlignedBlock() - checks aligned blocks only");
  234.         return false;
  235.     }
  236.     if (block->NSequences() != sequences->size()) {
  237.         ERRORMSG("MultipleAlignment::CheckAlignedBlock() - block size mismatch");
  238.         return false;
  239.     }
  240.     // make sure ranges are reasonable for each sequence
  241.     int row;
  242.     const Block
  243.         *prevBlock = GetBlockBefore(block),
  244.         *nextBlock = GetBlockAfter(block);
  245.     const Block::Range *range, *prevRange = NULL, *nextRange = NULL;
  246.     SequenceList::const_iterator sequence = sequences->begin();
  247.     for (row=0; row<block->NSequences(); ++row, ++sequence) {
  248.         range = block->GetRangeOfRow(row);
  249.         if (prevBlock) prevRange = prevBlock->GetRangeOfRow(row);
  250.         if (nextBlock) nextRange = nextBlock->GetRangeOfRow(row);
  251.         if (range->to - range->from + 1 != block->width ||       // check block width
  252.             (prevRange && range->from <= prevRange->to) ||          // check for range overlap
  253.             (nextRange && range->to >= nextRange->from) ||          // check for range overlap
  254.             range->from > range->to ||                              // check range values
  255.             range->to >= (*sequence)->Length()) {                   // check bounds of end
  256.             ERRORMSG("MultipleAlignment::CheckAlignedBlock() - range error");
  257.             return false;
  258.         }
  259.     }
  260.     return true;
  261. }
  262. bool BlockMultipleAlignment::AddAlignedBlockAtEnd(UngappedAlignedBlock *newBlock)
  263. {
  264.     blocks.push_back(newBlock);
  265.     return CheckAlignedBlock(newBlock);
  266. }
  267. UnalignedBlock * BlockMultipleAlignment::
  268.     CreateNewUnalignedBlockBetween(const Block *leftBlock, const Block *rightBlock)
  269. {
  270.     if ((leftBlock && !leftBlock->IsAligned()) ||
  271.         (rightBlock && !rightBlock->IsAligned())) {
  272.         ERRORMSG("CreateNewUnalignedBlockBetween() - passed an unaligned block");
  273.         return NULL;
  274.     }
  275.     int row, from, to, length;
  276.     SequenceList::const_iterator s, se = sequences->end();
  277.     UnalignedBlock *newBlock = new UnalignedBlock(this);
  278.     newBlock->width = 0;
  279.     for (row=0, s=sequences->begin(); s!=se; ++row, ++s) {
  280.         if (leftBlock)
  281.             from = leftBlock->GetRangeOfRow(row)->to + 1;
  282.         else
  283.             from = 0;
  284.         if (rightBlock)
  285.             to = rightBlock->GetRangeOfRow(row)->from - 1;
  286.         else
  287.             to = (*s)->Length() - 1;
  288.         newBlock->SetRangeOfRow(row, from, to);
  289.         length = to - from + 1;
  290.         if (length < 0) { // just to make sure...
  291.             ERRORMSG("CreateNewUnalignedBlockBetween() - unaligned length < 0");
  292.             return NULL;
  293.         }
  294.         if (length > newBlock->width) newBlock->width = length;
  295.     }
  296.     if (newBlock->width == 0) {
  297.         delete newBlock;
  298.         return NULL;
  299.     } else
  300.         return newBlock;
  301. }
  302. bool BlockMultipleAlignment::AddUnalignedBlocks(void)
  303. {
  304.     BlockList::iterator a, ae = blocks.end();
  305.     const Block *alignedBlock = NULL, *prevAlignedBlock = NULL;
  306.     Block *newUnalignedBlock;
  307.     // unaligned blocks to the left of each aligned block
  308.     for (a=blocks.begin(); a!=ae; ++a) {
  309.         alignedBlock = *a;
  310.         newUnalignedBlock = CreateNewUnalignedBlockBetween(prevAlignedBlock, alignedBlock);
  311.         if (newUnalignedBlock) blocks.insert(a, newUnalignedBlock);
  312.         prevAlignedBlock = alignedBlock;
  313.     }
  314.     // right tail
  315.     newUnalignedBlock = CreateNewUnalignedBlockBetween(alignedBlock, NULL);
  316.     if (newUnalignedBlock) {
  317.         blocks.insert(a, newUnalignedBlock);
  318.     }
  319.     return true;
  320. }
  321. bool BlockMultipleAlignment::UpdateBlockMapAndColors(bool clearRowInfo)
  322. {
  323.     int i = 0, j, n = 0;
  324.     BlockList::iterator b, be = blocks.end();
  325.     // reset old stuff, recalculate width
  326.     totalWidth = 0;
  327.     for (b=blocks.begin(); b!=be; ++b) totalWidth += (*b)->width;
  328. //    TESTMSG("alignment display size: " << totalWidth << " x " << NRows());
  329.     // fill out the block map
  330.     conservationColorer->Clear();
  331.     blockMap.resize(totalWidth);
  332.     UngappedAlignedBlock *aBlock;
  333.     for (b=blocks.begin(); b!=be; ++b) {
  334.         aBlock = dynamic_cast<UngappedAlignedBlock*>(*b);
  335.         if (aBlock) {
  336.             conservationColorer->AddBlock(aBlock);
  337.             ++n;
  338.         }
  339.         for (j=0; j<(*b)->width; ++j, ++i) {
  340.             blockMap[i].block = *b;
  341.             blockMap[i].blockColumn = j;
  342.             blockMap[i].alignedBlockNum = aBlock ? n : -1;
  343.         }
  344.     }
  345.     // if alignment changes, any pssm/scores/status/special colors become invalid
  346.     RemovePSSM();
  347.     if (clearRowInfo) ClearRowInfo();
  348.     ShowGeometryViolations(showGeometryViolations); // recalculate GV's
  349.     return true;
  350. }
  351. bool BlockMultipleAlignment::GetCharacterTraitsAt(
  352.     int alignmentColumn, int row, eUnalignedJustification justification,
  353.     char *character, Vector *color, bool *isHighlighted,
  354.     bool *drawBackground, Vector *cellBackgroundColor) const
  355. {
  356.     const Sequence *sequence;
  357.     int seqIndex;
  358.     bool isAligned;
  359.     if (!GetSequenceAndIndexAt(alignmentColumn, row, justification, &sequence, &seqIndex, &isAligned))
  360.         return false;
  361.     *character = (seqIndex >= 0) ? sequence->sequenceString[seqIndex] : '~';
  362.     if (isAligned)
  363.         *character = toupper(*character);
  364.     else
  365.         *character = tolower(*character);
  366.     // try to color by molecule first
  367.     if (sequence->molecule) {
  368.         *color = (seqIndex >= 0) ?
  369.             sequence->molecule->GetResidueColor(seqIndex) :
  370.             GlobalColors()->Get(Colors::eNoCoordinates);;
  371.     }
  372.     // otherwise (for unstructured sequence):
  373.     else {
  374.         StyleSettings::eColorScheme colorScheme =
  375.             sequence->parentSet->styleManager->GetGlobalStyle().proteinBackbone.colorScheme;
  376.         // color by hydrophobicity
  377.         if (sequence->isProtein && colorScheme == StyleSettings::eHydrophobicity) {
  378.             double hydrophobicity = GetHydrophobicity(toupper(*character));
  379.             *color = (hydrophobicity != UNKNOWN_HYDROPHOBICITY) ?
  380.                 GlobalColors()->Get(Colors::eHydrophobicityMap, hydrophobicity) :
  381.                 GlobalColors()->Get(Colors::eNoHydrophobicity);
  382.         }
  383.         // or color by charge
  384.         else if (sequence->isProtein && colorScheme == StyleSettings::eCharge) {
  385.             int charge = GetCharge(toupper(*character));
  386.             *color = GlobalColors()->Get(
  387.                 (charge > 0) ? Colors::ePositive : ((charge < 0) ? Colors::eNegative : Colors::eNeutral));;
  388.         }
  389.         // else, color by alignment color
  390.         else {
  391.             const Vector *aColor;
  392.             if (isAligned && (aColor = GetAlignmentColor(row, seqIndex, colorScheme)) != NULL) {
  393.                 *color = *aColor;
  394.             } else {
  395.                 *color = GlobalColors()->Get(Colors::eUnaligned);
  396.             }
  397.         }
  398.     }
  399.     if (seqIndex >= 0)
  400.         *isHighlighted = GlobalMessenger()->IsHighlighted(sequence, seqIndex);
  401.     else
  402.         *isHighlighted = false;
  403.     // add special alignment coloring (but don't override highlight)
  404.     *drawBackground = false;
  405.     if (!(*isHighlighted)) {
  406.         // check for block flagged for realignment
  407.         if (markBlocks.find(blockMap[alignmentColumn].block) != markBlocks.end()) {
  408.             *drawBackground = true;
  409.             *cellBackgroundColor = GlobalColors()->Get(Colors::eMarkBlock);
  410.         }
  411.         // optionally show geometry violations
  412.         if (showGeometryViolations && seqIndex >= 0 && geometryViolations[row][seqIndex]) {
  413.             *drawBackground = true;
  414.             *cellBackgroundColor = GlobalColors()->Get(Colors::eGeometryViolation);
  415.         }
  416.         // check for unmergeable alignment
  417.         const BlockMultipleAlignment *referenceAlignment = alignmentManager->GetCurrentMultipleAlignment();
  418.         if (referenceAlignment && referenceAlignment != this &&
  419. seqIndex >= 0 && GetMaster() == referenceAlignment->GetMaster()) {
  420.             bool unmergeable = false;
  421.             const Block *block = blockMap[alignmentColumn].block;
  422.             // case where master residues are aligned in multiple but not in this one
  423.             if (row == 0 && !isAligned && referenceAlignment->IsAligned(row, seqIndex))
  424.                 unmergeable = true;
  425.             // block boundaries in this inside aligned block of multiple
  426.             else if (row > 0 && isAligned) {
  427.                 const Block::Range *range = block->GetRangeOfRow(row);
  428.                 bool
  429.                     isLeftEdge = (seqIndex == range->from),
  430.                     isRightEdge = (seqIndex == range->to);
  431.                 if (isLeftEdge || isRightEdge) {
  432.                     // get corresponding block of multiple
  433.                     const Block::Range *masterRange = block->GetRangeOfRow(0);
  434.                     int masterIndex = masterRange->from + seqIndex - range->from;
  435.                     const Block *multipleBlock = referenceAlignment->GetBlock(0, masterIndex);
  436.                     masterRange = multipleBlock->GetRangeOfRow(0);
  437.                     if (multipleBlock->IsAligned() &&
  438.                         ((isLeftEdge && masterIndex > masterRange->from) ||
  439.                          (isRightEdge && masterIndex < masterRange->to)))
  440.                         unmergeable = true;
  441.                 }
  442.             }
  443.             // check for inserts relative to the multiple
  444.             else if (row > 0 && !isAligned) {
  445.                 const Block
  446.                     *blockBefore = GetBlockBefore(block),
  447.                     *blockAfter = GetBlockAfter(block),
  448.                     *refBlock;
  449.                 if (blockBefore && blockBefore->IsAligned() && blockAfter && blockAfter->IsAligned() &&
  450.                         referenceAlignment->GetBlock(0, blockBefore->GetRangeOfRow(0)->to) ==
  451.                         (refBlock=referenceAlignment->GetBlock(0, blockAfter->GetRangeOfRow(0)->from)) &&
  452.                         refBlock->IsAligned()
  453.                     )
  454.                     unmergeable = true;
  455.             }
  456.             if (unmergeable) {
  457.                 *drawBackground = true;
  458.                 *cellBackgroundColor = GlobalColors()->Get(Colors::eMergeFail);
  459.             }
  460.         }
  461.     }
  462.     return true;
  463. }
  464. bool BlockMultipleAlignment::GetSequenceAndIndexAt(
  465.     int alignmentColumn, int row, eUnalignedJustification requestedJustification,
  466.     const Sequence **sequence, int *index, bool *isAligned) const
  467. {
  468.     if (sequence) *sequence = (*sequences)[row];
  469.     const BlockInfo& blockInfo = blockMap[alignmentColumn];
  470.     if (!blockInfo.block->IsAligned()) {
  471.         if (isAligned) *isAligned = false;
  472.         // override requested justification for end blocks
  473.         if (blockInfo.block == blocks.back()) // also true if there's a single aligned block
  474.             requestedJustification = eLeft;
  475.         else if (blockInfo.block == blocks.front())
  476.             requestedJustification = eRight;
  477.     } else
  478.         if (isAligned) *isAligned = true;
  479.     if (index)
  480.         *index = blockInfo.block->GetIndexAt(blockInfo.blockColumn, row, requestedJustification);
  481.     return true;
  482. }
  483. int BlockMultipleAlignment::GetRowForSequence(const Sequence *sequence) const
  484. {
  485.     // this only works for structured sequences, since non-structure sequences can
  486.     // be repeated any number of times in the alignment; assumes repeated structures
  487.     // will each have a unique Sequence object
  488.     if (!sequence || !sequence->molecule) {
  489.         ERRORMSG("BlockMultipleAlignment::GetRowForSequence() - Sequence must have associated structure");
  490.         return -1;
  491.     }
  492.     if (cachePrevRow < 0 || cachePrevRow >= NRows() || sequence != (*sequences)[cachePrevRow]) {
  493.         int row;
  494.         for (row=0; row<NRows(); ++row) if ((*sequences)[row] == sequence) break;
  495.         if (row == NRows()) {
  496. //            ERRORMSG("BlockMultipleAlignment::GetRowForSequence() - can't find given Sequence");
  497.             return -1;
  498.         }
  499.         cachePrevRow = row;
  500.     }
  501.     return cachePrevRow;
  502. }
  503. const Vector * BlockMultipleAlignment::GetAlignmentColor(
  504.     int row, int seqIndex, StyleSettings::eColorScheme colorScheme) const
  505. {
  506.      const UngappedAlignedBlock *block = dynamic_cast<const UngappedAlignedBlock*>(GetBlock(row, seqIndex));
  507.      if (!block || !block->IsAligned()) {
  508.         WARNINGMSG("BlockMultipleAlignment::GetAlignmentColor() called on unaligned residue");
  509.         return NULL;
  510.     }
  511.     const Vector *alignedColor;
  512.     switch (colorScheme) {
  513.         case StyleSettings::eAligned:
  514.             alignedColor = GlobalColors()->Get(Colors::eConservationMap, Colors::nConservationMap - 1);
  515.             break;
  516.         case StyleSettings::eIdentity:
  517.             alignedColor = conservationColorer->
  518.                 GetIdentityColor(block, seqIndex - block->GetRangeOfRow(row)->from);
  519.             break;
  520.         case StyleSettings::eVariety:
  521.             alignedColor = conservationColorer->
  522.                 GetVarietyColor(block, seqIndex - block->GetRangeOfRow(row)->from);
  523.             break;
  524.         case StyleSettings::eWeightedVariety:
  525.             alignedColor = conservationColorer->
  526.                 GetWeightedVarietyColor(block, seqIndex - block->GetRangeOfRow(row)->from);
  527.             break;
  528.         case StyleSettings::eInformationContent:
  529.             alignedColor = conservationColorer->
  530.                 GetInformationContentColor(block, seqIndex - block->GetRangeOfRow(row)->from);
  531.             break;
  532.         case StyleSettings::eFit:
  533.             alignedColor = conservationColorer->
  534.                 GetFitColor(block, seqIndex - block->GetRangeOfRow(row)->from, row);
  535.             break;
  536.         case StyleSettings::eBlockFit:
  537.             alignedColor = conservationColorer->GetBlockFitColor(block, row);
  538.             break;
  539.         case StyleSettings::eBlockZFit:
  540.             alignedColor = conservationColorer->GetBlockZFitColor(block, row);
  541.             break;
  542.         case StyleSettings::eBlockRowFit:
  543.             alignedColor = conservationColorer->GetBlockRowFitColor(block, row);
  544.             break;
  545.         default:
  546.             alignedColor = NULL;
  547.     }
  548.     return alignedColor;
  549. }
  550. const Vector * BlockMultipleAlignment::GetAlignmentColor(
  551.     const Sequence *sequence, int seqIndex, StyleSettings::eColorScheme colorScheme) const
  552. {
  553.     int row = GetRowForSequence(sequence);
  554.     if (row < 0) return NULL;
  555.     return GetAlignmentColor(row, seqIndex, colorScheme);
  556. }
  557. bool BlockMultipleAlignment::IsAligned(int row, int seqIndex) const
  558. {
  559.     const Block *block = GetBlock(row, seqIndex);
  560.     return (block && block->IsAligned());
  561. }
  562. const Block * BlockMultipleAlignment::GetBlock(int row, int seqIndex) const
  563. {
  564.     // make sure we're in range for this sequence
  565.     if (row < 0 || seqIndex < 0 || row >= NRows() || seqIndex >= (*sequences)[row]->Length()) {
  566.         ERRORMSG("BlockMultipleAlignment::GetBlock() - coordinate out of range");
  567.         return NULL;
  568.     }
  569.     const Block::Range *range;
  570.     // first check to see if it's in the same block as last time.
  571.     if (cachePrevBlock) {
  572.         range = cachePrevBlock->GetRangeOfRow(row);
  573.         if (seqIndex >= range->from && seqIndex <= range->to) return cachePrevBlock;
  574.         ++cacheBlockIterator; // start search at next block
  575.     } else {
  576.         cacheBlockIterator = blocks.begin();
  577.     }
  578.     // otherwise, perform block search. This search is most efficient when queries
  579.     // happen in order from left to right along a given row.
  580.     do {
  581.         if (cacheBlockIterator == blocks.end()) cacheBlockIterator = blocks.begin();
  582.         range = (*cacheBlockIterator)->GetRangeOfRow(row);
  583.         if (seqIndex >= range->from && seqIndex <= range->to) {
  584.             cachePrevBlock = *cacheBlockIterator; // cache this block
  585.             return cachePrevBlock;
  586.         }
  587.         ++cacheBlockIterator;
  588.     } while (1);
  589. }
  590. int BlockMultipleAlignment::GetFirstAlignedBlockPosition(void) const
  591. {
  592.     BlockList::const_iterator b = blocks.begin();
  593.     if (blocks.size() > 0 && (*b)->IsAligned()) // first block is aligned
  594.         return 0;
  595.     else if (blocks.size() >= 2 && (*(++b))->IsAligned()) // second block is aligned
  596.         return blocks.front()->width;
  597.     else
  598.         return -1;
  599. }
  600. int BlockMultipleAlignment::GetAlignedSlaveIndex(int masterSeqIndex, int slaveRow) const
  601. {
  602.     const UngappedAlignedBlock
  603.         *aBlock = dynamic_cast<const UngappedAlignedBlock*>(GetBlock(0, masterSeqIndex));
  604.     if (!aBlock) return -1;
  605.     const Block::Range
  606.         *masterRange = aBlock->GetRangeOfRow(0),
  607.         *slaveRange = aBlock->GetRangeOfRow(slaveRow);
  608.     return (slaveRange->from + masterSeqIndex - masterRange->from);
  609. }
  610. void BlockMultipleAlignment::SelectedRange(int row, int from, int to,
  611.     eUnalignedJustification justification, bool toggle) const
  612. {
  613.     // translate from,to (alignment columns) into sequence indexes
  614.     const Sequence *sequence;
  615.     int fromIndex, toIndex;
  616.     bool ignored;
  617.     // trim selection area to size of this alignment
  618.     if (to >= totalWidth) to = totalWidth - 1;
  619.     // find first residue within range
  620.     while (from <= to) {
  621.         GetSequenceAndIndexAt(from, row, justification, &sequence, &fromIndex, &ignored);
  622.         if (fromIndex >= 0) break;
  623.         ++from;
  624.     }
  625.     if (from > to) return;
  626.     // find last residue within range
  627.     while (to >= from) {
  628.         GetSequenceAndIndexAt(to, row, justification, &sequence, &toIndex, &ignored);
  629.         if (toIndex >= 0) break;
  630.         --to;
  631.     }
  632.     if (toggle)
  633.         GlobalMessenger()->ToggleHighlights(sequence, fromIndex, toIndex);
  634.     else
  635.         GlobalMessenger()->AddHighlights(sequence, fromIndex, toIndex);
  636. }
  637. void BlockMultipleAlignment::GetAlignedBlockPosition(int alignmentIndex,
  638.     int *blockColumn, int *blockWidth) const
  639. {
  640.     *blockColumn = *blockWidth = -1;
  641.     const BlockInfo& info = blockMap[alignmentIndex];
  642.     if (info.block->IsAligned()) {
  643.         *blockColumn = info.blockColumn;
  644.         *blockWidth = info.block->width;
  645.     }
  646. }
  647. Block * BlockMultipleAlignment::GetBlockBefore(const Block *block) const
  648. {
  649.     Block *prevBlock = NULL;
  650.     BlockList::const_iterator b, be = blocks.end();
  651.     for (b=blocks.begin(); b!=be; ++b) {
  652.         if (*b == block) break;
  653.         prevBlock = *b;
  654.     }
  655.     return prevBlock;
  656. }
  657. const UnalignedBlock * BlockMultipleAlignment::GetUnalignedBlockBefore(
  658.     const UngappedAlignedBlock *aBlock) const
  659. {
  660.     const Block *prevBlock;
  661.     if (aBlock)
  662.         prevBlock = GetBlockBefore(aBlock);
  663.     else
  664.         prevBlock = blocks.back();
  665.     return dynamic_cast<const UnalignedBlock*>(prevBlock);
  666. }
  667. Block * BlockMultipleAlignment::GetBlockAfter(const Block *block) const
  668. {
  669.     BlockList::const_iterator b, be = blocks.end();
  670.     for (b=blocks.begin(); b!=be; ++b) {
  671.         if (*b == block) {
  672.             ++b;
  673.             if (b == be) break;
  674.             return *b;
  675.         }
  676.     }
  677.     return NULL;
  678. }
  679. void BlockMultipleAlignment::InsertBlockBefore(Block *newBlock, const Block *insertAt)
  680. {
  681.     BlockList::iterator b, be = blocks.end();
  682.     for (b=blocks.begin(); b!=be; ++b) {
  683.         if (*b == insertAt) {
  684.             blocks.insert(b, newBlock);
  685.             return;
  686.         }
  687.     }
  688.     WARNINGMSG("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
  689. }
  690. void BlockMultipleAlignment::InsertBlockAfter(const Block *insertAt, Block *newBlock)
  691. {
  692.     BlockList::iterator b, be = blocks.end();
  693.     for (b=blocks.begin(); b!=be; ++b) {
  694.         if (*b == insertAt) {
  695.             ++b;
  696.             blocks.insert(b, newBlock);
  697.             return;
  698.         }
  699.     }
  700.     WARNINGMSG("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
  701. }
  702. void BlockMultipleAlignment::RemoveBlock(Block *block)
  703. {
  704.     BlockList::iterator b, be = blocks.end();
  705.     for (b=blocks.begin(); b!=be; ++b) {
  706.         if (*b == block) {
  707.             delete *b;
  708.             blocks.erase(b);
  709.             InitCache();
  710.             return;
  711.         }
  712.     }
  713.     WARNINGMSG("BlockMultipleAlignment::RemoveBlock() - couldn't find block");
  714. }
  715. bool BlockMultipleAlignment::MoveBlockBoundary(int columnFrom, int columnTo)
  716. {
  717.     int blockColumn, blockWidth;
  718.     GetAlignedBlockPosition(columnFrom, &blockColumn, &blockWidth);
  719.     if (blockColumn < 0 || blockWidth <= 0) return false;
  720.     TRACEMSG("trying to move block boundary from " << columnFrom << " to " << columnTo);
  721.     const BlockInfo& info = blockMap[columnFrom];
  722.     int row, requestedShift = columnTo - columnFrom, actualShift = 0;
  723.     const Block::Range *range;
  724.     // shrink block from left
  725.     if (blockColumn == 0 && requestedShift > 0 && requestedShift < info.block->width) {
  726.         actualShift = requestedShift;
  727.         TRACEMSG("shrinking block from left");
  728.         for (row=0; row<NRows(); ++row) {
  729.             range = info.block->GetRangeOfRow(row);
  730.             info.block->SetRangeOfRow(row, range->from + requestedShift, range->to);
  731.         }
  732.         info.block->width -= requestedShift;
  733.         Block *prevBlock = GetBlockBefore(info.block);
  734.         if (prevBlock && !prevBlock->IsAligned()) {
  735.             for (row=0; row<NRows(); ++row) {
  736.                 range = prevBlock->GetRangeOfRow(row);
  737.                 prevBlock->SetRangeOfRow(row, range->from, range->to + requestedShift);
  738.             }
  739.             prevBlock->width += requestedShift;
  740.         } else {
  741.             Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(prevBlock, info.block);
  742.             if (newUnalignedBlock) InsertBlockBefore(newUnalignedBlock, info.block);
  743.             TRACEMSG("added new unaligned block");
  744.         }
  745.     }
  746.     // shrink block from right (requestedShift < 0)
  747.     else if (blockColumn == info.block->width - 1 &&
  748.                 requestedShift < 0 && requestedShift > -info.block->width) {
  749.         actualShift = requestedShift;
  750.         TRACEMSG("shrinking block from right");
  751.         for (row=0; row<NRows(); ++row) {
  752.             range = info.block->GetRangeOfRow(row);
  753.             info.block->SetRangeOfRow(row, range->from, range->to + requestedShift);
  754.         }
  755.         info.block->width += requestedShift;
  756.         Block *nextBlock = GetBlockAfter(info.block);
  757.         if (nextBlock && !nextBlock->IsAligned()) {
  758.             for (row=0; row<NRows(); ++row) {
  759.                 range = nextBlock->GetRangeOfRow(row);
  760.                 nextBlock->SetRangeOfRow(row, range->from + requestedShift, range->to);
  761.             }
  762.             nextBlock->width -= requestedShift;
  763.         } else {
  764.             Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(info.block, nextBlock);
  765.             if (newUnalignedBlock) InsertBlockAfter(info.block, newUnalignedBlock);
  766.             TRACEMSG("added new unaligned block");
  767.         }
  768.     }
  769.     // grow block to right
  770.     else if (blockColumn == info.block->width - 1 && requestedShift > 0) {
  771.         Block *nextBlock = GetBlockAfter(info.block);
  772.         if (nextBlock && !nextBlock->IsAligned()) {
  773.             int nRes;
  774.             actualShift = requestedShift;
  775.             for (row=0; row<NRows(); ++row) {
  776.                 range = nextBlock->GetRangeOfRow(row);
  777.                 nRes = range->to - range->from + 1;
  778.                 if (nRes < actualShift) actualShift = nRes;
  779.             }
  780.             if (actualShift) {
  781.                 TRACEMSG("growing block to right");
  782.                 for (row=0; row<NRows(); ++row) {
  783.                     range = info.block->GetRangeOfRow(row);
  784.                     info.block->SetRangeOfRow(row, range->from, range->to + actualShift);
  785.                     range = nextBlock->GetRangeOfRow(row);
  786.                     nextBlock->SetRangeOfRow(row, range->from + actualShift, range->to);
  787.                 }
  788.                 info.block->width += actualShift;
  789.                 nextBlock->width -= actualShift;
  790.                 if (nextBlock->width == 0) {
  791.                     RemoveBlock(nextBlock);
  792.                     TRACEMSG("removed empty block");
  793.                 }
  794.             }
  795.         }
  796.     }
  797.     // grow block to left (requestedShift < 0)
  798.     else if (blockColumn == 0 && requestedShift < 0) {
  799.         Block *prevBlock = GetBlockBefore(info.block);
  800.         if (prevBlock && !prevBlock->IsAligned()) {
  801.             int nRes;
  802.             actualShift = requestedShift;
  803.             for (row=0; row<NRows(); ++row) {
  804.                 range = prevBlock->GetRangeOfRow(row);
  805.                 nRes = range->to - range->from + 1;
  806.                 if (nRes < -actualShift) actualShift = -nRes;
  807.             }
  808.             if (actualShift) {
  809.                 TRACEMSG("growing block to left");
  810.                 for (row=0; row<NRows(); ++row) {
  811.                     range = info.block->GetRangeOfRow(row);
  812.                     info.block->SetRangeOfRow(row, range->from + actualShift, range->to);
  813.                     range = prevBlock->GetRangeOfRow(row);
  814.                     prevBlock->SetRangeOfRow(row, range->from, range->to + actualShift);
  815.                 }
  816.                 info.block->width -= actualShift;
  817.                 prevBlock->width += actualShift;
  818.                 if (prevBlock->width == 0) {
  819.                     RemoveBlock(prevBlock);
  820.                     TRACEMSG("removed empty block");
  821.                 }
  822.             }
  823.         }
  824.     }
  825.     if (actualShift != 0) {
  826.         UpdateBlockMapAndColors();
  827.         return true;
  828.     } else
  829.         return false;
  830. }
  831. bool BlockMultipleAlignment::ShiftRow(int row, int fromAlignmentIndex, int toAlignmentIndex,
  832.     eUnalignedJustification justification)
  833. {
  834.     if (fromAlignmentIndex == toAlignmentIndex) return false;
  835.     Block
  836.         *blockFrom = blockMap[fromAlignmentIndex].block,
  837.         *blockTo = blockMap[toAlignmentIndex].block;
  838.     // at least one end of the drag must be in an aligned block
  839.     UngappedAlignedBlock *ABlock =
  840.         dynamic_cast<UngappedAlignedBlock*>(blockFrom);
  841.     if (ABlock) {
  842.         if (blockTo != blockFrom && blockTo->IsAligned()) return false;
  843.     } else {
  844.         ABlock = dynamic_cast<UngappedAlignedBlock*>(blockTo);
  845.         if (!ABlock) return false;
  846.     }
  847.     // and the other must be in the same aligned block or an adjacent unaligned block
  848.     UnalignedBlock
  849.         *prevUABlock = dynamic_cast<UnalignedBlock*>(GetBlockBefore(ABlock)),
  850.         *nextUABlock = dynamic_cast<UnalignedBlock*>(GetBlockAfter(ABlock));
  851.     if (blockFrom != blockTo &&
  852.         ((ABlock == blockFrom && prevUABlock != blockTo && nextUABlock != blockTo) ||
  853.          (ABlock == blockTo && prevUABlock != blockFrom && nextUABlock != blockFrom)))
  854.         return false;
  855.     int requestedShift, actualShift = 0, width = 0;
  856.     // slightly different behaviour when dragging from unaligned to aligned...
  857.     if (!blockFrom->IsAligned()) {
  858.         int fromSeqIndex, toSeqIndex;
  859.         GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, NULL, &fromSeqIndex, NULL);
  860.         GetSequenceAndIndexAt(toAlignmentIndex, row, justification, NULL, &toSeqIndex, NULL);
  861.         if (fromSeqIndex < 0 || toSeqIndex < 0) return false;
  862.         requestedShift = toSeqIndex - fromSeqIndex;
  863.     }
  864.     // vs. dragging from aligned
  865.     else {
  866.         requestedShift = toAlignmentIndex - fromAlignmentIndex;
  867.     }
  868.     const Block::Range *prevRange = NULL, *nextRange = NULL,
  869.         *range = ABlock->GetRangeOfRow(row);
  870.     if (prevUABlock) prevRange = prevUABlock->GetRangeOfRow(row);
  871.     if (nextUABlock) nextRange = nextUABlock->GetRangeOfRow(row);
  872.     if (requestedShift > 0) {
  873.         if (prevUABlock) width = prevRange->to - prevRange->from + 1;
  874.         actualShift = (width > requestedShift) ? requestedShift : width;
  875.     } else {
  876.         if (nextUABlock) width = nextRange->to - nextRange->from + 1;
  877.         actualShift = (width > -requestedShift) ? requestedShift : -width;
  878.     }
  879.     if (actualShift == 0) return false;
  880.     TRACEMSG("shifting row " << row << " by " << actualShift);
  881.     ABlock->SetRangeOfRow(row, range->from - actualShift, range->to - actualShift);
  882.     if (prevUABlock) {
  883.         prevUABlock->SetRangeOfRow(row, prevRange->from, prevRange->to - actualShift);
  884.         prevUABlock->width = 0;
  885.         for (int i=0; i<NRows(); ++i) {
  886.             prevRange = prevUABlock->GetRangeOfRow(i);
  887.             width = prevRange->to - prevRange->from + 1;
  888.             if (width < 0) ERRORMSG("BlockMultipleAlignment::ShiftRow() - negative width on left");
  889.             if (width > prevUABlock->width) prevUABlock->width = width;
  890.         }
  891.         if (prevUABlock->width == 0) {
  892.             TRACEMSG("removing zero-width unaligned block on left");
  893.             RemoveBlock(prevUABlock);
  894.         }
  895.     } else {
  896.         TRACEMSG("creating unaligned block on left");
  897.         prevUABlock = CreateNewUnalignedBlockBetween(GetBlockBefore(ABlock), ABlock);
  898.         InsertBlockBefore(prevUABlock, ABlock);
  899.     }
  900.     if (nextUABlock) {
  901.         nextUABlock->SetRangeOfRow(row, nextRange->from - actualShift, nextRange->to);
  902.         nextUABlock->width = 0;
  903.         for (int i=0; i<NRows(); ++i) {
  904.             nextRange = nextUABlock->GetRangeOfRow(i);
  905.             width = nextRange->to - nextRange->from + 1;
  906.             if (width < 0) ERRORMSG("BlockMultipleAlignment::ShiftRow() - negative width on right");
  907.             if (width > nextUABlock->width) nextUABlock->width = width;
  908.         }
  909.         if (nextUABlock->width == 0) {
  910.             TRACEMSG("removing zero-width unaligned block on right");
  911.             RemoveBlock(nextUABlock);
  912.         }
  913.     } else {
  914.         TRACEMSG("creating unaligned block on right");
  915.         nextUABlock = CreateNewUnalignedBlockBetween(ABlock, GetBlockAfter(ABlock));
  916.         InsertBlockAfter(ABlock, nextUABlock);
  917.     }
  918.     if (!CheckAlignedBlock(ABlock))
  919.         ERRORMSG("BlockMultipleAlignment::ShiftRow() - shift failed to create valid aligned block");
  920.     UpdateBlockMapAndColors();
  921.     return true;
  922. }
  923. bool BlockMultipleAlignment::SplitBlock(int alignmentIndex)
  924. {
  925.     const BlockInfo& info = blockMap[alignmentIndex];
  926.     if (!info.block->IsAligned() || info.block->width < 2 || info.blockColumn == 0)
  927.         return false;
  928.     TRACEMSG("splitting block");
  929.     UngappedAlignedBlock *newAlignedBlock = new UngappedAlignedBlock(this);
  930.     newAlignedBlock->width = info.block->width - info.blockColumn;
  931.     info.block->width = info.blockColumn;
  932.     const Block::Range *range;
  933.     int oldTo;
  934.     for (int row=0; row<NRows(); ++row) {
  935.         range = info.block->GetRangeOfRow(row);
  936.         oldTo = range->to;
  937.         info.block->SetRangeOfRow(row, range->from, range->from + info.block->width - 1);
  938.         newAlignedBlock->SetRangeOfRow(row, oldTo - newAlignedBlock->width + 1, oldTo);
  939.     }
  940.     InsertBlockAfter(info.block, newAlignedBlock);
  941.     if (!CheckAlignedBlock(info.block) || !CheckAlignedBlock(newAlignedBlock))
  942.         ERRORMSG("BlockMultipleAlignment::SplitBlock() - split failed to create valid blocks");
  943.     UpdateBlockMapAndColors();
  944.     return true;
  945. }
  946. bool BlockMultipleAlignment::MergeBlocks(int fromAlignmentIndex, int toAlignmentIndex)
  947. {
  948.     Block
  949.         *expandedBlock = blockMap[fromAlignmentIndex].block,
  950.         *lastBlock = blockMap[toAlignmentIndex].block;
  951.     if (expandedBlock == lastBlock) return false;
  952.     int i;
  953.     for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i)
  954.         if (!blockMap[i].block->IsAligned()) return false;
  955.     TRACEMSG("merging block(s)");
  956.     for (i=0; i<NRows(); ++i)
  957.         expandedBlock->SetRangeOfRow(i,
  958.             expandedBlock->GetRangeOfRow(i)->from, lastBlock->GetRangeOfRow(i)->to);
  959.     expandedBlock->width = lastBlock->GetRangeOfRow(0)->to - expandedBlock->GetRangeOfRow(0)->from + 1;
  960.     Block *deletedBlock = NULL, *blockToDelete;
  961.     for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i) {
  962.         blockToDelete = blockMap[i].block;
  963.         if (blockToDelete == expandedBlock) continue;
  964.         if (blockToDelete != deletedBlock) {
  965.             deletedBlock = blockToDelete;
  966.             RemoveBlock(blockToDelete);
  967.         }
  968.     }
  969.     if (!CheckAlignedBlock(expandedBlock))
  970.         ERRORMSG("BlockMultipleAlignment::MergeBlocks() - merge failed to create valid block");
  971.     UpdateBlockMapAndColors();
  972.     return true;
  973. }
  974. bool BlockMultipleAlignment::CreateBlock(int fromAlignmentIndex, int toAlignmentIndex,
  975.     eUnalignedJustification justification)
  976. {
  977.     const BlockInfo& info = blockMap[fromAlignmentIndex];
  978.     UnalignedBlock *prevUABlock = dynamic_cast<UnalignedBlock*>(info.block);
  979.     if (!prevUABlock || info.block != blockMap[toAlignmentIndex].block) return false;
  980.     int row, seqIndexFrom, seqIndexTo,
  981.         newBlockWidth = toAlignmentIndex - fromAlignmentIndex + 1,
  982.         origWidth = prevUABlock->width;
  983.     vector < int > seqIndexesFrom(NRows()), seqIndexesTo(NRows());
  984.     const Sequence *seq;
  985. bool ignored;
  986.     for (row=0; row<NRows(); ++row) {
  987.         if (!GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, &seq, &seqIndexFrom, &ignored) ||
  988.             !GetSequenceAndIndexAt(toAlignmentIndex, row, justification, &seq, &seqIndexTo, &ignored) ||
  989.             seqIndexFrom < 0 || seqIndexTo < 0 ||
  990.             seqIndexTo - seqIndexFrom + 1 != newBlockWidth) return false;
  991.         seqIndexesFrom[row] = seqIndexFrom;
  992.         seqIndexesTo[row] = seqIndexTo;
  993.     }
  994.     TRACEMSG("creating new aligned and unaligned blocks");
  995.     UnalignedBlock *nextUABlock = new UnalignedBlock(this);
  996.     UngappedAlignedBlock *ABlock = new UngappedAlignedBlock(this);
  997.     prevUABlock->width = nextUABlock->width = 0;
  998.     bool deletePrevUABlock = true, deleteNextUABlock = true;
  999.     const Block::Range *prevRange;
  1000.     int rangeWidth;
  1001.     for (row=0; row<NRows(); ++row) {
  1002.         prevRange = prevUABlock->GetRangeOfRow(row);
  1003.         nextUABlock->SetRangeOfRow(row, seqIndexesTo[row] + 1, prevRange->to);
  1004.         rangeWidth = prevRange->to - seqIndexesTo[row];
  1005.         if (rangeWidth < 0)
  1006.             ERRORMSG("BlockMultipleAlignment::CreateBlock() - negative nextRange width");
  1007.         else if (rangeWidth > 0) {
  1008.             if (rangeWidth > nextUABlock->width) nextUABlock->width = rangeWidth;
  1009.             deleteNextUABlock = false;
  1010.         }
  1011.         prevUABlock->SetRangeOfRow(row, prevRange->from, seqIndexesFrom[row] - 1);
  1012.         rangeWidth = seqIndexesFrom[row] - prevRange->from;
  1013.         if (rangeWidth < 0)
  1014.             ERRORMSG("BlockMultipleAlignment::CreateBlock() - negative prevRange width");
  1015.         else if (rangeWidth > 0) {
  1016.             if (rangeWidth > prevUABlock->width) prevUABlock->width = rangeWidth;
  1017.             deletePrevUABlock = false;
  1018.         }
  1019.         ABlock->SetRangeOfRow(row, seqIndexesFrom[row], seqIndexesTo[row]);
  1020.     }
  1021.     ABlock->width = newBlockWidth;
  1022.     if (prevUABlock->width + ABlock->width + nextUABlock->width != origWidth)
  1023.         ERRORMSG("BlockMultipleAlignment::CreateBlock() - bad block widths sum");
  1024.     InsertBlockAfter(prevUABlock, ABlock);
  1025.     InsertBlockAfter(ABlock, nextUABlock);
  1026.     if (deletePrevUABlock) {
  1027.         TRACEMSG("deleting zero-width unaligned block on left");
  1028.         RemoveBlock(prevUABlock);
  1029.     }
  1030.     if (deleteNextUABlock) {
  1031.         TRACEMSG("deleting zero-width unaligned block on right");
  1032.         RemoveBlock(nextUABlock);
  1033.     }
  1034.     if (!CheckAlignedBlock(ABlock))
  1035.         ERRORMSG("BlockMultipleAlignment::CreateBlock() - failed to create valid block");
  1036.     UpdateBlockMapAndColors();
  1037.     return true;
  1038. }
  1039. bool BlockMultipleAlignment::DeleteBlock(int alignmentIndex)
  1040. {
  1041.     Block *block = blockMap[alignmentIndex].block;
  1042.     if (!block || !block->IsAligned()) return false;
  1043.     TRACEMSG("deleting block");
  1044.     Block
  1045.         *prevBlock = GetBlockBefore(block),
  1046.         *nextBlock = GetBlockAfter(block);
  1047.     // unaligned blocks on both sides - note that total alignment width can change!
  1048.     if (prevBlock && !prevBlock->IsAligned() && nextBlock && !nextBlock->IsAligned()) {
  1049.         const Block::Range *prevRange, *nextRange;
  1050.         int maxWidth = 0, width;
  1051.         for (int row=0; row<NRows(); ++row) {
  1052.             prevRange = prevBlock->GetRangeOfRow(row);
  1053.             nextRange = nextBlock->GetRangeOfRow(row);
  1054.             width = nextRange->to - prevRange->from + 1;
  1055.             prevBlock->SetRangeOfRow(row, prevRange->from, nextRange->to);
  1056.             if (width > maxWidth) maxWidth = width;
  1057.         }
  1058.         prevBlock->width = maxWidth;
  1059.         TRACEMSG("removing extra unaligned block");
  1060.         RemoveBlock(nextBlock);
  1061.     }
  1062.     // unaligned block on left only
  1063.     else if (prevBlock && !prevBlock->IsAligned()) {
  1064.         const Block::Range *prevRange, *range;
  1065.         for (int row=0; row<NRows(); ++row) {
  1066.             prevRange = prevBlock->GetRangeOfRow(row);
  1067.             range = block->GetRangeOfRow(row);
  1068.             prevBlock->SetRangeOfRow(row, prevRange->from, range->to);
  1069.         }
  1070.         prevBlock->width += block->width;
  1071.     }
  1072.     // unaligned block on right only
  1073.     else if (nextBlock && !nextBlock->IsAligned()) {
  1074.         const Block::Range *range, *nextRange;
  1075.         for (int row=0; row<NRows(); ++row) {
  1076.             range = block->GetRangeOfRow(row);
  1077.             nextRange = nextBlock->GetRangeOfRow(row);
  1078.             nextBlock->SetRangeOfRow(row, range->from, nextRange->to);
  1079.         }
  1080.         nextBlock->width += block->width;
  1081.     }
  1082.     // no adjacent unaligned blocks
  1083.     else {
  1084.         TRACEMSG("creating new unaligned block");
  1085.         Block *newBlock = CreateNewUnalignedBlockBetween(prevBlock, nextBlock);
  1086.         if (prevBlock)
  1087.             InsertBlockAfter(prevBlock, newBlock);
  1088.         else if (nextBlock)
  1089.             InsertBlockBefore(newBlock, nextBlock);
  1090.         else
  1091.             blocks.push_back(newBlock);
  1092.     }
  1093.     RemoveBlock(block);
  1094.     UpdateBlockMapAndColors();
  1095.     return true;
  1096. }
  1097. bool BlockMultipleAlignment::DeleteAllBlocks(void)
  1098. {
  1099.     if (blocks.size() == 0) return false;
  1100.     DELETE_ALL_AND_CLEAR(blocks, BlockList);
  1101.     InitCache();
  1102.     AddUnalignedBlocks();   // one single unaligned block for whole alignment
  1103.     UpdateBlockMapAndColors();
  1104.     return true;
  1105. }
  1106. bool BlockMultipleAlignment::DeleteRow(int row)
  1107. {
  1108.     if (row < 0 || row >= NRows()) {
  1109.         ERRORMSG("BlockMultipleAlignment::DeleteRow() - row out of range");
  1110.         return false;
  1111.     }
  1112.     // remove sequence from list
  1113.     SequenceList::iterator s = (const_cast<SequenceList*>(sequences))->begin();
  1114.     for (int i=0; i<row; ++i) ++s;
  1115.     (const_cast<SequenceList*>(sequences))->erase(s);
  1116.     // delete row from all blocks, removing any zero-width blocks
  1117.     BlockList::iterator b = blocks.begin(), br, be = blocks.end();
  1118.     while (b != be) {
  1119.         (*b)->DeleteRow(row);
  1120.         if ((*b)->width == 0) {
  1121.             br = b;
  1122.             ++b;
  1123.             TRACEMSG("deleting block resized to zero width");
  1124.             RemoveBlock(*br);
  1125.         } else
  1126.             ++b;
  1127.     }
  1128.     // update total alignment width
  1129.     UpdateBlockMapAndColors();
  1130.     InitCache();
  1131.     return true;
  1132. }
  1133. void BlockMultipleAlignment::GetUngappedAlignedBlocks(UngappedAlignedBlockList *uabs) const
  1134. {
  1135.     uabs->clear();
  1136.     uabs->reserve(blocks.size());
  1137.     BlockList::const_iterator b, be = blocks.end();
  1138.     for (b=blocks.begin(); b!=be; ++b) {
  1139.         UngappedAlignedBlock *uab = dynamic_cast<UngappedAlignedBlock*>(*b);
  1140.         if (uab) uabs->push_back(uab);
  1141.     }
  1142.     uabs->resize(uabs->size());
  1143. }
  1144. bool BlockMultipleAlignment::ExtractRows(
  1145.     const vector < int >& slavesToRemove, AlignmentList *pairwiseAlignments)
  1146. {
  1147.     if (slavesToRemove.size() == 0) return false;
  1148.     // make a bool list of rows to remove, also checking to make sure slave list items are in range
  1149.     int i;
  1150.     vector < bool > removeRows(NRows(), false);
  1151.     for (i=0; i<slavesToRemove.size(); ++i) {
  1152.         if (slavesToRemove[i] > 0 && slavesToRemove[i] < NRows()) {
  1153.             removeRows[slavesToRemove[i]] = true;
  1154.         } else {
  1155.             ERRORMSG("BlockMultipleAlignment::ExtractRows() - can't extract row "
  1156.                 << slavesToRemove[i]);
  1157.             return false;
  1158.         }
  1159.     }
  1160.     if (pairwiseAlignments) {
  1161.         TRACEMSG("creating new pairwise alignments");
  1162.         SetDiagPostLevel(eDiag_Warning);    // otherwise, info messages take a long time if lots of rows
  1163.         UngappedAlignedBlockList uaBlocks;
  1164.         GetUngappedAlignedBlocks(&uaBlocks);
  1165.         UngappedAlignedBlockList::const_iterator u, ue = uaBlocks.end();
  1166.         for (i=0; i<slavesToRemove.size(); ++i) {
  1167.             // redraw molecule associated with removed row
  1168.             const Molecule *molecule = GetSequenceOfRow(slavesToRemove[i])->molecule;
  1169.             if (molecule) GlobalMessenger()->PostRedrawMolecule(molecule);
  1170.             // create new pairwise alignment from each removed row
  1171.             SequenceList *newSeqs = new SequenceList(2);
  1172.             (*newSeqs)[0] = (*sequences)[0];
  1173.             (*newSeqs)[1] = (*sequences)[slavesToRemove[i]];
  1174.             BlockMultipleAlignment *newAlignment = new BlockMultipleAlignment(newSeqs, alignmentManager);
  1175.             for (u=uaBlocks.begin(); u!=ue; ++u) {
  1176.                 // only copy blocks that aren't flagged to be realigned
  1177.                 if (markBlocks.find(*u) == markBlocks.end()) {
  1178.                     UngappedAlignedBlock *newABlock = new UngappedAlignedBlock(newAlignment);
  1179.                     const Block::Range *range = (*u)->GetRangeOfRow(0);
  1180.                     newABlock->SetRangeOfRow(0, range->from, range->to);
  1181.                     range = (*u)->GetRangeOfRow(slavesToRemove[i]);
  1182.                     newABlock->SetRangeOfRow(1, range->from, range->to);
  1183.                     newABlock->width = range->to - range->from + 1;
  1184.                     newAlignment->AddAlignedBlockAtEnd(newABlock);
  1185.                 }
  1186.             }
  1187.             if (!newAlignment->AddUnalignedBlocks() ||
  1188.                 !newAlignment->UpdateBlockMapAndColors()) {
  1189.                 ERRORMSG("BlockMultipleAlignment::ExtractRows() - error creating new alignment");
  1190.                 return false;
  1191.             }
  1192.             // add aligned region info (for threader to use later on)
  1193.             if (uaBlocks.size() > 0) {
  1194.                 int excess = 0;
  1195.                 if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
  1196.                     WARNINGMSG("Can't get footprint excess residues from registry");
  1197.                 newAlignment->alignSlaveFrom =
  1198.                     uaBlocks.front()->GetRangeOfRow(slavesToRemove[i])->from - excess;
  1199.                 if (newAlignment->alignSlaveFrom < 0)
  1200.                     newAlignment->alignSlaveFrom = 0;
  1201.                 newAlignment->alignSlaveTo =
  1202.                     uaBlocks.back()->GetRangeOfRow(slavesToRemove[i])->to + excess;
  1203.                 if (newAlignment->alignSlaveTo >= (*newSeqs)[1]->Length())
  1204.                     newAlignment->alignSlaveTo = (*newSeqs)[1]->Length() - 1;
  1205.                 TRACEMSG((*newSeqs)[1]->identifier->ToString() << " aligned from "
  1206.                     << newAlignment->alignSlaveFrom << " to " << newAlignment->alignSlaveTo);
  1207.             }
  1208.             pairwiseAlignments->push_back(newAlignment);
  1209.         }
  1210.         SetDiagPostLevel(eDiag_Info);
  1211.     }
  1212.     // remove sequences
  1213.     TRACEMSG("deleting sequences");
  1214.     VectorRemoveElements(*(const_cast<SequenceList*>(sequences)), removeRows, slavesToRemove.size());
  1215.     VectorRemoveElements(rowDoubles, removeRows, slavesToRemove.size());
  1216.     VectorRemoveElements(rowStrings, removeRows, slavesToRemove.size());
  1217.     // delete row from all blocks, removing any zero-width blocks
  1218.     TRACEMSG("deleting alignment rows from blocks");
  1219.     BlockList::const_iterator b = blocks.begin(), br, be = blocks.end();
  1220.     while (b != be) {
  1221.         (*b)->DeleteRows(removeRows, slavesToRemove.size());
  1222.         if ((*b)->width == 0) {
  1223.             br = b;
  1224.             ++b;
  1225.             TRACEMSG("deleting block resized to zero width");
  1226.             RemoveBlock(*br);
  1227.         } else
  1228.             ++b;
  1229.     }
  1230.     // update total alignment width
  1231.     UpdateBlockMapAndColors();
  1232.     InitCache();
  1233.     return true;
  1234. }
  1235. bool BlockMultipleAlignment::MergeAlignment(const BlockMultipleAlignment *newAlignment)
  1236. {
  1237.     // check to see if the alignment is compatible - must have same master sequence
  1238.     // and blocks of new alignment must contain blocks of this alignment; at same time,
  1239.     // build up map of aligned blocks in new alignment that correspond to aligned blocks
  1240.     // of this object, for convenient lookup later
  1241.     if (newAlignment->GetMaster() != GetMaster()) return false;
  1242.     const Block::Range *newRange, *thisRange;
  1243.     BlockList::const_iterator nb, nbe = newAlignment->blocks.end();
  1244.     BlockList::iterator b, be = blocks.end();
  1245.     typedef map < UngappedAlignedBlock *, const UngappedAlignedBlock * > AlignedBlockMap;
  1246.     AlignedBlockMap correspondingNewBlocks;
  1247.     for (b=blocks.begin(); b!=be; ++b) {
  1248.         if (!(*b)->IsAligned()) continue;
  1249.         for (nb=newAlignment->blocks.begin(); nb!=nbe; ++nb) {
  1250.             if (!(*nb)->IsAligned()) continue;
  1251.             newRange = (*nb)->GetRangeOfRow(0);
  1252.             thisRange = (*b)->GetRangeOfRow(0);
  1253.             if (newRange->from <= thisRange->from && newRange->to >= thisRange->to) {
  1254.                 correspondingNewBlocks[dynamic_cast<UngappedAlignedBlock*>(*b)] =
  1255.                     dynamic_cast<const UngappedAlignedBlock*>(*nb);
  1256.                 break;
  1257.             }
  1258.         }
  1259.         if (nb == nbe) return false;    // no corresponding block found
  1260.     }
  1261.     // add slave sequences from new alignment; also copy scores/status
  1262.     SequenceList *modSequences = const_cast<SequenceList*>(sequences);
  1263.     int i, nNewRows = newAlignment->sequences->size() - 1;
  1264.     modSequences->resize(modSequences->size() + nNewRows);
  1265.     rowDoubles.resize(rowDoubles.size() + nNewRows);
  1266.     rowStrings.resize(rowStrings.size() + nNewRows);
  1267.     for (i=0; i<nNewRows; ++i) {
  1268.         (*modSequences)[modSequences->size() + i - nNewRows] = (*(newAlignment->sequences))[i + 1];
  1269.         SetRowDouble(NRows() + i - nNewRows, newAlignment->GetRowDouble(i + 1));
  1270.         SetRowStatusLine(NRows() + i - nNewRows, newAlignment->GetRowStatusLine(i + 1));
  1271.     }
  1272.     // now that we know blocks are compatible, add new rows at end of this alignment, containing
  1273.     // all rows (except master) from new alignment; only that part of the aligned blocks from
  1274.     // the new alignment that intersect with the aligned blocks from this object are added, so
  1275.     // that this object's block structure is unchanged
  1276.     // resize all aligned blocks
  1277.     for (b=blocks.begin(); b!=be; ++b) (*b)->AddRows(nNewRows);
  1278.     // set ranges of aligned blocks, and add them to the list
  1279.     AlignedBlockMap::const_iterator ab, abe = correspondingNewBlocks.end();
  1280.     const Block::Range *thisMaster, *newMaster;
  1281.     for (ab=correspondingNewBlocks.begin(); ab!=abe; ++ab) {
  1282.         thisMaster = ab->first->GetRangeOfRow(0);
  1283.         newMaster = ab->second->GetRangeOfRow(0);
  1284.         for (i=0; i<nNewRows; ++i) {
  1285.             newRange = ab->second->GetRangeOfRow(i + 1);
  1286.             ab->first->SetRangeOfRow(NRows() + i - nNewRows,
  1287.                 newRange->from + thisMaster->from - newMaster->from,
  1288.                 newRange->to + thisMaster->to - newMaster->to);
  1289.         }
  1290.     }
  1291.     // delete then recreate the unaligned blocks, in case a merge requires the
  1292.     // creation of a new unaligned block
  1293.     for (b=blocks.begin(); b!=be; ) {
  1294.         if (!(*b)->IsAligned()) {
  1295.             BlockList::iterator bb(b);
  1296.             ++bb;
  1297.             delete *b;
  1298.             blocks.erase(b);
  1299.             b = bb;
  1300.         } else
  1301.             ++b;
  1302.     }
  1303.     InitCache();
  1304.     // update this alignment, but leave row scores/status alone
  1305.     if (!AddUnalignedBlocks() || !UpdateBlockMapAndColors(false)) {
  1306.         ERRORMSG("BlockMultipleAlignment::MergeAlignment() - internal update after merge failed");
  1307.         return false;
  1308.     }
  1309.     return true;
  1310. }
  1311. int BlockMultipleAlignment::ShowGeometryViolations(bool showGV)
  1312. {
  1313.     geometryViolations.clear();
  1314.     if (!showGV || !GetMaster()->molecule || GetMaster()->molecule->parentSet->isAlphaOnly) {
  1315.         showGeometryViolations = false;
  1316.         return 0;
  1317.     }
  1318.     Threader::GeometryViolationsForRow violations;
  1319.     int nViolations = alignmentManager->threader->GetGeometryViolations(this, &violations);
  1320.     if (violations.size() != NRows()) {
  1321.         ERRORMSG("BlockMultipleAlignment::ShowGeometryViolations() - wrong size list");
  1322.         showGeometryViolations = false;
  1323.         return 0;
  1324.     }
  1325.     geometryViolations.resize(NRows());
  1326.     for (int row=0; row<NRows(); ++row) {
  1327.         geometryViolations[row].resize(GetSequenceOfRow(row)->Length(), false);
  1328.         Threader::IntervalList::const_iterator i, ie = violations[row].end();
  1329.         for (i=violations[row].begin(); i!=ie; ++i)
  1330.             for (int l=i->first; l<=i->second; ++l)
  1331.                 geometryViolations[row][l] = true;
  1332.     }
  1333.     showGeometryViolations = true;
  1334.     return nViolations;
  1335. }
  1336. CSeq_align * CreatePairwiseSeqAlignFromMultipleRow(const BlockMultipleAlignment *multiple,
  1337.     const BlockMultipleAlignment::UngappedAlignedBlockList& blocks, int slaveRow)
  1338. {
  1339.     if (!multiple || slaveRow < 0 || slaveRow >= multiple->NRows()) {
  1340.         ERRORMSG("CreatePairwiseSeqAlignFromMultipleRow() - bad parameters");
  1341.         return NULL;
  1342.     }
  1343.     CSeq_align *seqAlign = new CSeq_align();
  1344.     seqAlign->SetType(CSeq_align::eType_partial);
  1345.     seqAlign->SetDim(2);
  1346.     CSeq_align::C_Segs::TDendiag& denDiags = seqAlign->SetSegs().SetDendiag();
  1347.     denDiags.resize((blocks.size() > 0) ? blocks.size() : 1);
  1348.     CSeq_align::C_Segs::TDendiag::iterator d, de = denDiags.end();
  1349.     BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b = blocks.begin();
  1350.     const Block::Range *range;
  1351.     for (d=denDiags.begin(); d!=de; ++d, ++b) {
  1352.         CDense_diag *denDiag = new CDense_diag();
  1353.         d->Reset(denDiag);
  1354.         denDiag->SetDim(2);
  1355.         denDiag->SetIds().resize(2);
  1356.         // master row
  1357.         denDiag->SetIds().front().Reset(multiple->GetSequenceOfRow(0)->CreateSeqId());
  1358.         if (blocks.size() > 0) {
  1359.             range = (*b)->GetRangeOfRow(0);
  1360.             denDiag->SetStarts().push_back(range->from);
  1361.         } else
  1362.             denDiag->SetStarts().push_back(0);
  1363.         // slave row
  1364.         denDiag->SetIds().back().Reset(multiple->GetSequenceOfRow(slaveRow)->CreateSeqId());
  1365.         if (blocks.size() > 0) {
  1366.             range = (*b)->GetRangeOfRow(slaveRow);
  1367.             denDiag->SetStarts().push_back(range->from);
  1368.         } else
  1369.             denDiag->SetStarts().push_back(0);
  1370.         // block width
  1371.         denDiag->SetLen((blocks.size() > 0) ? (*b)->width : 0);
  1372.     }
  1373.     return seqAlign;
  1374. }
  1375. bool BlockMultipleAlignment::MarkBlock(int column)
  1376. {
  1377.     if (column >= 0 && column < blockMap.size() && blockMap[column].block->IsAligned()) {
  1378.         TRACEMSG("marked block for realignment");
  1379.         markBlocks[blockMap[column].block] = true;
  1380.         return true;
  1381.     }
  1382.     return false;
  1383. }
  1384. bool BlockMultipleAlignment::ClearMarks(void)
  1385. {
  1386.     if (markBlocks.size() == 0) return false;
  1387.     markBlocks.clear();
  1388.     return true;
  1389. }
  1390. bool BlockMultipleAlignment::HighlightAlignedColumnsOfMasterRange(int from, int to) const
  1391. {
  1392.     const Sequence *master = GetMaster();
  1393.     // must do one column at a time, rather than range, in case there are inserts wrt the master
  1394.     bool anyError = false;
  1395.     for (int i=from; i<=to; ++i) {
  1396.         // sanity check
  1397.         if (i < 0 || i >= master->Length() || !IsAligned(0, i)) {
  1398.             WARNINGMSG("Can't highlight alignment at master residue " << (i+1));
  1399.             anyError = true;
  1400.             // highlight unaligned residues, but master only
  1401.             if (i >= 0 && i < master->Length())
  1402.                 GlobalMessenger()->AddHighlights(GetSequenceOfRow(0), i, i);
  1403.             continue;
  1404.         }
  1405.         // get block and offset
  1406.         const Block *block = GetBlock(0, i);
  1407.         int blockOffset = i - block->GetRangeOfRow(0)->from;
  1408.         // highlight aligned residue in each row
  1409.         for (int row=0; row<NRows(); ++row) {
  1410.             int slaveIndex = block->GetRangeOfRow(row)->from + blockOffset;
  1411.             GlobalMessenger()->AddHighlights(GetSequenceOfRow(row), slaveIndex, slaveIndex);
  1412.         }
  1413.     }
  1414.     return !anyError;
  1415. }
  1416. int BlockMultipleAlignment::NAlignedBlocks(void) const
  1417. {
  1418.     int n = 0;
  1419.     BlockList::const_iterator b, be = blocks.end();
  1420.     for (b=blocks.begin(); b!=be; ++b) if ((*b)->IsAligned()) ++n;
  1421.     return n;
  1422. }
  1423. int BlockMultipleAlignment::GetAlignmentIndex(int row, int seqIndex, eUnalignedJustification justification)
  1424. {
  1425.     if (row < 0 || row >= NRows() || seqIndex < 0 || seqIndex >= GetSequenceOfRow(row)->Length()) {
  1426.         ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - coordinate out of range");
  1427.         return -1;
  1428.     }
  1429.     int alignmentIndex, blockColumn;
  1430.     const Block *block = NULL;
  1431.     const Block::Range *range;
  1432.     for (alignmentIndex=0; alignmentIndex<totalWidth; ++alignmentIndex) {
  1433.         // check each block to see if index is in range
  1434.         if (block != blockMap[alignmentIndex].block) {
  1435.             block = blockMap[alignmentIndex].block;
  1436.             range = block->GetRangeOfRow(row);
  1437.             if (seqIndex >= range->from && seqIndex <= range->to) {
  1438.                 // override requested justification for end blocks
  1439.                 if (block == blocks.back()) // also true if there's a single aligned block
  1440.                     justification = eLeft;
  1441.                 else if (block == blocks.front())
  1442.                     justification = eRight;
  1443.                 // search block columns to find index (inefficient, but avoids having to write
  1444.                 // inverse functions of Block::GetIndexAt()
  1445.                 for (blockColumn=0; blockColumn<block->width; ++blockColumn) {
  1446.                     if (seqIndex == block->GetIndexAt(blockColumn, row, justification))
  1447.                         return alignmentIndex + blockColumn;
  1448.                 }
  1449.                 ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - can't find index in block");
  1450.                 return -1;
  1451.             }
  1452.         }
  1453.     }
  1454.     // should never get here...
  1455.     ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - confused");
  1456.     return -1;
  1457. }
  1458. // creates a SeqAlign from a BlockMultipleAlignment
  1459. SeqAlignPtr BlockMultipleAlignment::CreateCSeqAlign(void) const
  1460. {
  1461.     // one SeqAlign (chained into a linked list) for each slave row
  1462.     SeqAlignPtr prevSap = NULL, firstSap = NULL;
  1463.     for (int row=1; row<NRows(); ++row) {
  1464.         SeqAlignPtr sap = SeqAlignNew();
  1465.         if (prevSap) prevSap->next = sap;
  1466.         prevSap = sap;
  1467. if (!firstSap) firstSap = sap;
  1468.         sap->type = SAT_PARTIAL;
  1469.         sap->dim = 2;
  1470.         sap->segtype = SAS_DENDIAG;
  1471.         DenseDiagPtr prevDd = NULL;
  1472.         UngappedAlignedBlockList blocks;
  1473.         GetUngappedAlignedBlocks(&blocks);
  1474.         UngappedAlignedBlockList::const_iterator b, be = blocks.end();
  1475.         for (b=blocks.begin(); b!=be; ++b) {
  1476.             DenseDiagPtr dd = DenseDiagNew();
  1477.             if (prevDd) prevDd->next = dd;
  1478.             prevDd = dd;
  1479.             if (b == blocks.begin()) sap->segs = dd;
  1480.             dd->dim = 2;
  1481.             GetSequenceOfRow(0)->AddCSeqId(&(dd->id), false);      // master
  1482.             GetSequenceOfRow(row)->AddCSeqId(&(dd->id), false);    // slave
  1483.             dd->len = (*b)->width;
  1484.             dd->starts = (Int4Ptr) MemNew(2 * sizeof(Int4));
  1485.             const Block::Range *range = (*b)->GetRangeOfRow(0);
  1486.             dd->starts[0] = range->from;
  1487.             range = (*b)->GetRangeOfRow(row);
  1488.             dd->starts[1] = range->from;
  1489.         }
  1490.     }
  1491. return firstSap;
  1492. }
  1493. ///// UngappedAlignedBlock methods /////
  1494. char UngappedAlignedBlock::GetCharacterAt(int blockColumn, int row) const
  1495. {
  1496.     return (*(parentAlignment->sequences))[row]->sequenceString[GetIndexAt(blockColumn, row)];
  1497. }
  1498. Block * UngappedAlignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
  1499. {
  1500.     UngappedAlignedBlock *copy = new UngappedAlignedBlock(newMultiple);
  1501.     const Block::Range *range;
  1502.     for (int row=0; row<NSequences(); ++row) {
  1503.         range = GetRangeOfRow(row);
  1504.         copy->SetRangeOfRow(row, range->from, range->to);
  1505.     }
  1506.     copy->width = width;
  1507.     return copy;
  1508. }
  1509. void UngappedAlignedBlock::DeleteRow(int row)
  1510. {
  1511.     RangeList::iterator r = ranges.begin();
  1512.     for (int i=0; i<row; ++i) ++r;
  1513.     ranges.erase(r);
  1514. }
  1515. void UngappedAlignedBlock::DeleteRows(vector < bool >& removeRows, int nToRemove)
  1516. {
  1517.     VectorRemoveElements(ranges, removeRows, nToRemove);
  1518. }
  1519. ///// UnalignedBlock methods /////
  1520. int UnalignedBlock::GetIndexAt(int blockColumn, int row,
  1521.         BlockMultipleAlignment::eUnalignedJustification justification) const
  1522. {
  1523.     const Block::Range *range = GetRangeOfRow(row);
  1524.     int seqIndex, rangeWidth, rangeMiddle, extraSpace;
  1525.     switch (justification) {
  1526.         case BlockMultipleAlignment::eLeft:
  1527.             seqIndex = range->from + blockColumn;
  1528.             break;
  1529.         case BlockMultipleAlignment::eRight:
  1530.             seqIndex = range->to - width + blockColumn + 1;
  1531.             break;
  1532.         case BlockMultipleAlignment::eCenter:
  1533.             rangeWidth = (range->to - range->from + 1);
  1534.             extraSpace = (width - rangeWidth) / 2;
  1535.             if (blockColumn < extraSpace || blockColumn >= extraSpace + rangeWidth)
  1536.                 seqIndex = -1;
  1537.             else
  1538.                 seqIndex = range->from + blockColumn - extraSpace;
  1539.             break;
  1540.         case BlockMultipleAlignment::eSplit:
  1541.             rangeWidth = (range->to - range->from + 1);
  1542.             rangeMiddle = (rangeWidth / 2) + (rangeWidth % 2);
  1543.             extraSpace = width - rangeWidth;
  1544.             if (blockColumn < rangeMiddle)
  1545.                 seqIndex = range->from + blockColumn;
  1546.             else if (blockColumn >= extraSpace + rangeMiddle)
  1547.                 seqIndex = range->to - width + blockColumn + 1;
  1548.             else
  1549.                 seqIndex = -1;
  1550.             break;
  1551.     }
  1552.     if (seqIndex < range->from || seqIndex > range->to) seqIndex = -1;
  1553.     return seqIndex;
  1554. }
  1555. void UnalignedBlock::Resize(void)
  1556. {
  1557.     width = 0;
  1558.     for (int i=0; i<NSequences(); ++i) {
  1559.         int blockWidth = ranges[i].to - ranges[i].from + 1;
  1560.         if (blockWidth > width) width = blockWidth;
  1561.     }
  1562. }
  1563. void UnalignedBlock::DeleteRow(int row)
  1564. {
  1565.     RangeList::iterator r = ranges.begin();
  1566.     for (int i=0; i<row; ++i) ++r;
  1567.     ranges.erase(r);
  1568.     Resize();
  1569. }
  1570. void UnalignedBlock::DeleteRows(vector < bool >& removeRows, int nToRemove)
  1571. {
  1572.     VectorRemoveElements(ranges, removeRows, nToRemove);
  1573.     Resize();
  1574. }
  1575. Block * UnalignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
  1576. {
  1577.     UnalignedBlock *copy = new UnalignedBlock(newMultiple);
  1578.     const Block::Range *range;
  1579.     for (int row=0; row<NSequences(); ++row) {
  1580.         range = GetRangeOfRow(row);
  1581.         copy->SetRangeOfRow(row, range->from, range->to);
  1582.     }
  1583.     copy->width = width;
  1584.     return copy;
  1585. }
  1586. END_SCOPE(Cn3D)
  1587. /*
  1588. * ---------------------------------------------------------------------------
  1589. * $Log: block_multiple_alignment.cpp,v $
  1590. * Revision 1000.3  2004/06/01 18:27:49  gouriano
  1591. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.57
  1592. *
  1593. * Revision 1.57  2004/05/26 14:18:29  thiessen
  1594. * remove status stuff in ExtractRows
  1595. *
  1596. * Revision 1.56  2004/05/21 21:41:38  gorelenk
  1597. * Added PCH ncbi_pch.hpp
  1598. *
  1599. * Revision 1.55  2004/03/15 17:27:02  thiessen
  1600. * prefer prefix vs. postfix ++/-- operators
  1601. *
  1602. * Revision 1.54  2004/02/19 17:04:43  thiessen
  1603. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  1604. *
  1605. * Revision 1.53  2003/11/06 18:52:31  thiessen
  1606. * make geometry violations shown on/off; allow multiple pmid entry in ref dialog
  1607. *
  1608. * Revision 1.52  2003/07/14 18:37:07  thiessen
  1609. * change GetUngappedAlignedBlocks() param types; other syntax changes
  1610. *
  1611. * Revision 1.51  2003/06/13 14:46:58  thiessen
  1612. * fix row cache bug
  1613. *
  1614. * Revision 1.50  2003/06/12 14:21:11  thiessen
  1615. * when saving single-row alignment, make master-master alignment in asn
  1616. *
  1617. * Revision 1.49  2003/02/06 16:39:52  thiessen
  1618. * add block row fit coloring
  1619. *
  1620. * Revision 1.48  2003/02/03 19:20:00  thiessen
  1621. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  1622. *
  1623. * Revision 1.47  2003/01/31 17:18:58  thiessen
  1624. * many small additions and changes...
  1625. *
  1626. * Revision 1.46  2003/01/30 14:00:23  thiessen
  1627. * add Block Z Fit coloring
  1628. *
  1629. * Revision 1.45  2003/01/28 21:07:56  thiessen
  1630. * add block fit coloring algorithm; tweak row dragging; fix style bug
  1631. *
  1632. * Revision 1.44  2003/01/23 20:03:05  thiessen
  1633. * add BLAST Neighbor algorithm
  1634. *
  1635. * Revision 1.43  2003/01/22 14:47:30  thiessen
  1636. * cache PSSM in BlockMultipleAlignment
  1637. *
  1638. * Revision 1.42  2002/12/13 15:33:37  thiessen
  1639. * tweak to HighlightAlignedColumnsOfMasterRange
  1640. *
  1641. * Revision 1.41  2002/12/12 23:07:18  thiessen
  1642. * improved handling of alignment annotation errors
  1643. *
  1644. * Revision 1.40  2002/11/18 20:49:11  thiessen
  1645. * move unaligned/no-coord colors into Colors class
  1646. *
  1647. * Revision 1.39  2002/10/23 01:29:25  thiessen
  1648. * fix block cache bug
  1649. *
  1650. * Revision 1.38  2002/08/28 20:30:33  thiessen
  1651. * fix proximity sort bug
  1652. *
  1653. * Revision 1.37  2002/07/26 15:28:44  thiessen
  1654. * add Alejandro's block alignment algorithm
  1655. *
  1656. * Revision 1.36  2002/03/28 14:06:02  thiessen
  1657. * preliminary BLAST/PSSM ; new CD startup style
  1658. *
  1659. * Revision 1.35  2002/02/21 22:01:49  thiessen
  1660. * remember alignment range on demotion
  1661. *
  1662. * Revision 1.34  2002/02/19 14:59:38  thiessen
  1663. * add CDD reject and purge sequence
  1664. *
  1665. * Revision 1.33  2002/02/05 18:53:24  thiessen
  1666. * scroll to residue in sequence windows when selected in structure window
  1667. *
  1668. * Revision 1.32  2001/12/06 23:13:44  thiessen
  1669. * finish import/align new sequences into single-structure data; many small tweaks
  1670. *
  1671. * Revision 1.31  2001/11/27 16:26:06  thiessen
  1672. * major update to data management system
  1673. *
  1674. * Revision 1.30  2001/10/02 23:26:47  thiessen
  1675. * fix bug in highlight columns in master range
  1676. *
  1677. * Revision 1.29  2001/09/27 15:37:57  thiessen
  1678. * decouple sequence import and BLAST
  1679. *
  1680. * Revision 1.28  2001/09/04 14:40:19  thiessen
  1681. * add rainbow and charge coloring
  1682. *
  1683. * Revision 1.27  2001/08/24 00:41:35  thiessen
  1684. * tweak conservation colors and opengl font handling
  1685. *
  1686. * Revision 1.26  2001/08/09 19:07:13  thiessen
  1687. * add temperature and hydrophobicity coloring
  1688. *
  1689. * Revision 1.25  2001/07/19 19:14:38  thiessen
  1690. * working CDD alignment annotator ; misc tweaks
  1691. *
  1692. * Revision 1.24  2001/07/10 16:39:54  thiessen
  1693. * change selection control keys; add CDD name/notes dialogs
  1694. *
  1695. * Revision 1.23  2001/06/21 02:02:33  thiessen
  1696. * major update to molecule identification and highlighting ; add toggle highlight (via alt)
  1697. *
  1698. * Revision 1.22  2001/06/02 17:22:45  thiessen
  1699. * fixes for GCC
  1700. *
  1701. * Revision 1.21  2001/06/01 13:35:58  thiessen
  1702. * add aligned block number to status line
  1703. *
  1704. * Revision 1.20  2001/05/31 18:47:06  thiessen
  1705. * add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
  1706. *
  1707. * Revision 1.19  2001/05/23 17:45:00  thiessen
  1708. * fix merge bug where unaligned block needs to be added
  1709. *
  1710. * Revision 1.18  2001/05/11 02:10:41  thiessen
  1711. * add better merge fail indicators; tweaks to windowing/taskbar
  1712. *
  1713. * Revision 1.17  2001/05/09 17:15:06  thiessen
  1714. * add automatic block removal upon demotion
  1715. *
  1716. * Revision 1.16  2001/05/02 13:46:27  thiessen
  1717. * major revision of stuff relating to saving of updates; allow stored null-alignments
  1718. *
  1719. * Revision 1.15  2001/04/19 14:24:05  thiessen
  1720. * fix row selection bug when alignments are of different size
  1721. *
  1722. * Revision 1.14  2001/04/19 12:58:32  thiessen
  1723. * allow merge and delete of individual updates
  1724. *
  1725. * Revision 1.13  2001/04/05 22:55:34  thiessen
  1726. * change bg color handling ; show geometry violations
  1727. *
  1728. * Revision 1.12  2001/04/04 00:27:14  thiessen
  1729. * major update - add merging, threader GUI controls
  1730. *
  1731. * Revision 1.11  2001/03/30 03:07:33  thiessen
  1732. * add threader score calculation & sorting
  1733. *
  1734. * Revision 1.10  2001/03/22 19:11:09  thiessen
  1735. * don't allow drag of gaps
  1736. *
  1737. * Revision 1.9  2001/03/22 00:33:16  thiessen
  1738. * initial threading working (PSSM only); free color storage in undo stack
  1739. *
  1740. * Revision 1.8  2001/03/09 15:49:03  thiessen
  1741. * major changes to add initial update viewer
  1742. *
  1743. * Revision 1.7  2001/03/06 20:20:50  thiessen
  1744. * progress towards >1 alignment in a SequenceDisplay ; misc minor fixes
  1745. *
  1746. * Revision 1.6  2001/03/01 20:15:50  thiessen
  1747. * major rearrangement of sequence viewer code into base and derived classes
  1748. *
  1749. * Revision 1.5  2001/02/13 20:33:49  thiessen
  1750. * add information content coloring
  1751. *
  1752. * Revision 1.4  2001/02/13 01:03:56  thiessen
  1753. * backward-compatible domain ID's in output; add ability to delete rows
  1754. *
  1755. * Revision 1.3  2001/02/02 20:17:33  thiessen
  1756. * can read in CDD with multi-structure but no struct. alignments
  1757. *
  1758. * Revision 1.2  2001/01/26 19:29:59  thiessen
  1759. * limit undo stack size ; fix memory leak
  1760. *
  1761. * Revision 1.1  2000/11/30 15:49:35  thiessen
  1762. * add show/hide rows; unpack sec. struc. and domain features
  1763. *
  1764. */