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

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

生物技术

开发平台：

C/C++

su_block_multiple_alignment.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: su_block_multiple_alignment.cpp,v $
  4.  * PRODUCTION Revision 1000.0  2004/06/01 18:14:28  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.9
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: su_block_multiple_alignment.cpp,v 1000.0 2004/06/01 18:14:28 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors:  Paul Thiessen
  36. *
  37. * File Description:
  38. *      Classes to hold alignment data
  39. *
  40. * ===========================================================================
  41. */
  43. #include <ncbi_pch.hpp>
  44. #include <corelib/ncbistd.hpp>
  45. #include <corelib/ncbi_limits.hpp>
  46. #include <util/tables/raw_scoremat.h>
  48. #include <objects/seqalign/Dense_diag.hpp>
  50. #include "su_block_multiple_alignment.hpp"
  51. #include "su_sequence_set.hpp"
  52. #include "su_private.hpp"
  54. // C-toolkit stuff for PSSM calculation
  55. #include <objalign.h>
  56. #include <blast.h>
  57. #include <blastkar.h>
  58. #include <cddutil.h>
  59. #include <thrdatd.h>
  60. #include <thrddecl.h>
  62. USING_NCBI_SCOPE;
  63. USING_SCOPE(objects);
  66. BEGIN_SCOPE(struct_util)
  68. // from su_sequence_set.cpp
  69. extern BioseqPtr GetOrCreateBioseq(const Sequence *sequence);
  70. extern void AddCSeqId(SeqIdPtr *sid, const ncbi::objects::CSeq_id& cppid);
  72. static const int SCALING_FACTOR = 1000000;
  73. static const string ThreaderResidues = "ARNDCQEGHILKMFPSTWYV";
  74. static const string BLASTResidues = "-ABCDEFGHIKLMNPQRSTVWXYZU*";
  76. BlockMultipleAlignment::BlockMultipleAlignment(const SequenceList& sequenceList)
  77. {
  78.     m_sequences = sequenceList;
  79.     m_pssm = NULL;
  81.     InitCache();
  82.     m_rowDoubles.resize(m_sequences.size(), 0.0);
  83.     m_rowStrings.resize(m_sequences.size());
  84. }
  86. void BlockMultipleAlignment::InitCache(void)
  87. {
  88.     m_cachePrevRow = eUndefined;
  89.     m_cachePrevBlock = NULL;
  90.     m_cacheBlockIterator = m_blocks.begin();
  91. }
  93. BlockMultipleAlignment::~BlockMultipleAlignment(void)
  94. {
  95.     RemovePSSM();
  96. }
  98. BlockMultipleAlignment * BlockMultipleAlignment::Clone(void) const
  99. {
  100.     BlockMultipleAlignment *copy = new BlockMultipleAlignment(m_sequences);
  101.     BlockList::const_iterator b, be = m_blocks.end();
  102.     for (b=m_blocks.begin(); b!=be; ++b)
  103.         copy->m_blocks.push_back(CRef<Block>((*b)->Clone(copy)));
  104.     copy->UpdateBlockMap();
  105.     copy->m_rowDoubles = m_rowDoubles;
  106.     copy->m_rowStrings = m_rowStrings;
  107.     return copy;
  108. }
  110. static inline unsigned int ScreenResidueCharacter(char original)
  111. {
  112.     char ch = toupper(original);
  113.     switch (ch) {
  114.         case 'A': case 'R': case 'N': case 'D': case 'C':
  115.         case 'Q': case 'E': case 'G': case 'H': case 'I':
  116.         case 'L': case 'K': case 'M': case 'F': case 'P':
  117.         case 'S': case 'T': case 'W': case 'Y': case 'V':
  118.         case 'B': case 'Z':
  119.             break;
  120.         default:
  121.             ch = 'X'; // make all but natural aa's just 'X'
  122.     }
  123.     return ch;
  124. }
  126. static int GetBLOSUM62Score(char a, char b)
  127. {
  128.     static SNCBIFullScoreMatrix Blosum62Matrix;
  129.     static bool unpacked = false;
  131.     if (!unpacked) {
  132.         NCBISM_Unpack(&NCBISM_Blosum62, &Blosum62Matrix);
  133.         unpacked = true;
  134.     }
  136.     return Blosum62Matrix.s[ScreenResidueCharacter(a)][ScreenResidueCharacter(b)];
  137. }
  139. // creates a SeqAlign from a BlockMultipleAlignment
  140. static SeqAlignPtr CreateCSeqAlign(const BlockMultipleAlignment& multiple)
  141. {
  142.     // one SeqAlign (chained into a linked list) for each slave row
  143.     SeqAlignPtr prevSap = NULL, firstSap = NULL;
  144.     for (unsigned int row=1; row<multiple.NRows(); ++row) {
  146.         SeqAlignPtr sap = SeqAlignNew();
  147.         if (prevSap) prevSap->next = sap;
  148.         prevSap = sap;
  149.         if (!firstSap) firstSap = sap;
  151.         sap->type = SAT_PARTIAL;
  152.         sap->dim = 2;
  153.         sap->segtype = SAS_DENDIAG;
  155.         DenseDiagPtr prevDd = NULL;
  156.         BlockMultipleAlignment::UngappedAlignedBlockList blocks;
  157.         multiple.GetUngappedAlignedBlocks(&blocks);
  158.         BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
  160.         for (b=blocks.begin(); b!=be; ++b) {
  161.             DenseDiagPtr dd = DenseDiagNew();
  162.             if (prevDd) prevDd->next = dd;
  163.             prevDd = dd;
  164.             if (b == blocks.begin()) sap->segs = dd;
  166.             dd->dim = 2;
  167.             AddCSeqId(&(dd->id), multiple.GetSequenceOfRow(0)->GetPreferredIdentifier());
  168.             AddCSeqId(&(dd->id), multiple.GetSequenceOfRow(row)->GetPreferredIdentifier());
  170.             dd->len = (*b)->m_width;
  172.             dd->starts = (Int4Ptr) MemNew(2 * sizeof(Int4));
  173.             const Block::Range *range = (*b)->GetRangeOfRow(0);
  174.             dd->starts[0] = range->from;
  175.             range = (*b)->GetRangeOfRow(row);
  176.             dd->starts[1] = range->from;
  177.         }
  178.     }
  180.     return firstSap;
  181. }
  183. static void CreateAllBioseqs(const BlockMultipleAlignment& multiple)
  184. {
  185.     for (unsigned int row=0; row<multiple.NRows(); ++row)
  186.         GetOrCreateBioseq(multiple.GetSequenceOfRow(row));
  187. }
  189. static Seq_Mtf * CreateSeqMtf(const BlockMultipleAlignment& multiple,
  190.     double weightPSSM, BLAST_KarlinBlkPtr karlinBlock)
  191. {
  192.     // special case for "PSSM" of single-row "alignment" - just use BLOSUM62 score
  193.     if (multiple.NRows() == 1) {
  194.         Seq_Mtf *seqMtf = NewSeqMtf(multiple.GetMaster()->Length(), ThreaderResidues.size());
  195.         for (unsigned int res=0; res<multiple.GetMaster()->Length(); ++res)
  196.             for (unsigned int aa=0; aa<ThreaderResidues.size(); ++aa)
  197.                 seqMtf->ww[res][aa] = ThrdRound(weightPSSM * SCALING_FACTOR *
  198.                     GetBLOSUM62Score(multiple.GetMaster()->m_sequenceString[res], ThreaderResidues[aa]));
  199.         WARNING_MESSAGE("Created Seq_Mtf (PSSM) from BLOSUM62 scores");
  200.         return seqMtf;
  201.     }
  203.     // convert all sequences to Bioseqs
  204.     CreateAllBioseqs(multiple);
  206.     // create SeqAlign from this BlockMultipleAlignment
  207.     SeqAlignPtr seqAlign = CreateCSeqAlign(multiple);
  209.     // "spread" unaligned residues between aligned blocks, for PSSM construction
  210.     CddDegapSeqAlign(seqAlign);
  212.     Seq_Mtf *seqMtf = NULL;
  213.     for (int i=11; i>=1; --i) {
  214.         // first try auto-determined pseudocount (-1); if fails, find higest <= 10 that works
  215.         int pseudocount = (i == 11) ? -1 : i;
  216.         seqMtf = CddDenDiagCposComp2KBP(
  217.             GetOrCreateBioseq(multiple.GetMaster()),
  218.             pseudocount,
  219.             seqAlign,
  220.             NULL,
  221.             NULL,
  222.             weightPSSM,
  223.             SCALING_FACTOR,
  224.             NULL,
  225.             karlinBlock
  226.         );
  227.         if (seqMtf)
  228.             break;
  229.         else
  230.             WARNING_MESSAGE("Cannot use " << ((pseudocount == -1) ? "(empirical) " : "")
  231.                 << "pseudocount of " << pseudocount);
  232.     }
  234.     if (seqMtf)
  235.         TRACE_MESSAGE("created Seq_Mtf (PSSM)");
  236.     else
  237.         ERROR_MESSAGE("Cannot find any pseudocount that yields an acceptable PSSM!");
  239.     SeqAlignSetFree(seqAlign);
  240.     return seqMtf;
  241. }
  243. // gives BLAST residue number for a character (or # for 'X' if char not found)
  244. int LookupBLASTResidueNumberFromCharacter(unsigned char r)
  245. {
  246.     typedef map < unsigned char, int > Char2Int;
  247.     static Char2Int charMap;
  249.     if (charMap.size() == 0) {
  250.         for (unsigned int i=0; i<BLASTResidues.size(); ++i)
  251.             charMap[BLASTResidues[i]] = i;
  252.     }
  254.     Char2Int::const_iterator n = charMap.find(toupper(r));
  255.     if (n != charMap.end())
  256.         return n->second;
  257.     else
  258.         return charMap.find('X')->second;
  259. }
  261. // gives BLAST residue number for a threader residue number (or # for 'X' if char == -1)
  262. static int LookupBLASTResidueNumberFromThreaderResidueNumber(unsigned char r)
  263. {
  264.     r = toupper(r);
  265.     return LookupBLASTResidueNumberFromCharacter(
  266.             (r < ThreaderResidues.size()) ? ThreaderResidues[r] : 'X');
  267. }
  269. static Int4 Round(double d)
  270. {
  271.     if (d >= 0.0)
  272.         return (Int4) (d + 0.5);
  273.     else
  274.         return (Int4) (d - 0.5);
  275. }
  277. const BLAST_Matrix * BlockMultipleAlignment::GetPSSM(void) const
  278. {
  279.     if (m_pssm) return m_pssm;
  281.     // for now, use threader's SeqMtf
  282.     BLAST_KarlinBlkPtr karlinBlock = BlastKarlinBlkCreate();
  283.     Seq_Mtf *seqMtf = CreateSeqMtf(*this, 1.0, karlinBlock);
  285.     m_pssm = (BLAST_Matrix *) MemNew(sizeof(BLAST_Matrix));
  286.     m_pssm->is_prot = TRUE;
  287.     m_pssm->name = StringSave("BLOSUM62");
  288.     m_pssm->karlinK = karlinBlock->K;
  289.     m_pssm->rows = seqMtf->n + 1;
  290.     m_pssm->columns = 26;
  292.     int i, j;
  293.     m_pssm->matrix = (Int4 **) MemNew(m_pssm->rows * sizeof(Int4 *));
  294.     for (i=0; i<m_pssm->rows; ++i) {
  295.         m_pssm->matrix[i] = (Int4 *) MemNew(m_pssm->columns * sizeof(Int4));
  297.         // set scores from threader matrix
  298.         if (i < seqMtf->n) {
  299.             // initialize all rows with custom score, or BLAST_SCORE_MIN; to match what Aron's function creates
  300.             for (j=0; j<m_pssm->columns; ++j)
  301.                 m_pssm->matrix[i][j] = (j == 21 ? -1 : (j == 25 ? -4 : BLAST_SCORE_MIN));
  303.             for (j=0; j<seqMtf->AlphabetSize; ++j) {
  304.                 m_pssm->matrix[i][LookupBLASTResidueNumberFromThreaderResidueNumber(j)] =
  305.                     Round(((double) seqMtf->ww[i][j]) / SCALING_FACTOR);
  306.             }
  307.         } else {
  308.             // initialize last row with BLAST_SCORE_MIN
  309.             for (j=0; j<m_pssm->columns; ++j)
  310.                 m_pssm->matrix[i][j] = BLAST_SCORE_MIN;
  311.         }
  312.     }
  314.     m_pssm->posFreqs = NULL;
  316.     FreeSeqMtf(seqMtf);
  317.     BlastKarlinBlkDestruct(karlinBlock);
  318.     return m_pssm;
  319. }
  321. void BlockMultipleAlignment::RemovePSSM(void) const
  322. {
  323.     if (m_pssm) {
  324.         BLAST_MatrixDestruct(m_pssm);
  325.         m_pssm = NULL;
  326.     }
  327. }
  329. bool BlockMultipleAlignment::CheckAlignedBlock(const Block *block) const
  330. {
  331.     if (!block || !block->IsAligned()) {
  332.         ERROR_MESSAGE("CheckAlignedBlock() - checks aligned blocks only");
  333.         return false;
  334.     }
  335.     if (block->NSequences() != m_sequences.size()) {
  336.         ERROR_MESSAGE("CheckAlignedBlock() - block size mismatch");
  337.         return false;
  338.     }
  340.     // make sure ranges are reasonable for each sequence
  341.     unsigned int row;
  342.     const Block
  343.         *prevBlock = GetBlockBefore(block),
  344.         *nextBlock = GetBlockAfter(block);
  345.     const Block::Range *range, *prevRange = NULL, *nextRange = NULL;
  346.     SequenceList::const_iterator sequence = m_sequences.begin();
  347.     for (row=0; row<block->NSequences(); ++row, ++sequence) {
  348.         range = block->GetRangeOfRow(row);
  349.         if (prevBlock) prevRange = prevBlock->GetRangeOfRow(row);
  350.         if (nextBlock) nextRange = nextBlock->GetRangeOfRow(row);
  351.         if (range->to - range->from + 1 != block->m_width ||       // check block m_width
  352.             (prevRange && range->from <= prevRange->to) ||          // check for range overlap
  353.             (nextRange && range->to >= nextRange->from) ||          // check for range overlap
  354.             range->from > range->to ||                              // check range values
  355.             range->to >= (*sequence)->Length())                     // check bounds of end
  356.         {
  357.             ERROR_MESSAGE("CheckAlignedBlock() - range error");
  358.             return false;
  359.         }
  360.     }
  362.     return true;
  363. }
  365. bool BlockMultipleAlignment::AddAlignedBlockAtEnd(UngappedAlignedBlock *newBlock)
  366. {
  367.     m_blocks.push_back(CRef<Block>(newBlock));
  368.     return CheckAlignedBlock(newBlock);
  369. }
  371. UnalignedBlock * BlockMultipleAlignment::
  372.     CreateNewUnalignedBlockBetween(const Block *leftBlock, const Block *rightBlock)
  373. {
  374.     if ((leftBlock && !leftBlock->IsAligned()) ||
  375.         (rightBlock && !rightBlock->IsAligned())) {
  376.         ERROR_MESSAGE("CreateNewUnalignedBlockBetween() - passed an unaligned block");
  377.         return NULL;
  378.     }
  380.     unsigned int row, from, to, length;
  381.     SequenceList::const_iterator s, se = m_sequences.end();
  383.     UnalignedBlock *newBlock = new UnalignedBlock(this);
  384.     newBlock->m_width = 0;
  386.     for (row=0, s=m_sequences.begin(); s!=se; ++row, ++s) {
  388.         if (leftBlock)
  389.             from = leftBlock->GetRangeOfRow(row)->to + 1;
  390.         else
  391.             from = 0;
  393.         if (rightBlock)
  394.             to = rightBlock->GetRangeOfRow(row)->from - 1;
  395.         else
  396.             to = (*s)->Length() - 1;
  398.         newBlock->SetRangeOfRow(row, from, to);
  400.         length = to - from + 1;
  401.         if ((to - from + 1) < 0) { // just to make sure...
  402.             ERROR_MESSAGE("CreateNewUnalignedBlockBetween() - unaligned length < 0");
  403.             return NULL;
  404.         }
  405.         if (length > newBlock->m_width)
  406.             newBlock->m_width = length;
  407.     }
  409.     if (newBlock->m_width == 0) {
  410.         delete newBlock;
  411.         return NULL;
  412.     } else
  413.         return newBlock;
  414. }
  416. bool BlockMultipleAlignment::AddUnalignedBlocks(void)
  417. {
  418.     BlockList::iterator a, ae = m_blocks.end();
  419.     const Block *alignedBlock = NULL, *prevAlignedBlock = NULL;
  420.     Block *newUnalignedBlock;
  422.     // unaligned m_blocks to the left of each aligned block
  423.     for (a=m_blocks.begin(); a!=ae; ++a) {
  424.         alignedBlock = *a;
  425.         newUnalignedBlock = CreateNewUnalignedBlockBetween(prevAlignedBlock, alignedBlock);
  426.         if (newUnalignedBlock)
  427.             m_blocks.insert(a, CRef<Block>(newUnalignedBlock));
  428.         prevAlignedBlock = alignedBlock;
  429.     }
  431.     // right tail
  432.     newUnalignedBlock = CreateNewUnalignedBlockBetween(alignedBlock, NULL);
  433.     if (newUnalignedBlock) {
  434.         m_blocks.insert(a, CRef<Block>(newUnalignedBlock));
  435.     }
  437.     return true;
  438. }
  440. bool BlockMultipleAlignment::UpdateBlockMap(bool clearRowInfo)
  441. {
  442.     unsigned int i = 0, j, n = 0;
  443.     BlockList::iterator b, be = m_blocks.end();
  445.     // reset old stuff, recalculate m_width
  446.     m_totalWidth = 0;
  447.     for (b=m_blocks.begin(); b!=be; ++b)
  448.         m_totalWidth += (*b)->m_width;
  450.     // fill out the block map
  451.     m_blockMap.resize(m_totalWidth);
  452.     UngappedAlignedBlock *aBlock;
  453.     for (b=m_blocks.begin(); b!=be; ++b) {
  454.         aBlock = dynamic_cast<UngappedAlignedBlock*>(b->GetPointer());
  455.         if (aBlock)
  456.             ++n;
  457.         for (j=0; j<(*b)->m_width; ++j, ++i) {
  458.             m_blockMap[i].block = *b;
  459.             m_blockMap[i].blockColumn = j;
  460.             m_blockMap[i].alignedBlockNum = aBlock ? n : eUndefined;
  461.         }
  462.     }
  464.     // if alignment changes, any pssm/scores/status become invalid
  465.     RemovePSSM();
  466.     if (clearRowInfo) ClearRowInfo();
  468.     return true;
  469. }
  471. bool BlockMultipleAlignment::GetCharacterAt(
  472.     unsigned int alignmentColumn, unsigned int row, eUnalignedJustification justification,
  473.     char *character) const
  474. {
  475.     const Sequence *sequence;
  476.     unsigned int seqIndex;
  477.     bool isAligned;
  479.     if (!GetSequenceAndIndexAt(alignmentColumn, row, justification, &sequence, &seqIndex, &isAligned))
  480.         return false;
  482.     *character = (seqIndex >= 0) ? sequence->m_sequenceString[seqIndex] : '~';
  483.     if (isAligned)
  484.         *character = toupper(*character);
  485.     else
  486.         *character = tolower(*character);
  488.     return true;
  489. }
  491. bool BlockMultipleAlignment::GetSequenceAndIndexAt(
  492.     unsigned int alignmentColumn, unsigned int row, eUnalignedJustification requestedJustification,
  493.     const Sequence **sequence, unsigned int *index, bool *isAligned) const
  494. {
  495.     if (sequence)
  496.         *sequence = m_sequences[row];
  498.     const BlockInfo& blockInfo = m_blockMap[alignmentColumn];
  500.     if (!blockInfo.block->IsAligned()) {
  501.         if (isAligned)
  502.             *isAligned = false;
  503.         // override requested justification for end m_blocks
  504.         if (blockInfo.block == m_blocks.back().GetPointer()) // also true if there's a single aligned block
  505.             requestedJustification = eLeft;
  506.         else if (blockInfo.block == m_blocks.front().GetPointer())
  507.             requestedJustification = eRight;
  508.     } else
  509.         if (isAligned)
  510.             *isAligned = true;
  512.     if (index)
  513.         *index = blockInfo.block->GetIndexAt(blockInfo.blockColumn, row, requestedJustification);
  515.     return true;
  516. }
  518. bool BlockMultipleAlignment::IsAligned(unsigned int row, unsigned int seqIndex) const
  519. {
  520.     const Block *block = GetBlock(row, seqIndex);
  521.     return (block && block->IsAligned());
  522. }
  524. const Block * BlockMultipleAlignment::GetBlock(unsigned int row, unsigned int seqIndex) const
  525. {
  526.     // make sure we're in range for this sequence
  527.     if (row < 0 || seqIndex < 0 || row >= NRows() || seqIndex >= m_sequences[row]->Length()) {
  528.         ERROR_MESSAGE("BlockMultipleAlignment::GetBlock() - coordinate out of range");
  529.         return NULL;
  530.     }
  532.     const Block::Range *range;
  534.     // first check to see if it's in the same block as last time.
  535.     if (m_cachePrevBlock) {
  536.         range = m_cachePrevBlock->GetRangeOfRow(row);
  537.         if (seqIndex >= range->from && seqIndex <= range->to)
  538.             return m_cachePrevBlock;
  539.         ++m_cacheBlockIterator; // start search at next block
  540.     } else {
  541.         m_cacheBlockIterator = m_blocks.begin();
  542.     }
  544.     // otherwise, perform block search. This search is most efficient when queries
  545.     // happen in order from left to right along a given row.
  546.     do {
  547.         if (m_cacheBlockIterator == m_blocks.end())
  548.             m_cacheBlockIterator = m_blocks.begin();
  549.         range = (*m_cacheBlockIterator)->GetRangeOfRow(row);
  550.         if (seqIndex >= range->from && seqIndex <= range->to) {
  551.             m_cachePrevBlock = *m_cacheBlockIterator; // cache this block
  552.             return m_cachePrevBlock;
  553.         }
  554.         ++m_cacheBlockIterator;
  555.     } while (1);
  556. }
  558. unsigned int BlockMultipleAlignment::GetFirstAlignedBlockPosition(void) const
  559. {
  560.     BlockList::const_iterator b = m_blocks.begin();
  561.     if (m_blocks.size() > 0 && (*b)->IsAligned()) // first block is aligned
  562.         return 0;
  563.     else if (m_blocks.size() >= 2 && (*(++b))->IsAligned()) // second block is aligned
  564.         return m_blocks.front()->m_width;
  565.     else
  566.         return eUndefined;
  567. }
  569. unsigned int BlockMultipleAlignment::GetAlignedSlaveIndex(unsigned int masterSeqIndex, unsigned int slaveRow) const
  570. {
  571.     const UngappedAlignedBlock
  572.         *aBlock = dynamic_cast<const UngappedAlignedBlock*>(GetBlock(0, masterSeqIndex));
  573.     if (!aBlock)
  574.         return eUndefined;
  576.     const Block::Range
  577.         *masterRange = aBlock->GetRangeOfRow(0),
  578.         *slaveRange = aBlock->GetRangeOfRow(slaveRow);
  579.     return (slaveRange->from + masterSeqIndex - masterRange->from);
  580. }
  582. void BlockMultipleAlignment::GetAlignedBlockPosition(unsigned int alignmentIndex,
  583.     unsigned int *blockColumn, unsigned int *blockWidth) const
  584. {
  585.     *blockColumn = *blockWidth = eUndefined;
  586.     const BlockInfo& info = m_blockMap[alignmentIndex];
  587.     if (info.block->IsAligned()) {
  588.         *blockColumn = info.blockColumn;
  589.         *blockWidth = info.block->m_width;
  590.     }
  591. }
  593. Block * BlockMultipleAlignment::GetBlockBefore(const Block *block)
  594. {
  595.     Block *prevBlock = NULL;
  596.     BlockList::iterator b, be = m_blocks.end();
  597.     for (b=m_blocks.begin(); b!=be; ++b) {
  598.         if (*b == block) break;
  599.         prevBlock = b->GetPointer();
  600.     }
  601.     return prevBlock;
  602. }
  604. const Block * BlockMultipleAlignment::GetBlockBefore(const Block *block) const
  605. {
  606.     const Block *prevBlock = NULL;
  607.     BlockList::const_iterator b, be = m_blocks.end();
  608.     for (b=m_blocks.begin(); b!=be; ++b) {
  609.         if (*b == block) break;
  610.         prevBlock = b->GetPointer();
  611.     }
  612.     return prevBlock;
  613. }
  615. Block * BlockMultipleAlignment::GetBlockAfter(const Block *block)
  616. {
  617.     BlockList::iterator b, be = m_blocks.end();
  618.     for (b=m_blocks.begin(); b!=be; ++b) {
  619.         if (*b == block) {
  620.             ++b;
  621.             if (b == be) break;
  622.             return *b;
  623.         }
  624.     }
  625.     return NULL;
  626. }
  628. const Block * BlockMultipleAlignment::GetBlockAfter(const Block *block) const
  629. {
  630.     BlockList::const_iterator b, be = m_blocks.end();
  631.     for (b=m_blocks.begin(); b!=be; ++b) {
  632.         if (*b == block) {
  633.             ++b;
  634.             if (b == be) break;
  635.             return *b;
  636.         }
  637.     }
  638.     return NULL;
  639. }
  641. const UnalignedBlock * BlockMultipleAlignment::GetUnalignedBlockBefore(
  642.     const UngappedAlignedBlock *aBlock) const
  643. {
  644.     const Block *prevBlock;
  645.     if (aBlock)
  646.         prevBlock = GetBlockBefore(aBlock);
  647.     else
  648.         prevBlock = m_blocks.back().GetPointer();
  649.     return dynamic_cast<const UnalignedBlock*>(prevBlock);
  650. }
  652. void BlockMultipleAlignment::InsertBlockBefore(Block *newBlock, const Block *insertAt)
  653. {
  654.     BlockList::iterator b, be = m_blocks.end();
  655.     for (b=m_blocks.begin(); b!=be; ++b) {
  656.         if (*b == insertAt) {
  657.             m_blocks.insert(b, CRef<Block>(newBlock));
  658.             return;
  659.         }
  660.     }
  661.     WARNING_MESSAGE("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
  662. }
  664. void BlockMultipleAlignment::InsertBlockAfter(const Block *insertAt, Block *newBlock)
  665. {
  666.     BlockList::iterator b, be = m_blocks.end();
  667.     for (b=m_blocks.begin(); b!=be; ++b) {
  668.         if (*b == insertAt) {
  669.             ++b;
  670.             m_blocks.insert(b, CRef<Block>(newBlock));
  671.             return;
  672.         }
  673.     }
  674.     WARNING_MESSAGE("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
  675. }
  677. void BlockMultipleAlignment::RemoveBlock(Block *block)
  678. {
  679.     BlockList::iterator b, be = m_blocks.end();
  680.     for (b=m_blocks.begin(); b!=be; ++b) {
  681.         if (*b == block) {
  682.             m_blocks.erase(b);
  683.             InitCache();
  684.             return;
  685.         }
  686.     }
  687.     WARNING_MESSAGE("BlockMultipleAlignment::RemoveBlock() - couldn't find block");
  688. }
  690. bool BlockMultipleAlignment::MoveBlockBoundary(unsigned int columnFrom, unsigned int columnTo)
  691. {
  692.     unsigned int blockColumn, blockWidth;
  693.     GetAlignedBlockPosition(columnFrom, &blockColumn, &blockWidth);
  694.     if (blockColumn == eUndefined || blockWidth == 0 || blockWidth == eUndefined) return false;
  696.     TRACE_MESSAGE("trying to move block boundary from " << columnFrom << " to " << columnTo);
  698.     const BlockInfo& info = m_blockMap[columnFrom];
  699.     unsigned int row;
  700.     int requestedShift = columnTo - columnFrom, actualShift = 0;
  701.     const Block::Range *range;
  703.     // shrink block from left
  704.     if (blockColumn == 0 && requestedShift > 0 && requestedShift < (int) info.block->m_width) {
  705.         actualShift = requestedShift;
  706.         TRACE_MESSAGE("shrinking block from left");
  707.         for (row=0; row<NRows(); ++row) {
  708.             range = info.block->GetRangeOfRow(row);
  709.             info.block->SetRangeOfRow(row, range->from + requestedShift, range->to);
  710.         }
  711.         info.block->m_width -= requestedShift;
  712.         Block *prevBlock = GetBlockBefore(info.block);
  713.         if (prevBlock && !prevBlock->IsAligned()) {
  714.             for (row=0; row<NRows(); ++row) {
  715.                 range = prevBlock->GetRangeOfRow(row);
  716.                 prevBlock->SetRangeOfRow(row, range->from, range->to + requestedShift);
  717.             }
  718.             prevBlock->m_width += requestedShift;
  719.         } else {
  720.             Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(prevBlock, info.block);
  721.             if (newUnalignedBlock)
  722.                 InsertBlockBefore(newUnalignedBlock, info.block);
  723.             TRACE_MESSAGE("added new unaligned block");
  724.         }
  725.     }
  727.     // shrink block from right (requestedShift < 0)
  728.     else if (blockColumn == info.block->m_width - 1 &&
  729.                 requestedShift < 0 && ((unsigned int) -requestedShift) < info.block->m_width) {
  730.         actualShift = requestedShift;
  731.         TRACE_MESSAGE("shrinking block from right");
  732.         for (row=0; row<NRows(); ++row) {
  733.             range = info.block->GetRangeOfRow(row);
  734.             info.block->SetRangeOfRow(row, range->from, range->to + requestedShift);
  735.         }
  736.         info.block->m_width += requestedShift;
  737.         Block *nextBlock = GetBlockAfter(info.block);
  738.         if (nextBlock && !nextBlock->IsAligned()) {
  739.             for (row=0; row<NRows(); ++row) {
  740.                 range = nextBlock->GetRangeOfRow(row);
  741.                 nextBlock->SetRangeOfRow(row, range->from + requestedShift, range->to);
  742.             }
  743.             nextBlock->m_width -= requestedShift;
  744.         } else {
  745.             Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(info.block, nextBlock);
  746.             if (newUnalignedBlock)
  747.                 InsertBlockAfter(info.block, newUnalignedBlock);
  748.             TRACE_MESSAGE("added new unaligned block");
  749.         }
  750.     }
  752.     // grow block to right
  753.     else if (blockColumn == info.block->m_width - 1 && requestedShift > 0) {
  754.         Block *nextBlock = GetBlockAfter(info.block);
  755.         if (nextBlock && !nextBlock->IsAligned()) {
  756.             int nRes;
  757.             actualShift = requestedShift;
  758.             for (row=0; row<NRows(); ++row) {
  759.                 range = nextBlock->GetRangeOfRow(row);
  760.                 nRes = range->to - range->from + 1;
  761.                 if (nRes < actualShift)
  762.                     actualShift = nRes;
  763.             }
  764.             if (actualShift) {
  765.                 TRACE_MESSAGE("growing block to right");
  766.                 for (row=0; row<NRows(); ++row) {
  767.                     range = info.block->GetRangeOfRow(row);
  768.                     info.block->SetRangeOfRow(row, range->from, range->to + actualShift);
  769.                     range = nextBlock->GetRangeOfRow(row);
  770.                     nextBlock->SetRangeOfRow(row, range->from + actualShift, range->to);
  771.                 }
  772.                 info.block->m_width += actualShift;
  773.                 nextBlock->m_width -= actualShift;
  774.                 if (nextBlock->m_width == 0) {
  775.                     RemoveBlock(nextBlock);
  776.                     TRACE_MESSAGE("removed empty block");
  777.                 }
  778.             }
  779.         }
  780.     }
  782.     // grow block to left (requestedShift < 0)
  783.     else if (blockColumn == 0 && requestedShift < 0) {
  784.         Block *prevBlock = GetBlockBefore(info.block);
  785.         if (prevBlock && !prevBlock->IsAligned()) {
  786.             int nRes;
  787.             actualShift = requestedShift;
  788.             for (row=0; row<NRows(); ++row) {
  789.                 range = prevBlock->GetRangeOfRow(row);
  790.                 nRes = range->to - range->from + 1;
  791.                 if (nRes < -actualShift) actualShift = -nRes;
  792.             }
  793.             if (actualShift) {
  794.                 TRACE_MESSAGE("growing block to left");
  795.                 for (row=0; row<NRows(); ++row) {
  796.                     range = info.block->GetRangeOfRow(row);
  797.                     info.block->SetRangeOfRow(row, range->from + actualShift, range->to);
  798.                     range = prevBlock->GetRangeOfRow(row);
  799.                     prevBlock->SetRangeOfRow(row, range->from, range->to + actualShift);
  800.                 }
  801.                 info.block->m_width -= actualShift;
  802.                 prevBlock->m_width += actualShift;
  803.                 if (prevBlock->m_width == 0) {
  804.                     RemoveBlock(prevBlock);
  805.                     TRACE_MESSAGE("removed empty block");
  806.                 }
  807.             }
  808.         }
  809.     }
  811.     if (actualShift != 0) {
  812.         UpdateBlockMap();
  813.         return true;
  814.     } else
  815.         return false;
  816. }
  818. bool BlockMultipleAlignment::ShiftRow(unsigned int row, unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex,
  819.     eUnalignedJustification justification)
  820. {
  821.     if (fromAlignmentIndex == toAlignmentIndex)
  822.         return false;
  824.     Block
  825.         *blockFrom = m_blockMap[fromAlignmentIndex].block,
  826.         *blockTo = m_blockMap[toAlignmentIndex].block;
  828.     // at least one end of the drag must be in an aligned block
  829.     UngappedAlignedBlock *ABlock =
  830.         dynamic_cast<UngappedAlignedBlock*>(blockFrom);
  831.     if (ABlock) {
  832.         if (blockTo != blockFrom && blockTo->IsAligned())
  833.             return false;
  834.     } else {
  835.         ABlock = dynamic_cast<UngappedAlignedBlock*>(blockTo);
  836.         if (!ABlock)
  837.             return false;
  838.     }
  840.     // and the other must be in the same aligned block or an adjacent unaligned block
  841.     UnalignedBlock
  842.         *prevUABlock = dynamic_cast<UnalignedBlock*>(GetBlockBefore(ABlock)),
  843.         *nextUABlock = dynamic_cast<UnalignedBlock*>(GetBlockAfter(ABlock));
  844.     if (blockFrom != blockTo &&
  845.         ((ABlock == blockFrom && prevUABlock != blockTo && nextUABlock != blockTo) ||
  846.          (ABlock == blockTo && prevUABlock != blockFrom && nextUABlock != blockFrom)))
  847.         return false;
  849.     int requestedShift, actualShift = 0, width = 0;
  851.     // slightly different behaviour when dragging from unaligned to aligned...
  852.     if (!blockFrom->IsAligned()) {
  853.         unsigned int fromSeqIndex, toSeqIndex;
  854.         GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, NULL, &fromSeqIndex, NULL);
  855.         GetSequenceAndIndexAt(toAlignmentIndex, row, justification, NULL, &toSeqIndex, NULL);
  856.         if (fromSeqIndex == eUndefined || toSeqIndex == eUndefined)
  857.             return false;
  858.         requestedShift = toSeqIndex - fromSeqIndex;
  859.     }
  861.     // vs. dragging from aligned
  862.     else {
  863.         requestedShift = toAlignmentIndex - fromAlignmentIndex;
  864.     }
  866.     const Block::Range *prevRange = NULL, *nextRange = NULL, *range = ABlock->GetRangeOfRow(row);
  867.     if (prevUABlock) prevRange = prevUABlock->GetRangeOfRow(row);
  868.     if (nextUABlock) nextRange = nextUABlock->GetRangeOfRow(row);
  869.     if (requestedShift > 0) {
  870.         if (prevUABlock)
  871.             width = prevRange->to - prevRange->from + 1;
  872.         actualShift = (width > requestedShift) ? requestedShift : width;
  873.     } else {
  874.         if (nextUABlock)
  875.             width = nextRange->to - nextRange->from + 1;
  876.         actualShift = (width > -requestedShift) ? requestedShift : -width;
  877.     }
  878.     if (actualShift == 0) return false;
  880.     TRACE_MESSAGE("shifting row " << row << " by " << actualShift);
  882.     ABlock->SetRangeOfRow(row, range->from - actualShift, range->to - actualShift);
  884.     if (prevUABlock) {
  885.         prevUABlock->SetRangeOfRow(row, prevRange->from, prevRange->to - actualShift);
  886.         prevUABlock->m_width = 0;
  887.         for (unsigned int i=0; i<NRows(); ++i) {
  888.             prevRange = prevUABlock->GetRangeOfRow(i);
  889.             width = prevRange->to - prevRange->from + 1;
  890.             if (width < 0)
  891.                 ERROR_MESSAGE("BlockMultipleAlignment::ShiftRow() - negative width on left");
  892.             if ((unsigned int) width > prevUABlock->m_width)
  893.                 prevUABlock->m_width = width;
  894.         }
  895.         if (prevUABlock->m_width == 0) {
  896.             TRACE_MESSAGE("removing zero-m_width unaligned block on left");
  897.             RemoveBlock(prevUABlock);
  898.         }
  899.     } else {
  900.         TRACE_MESSAGE("creating unaligned block on left");
  901.         prevUABlock = CreateNewUnalignedBlockBetween(GetBlockBefore(ABlock), ABlock);
  902.         InsertBlockBefore(prevUABlock, ABlock);
  903.     }
  905.     if (nextUABlock) {
  906.         nextUABlock->SetRangeOfRow(row, nextRange->from - actualShift, nextRange->to);
  907.         nextUABlock->m_width = 0;
  908.         for (unsigned int i=0; i<NRows(); ++i) {
  909.             nextRange = nextUABlock->GetRangeOfRow(i);
  910.             width = nextRange->to - nextRange->from + 1;
  911.             if (width < 0)
  912.                 ERROR_MESSAGE("BlockMultipleAlignment::ShiftRow() - negative width on right");
  913.             if ((unsigned int) width > nextUABlock->m_width)
  914.                 nextUABlock->m_width = width;
  915.         }
  916.         if (nextUABlock->m_width == 0) {
  917.             TRACE_MESSAGE("removing zero-m_width unaligned block on right");
  918.             RemoveBlock(nextUABlock);
  919.         }
  920.     } else {
  921.         TRACE_MESSAGE("creating unaligned block on right");
  922.         nextUABlock = CreateNewUnalignedBlockBetween(ABlock, GetBlockAfter(ABlock));
  923.         InsertBlockAfter(ABlock, nextUABlock);
  924.     }
  926.     if (!CheckAlignedBlock(ABlock))
  927.         ERROR_MESSAGE("BlockMultipleAlignment::ShiftRow() - shift failed to create valid aligned block");
  929.     UpdateBlockMap();
  930.     return true;
  931. }
  933. bool BlockMultipleAlignment::SplitBlock(unsigned int alignmentIndex)
  934. {
  935.     const BlockInfo& info = m_blockMap[alignmentIndex];
  936.     if (!info.block->IsAligned() || info.block->m_width < 2 || info.blockColumn == 0)
  937.         return false;
  939.     TRACE_MESSAGE("splitting block");
  940.     UngappedAlignedBlock *newAlignedBlock = new UngappedAlignedBlock(this);
  941.     newAlignedBlock->m_width = info.block->m_width - info.blockColumn;
  942.     info.block->m_width = info.blockColumn;
  944.     const Block::Range *range;
  945.     unsigned int oldTo;
  946.     for (unsigned int row=0; row<NRows(); ++row) {
  947.         range = info.block->GetRangeOfRow(row);
  948.         oldTo = range->to;
  949.         info.block->SetRangeOfRow(row, range->from, range->from + info.block->m_width - 1);
  950.         newAlignedBlock->SetRangeOfRow(row, oldTo - newAlignedBlock->m_width + 1, oldTo);
  951.     }
  953.     InsertBlockAfter(info.block, newAlignedBlock);
  954.     if (!CheckAlignedBlock(info.block) || !CheckAlignedBlock(newAlignedBlock))
  955.         ERROR_MESSAGE("BlockMultipleAlignment::SplitBlock() - split failed to create valid m_blocks");
  957.     UpdateBlockMap();
  958.     return true;
  959. }
  961. bool BlockMultipleAlignment::MergeBlocks(unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex)
  962. {
  963.     Block
  964.         *expandedBlock = m_blockMap[fromAlignmentIndex].block,
  965.         *lastBlock = m_blockMap[toAlignmentIndex].block;
  966.     if (expandedBlock == lastBlock)
  967.         return false;
  968.     unsigned int i;
  969.     for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i)
  970.         if (!m_blockMap[i].block->IsAligned())
  971.             return false;
  973.     TRACE_MESSAGE("merging block(s)");
  974.     for (i=0; i<NRows(); ++i)
  975.         expandedBlock->SetRangeOfRow(i, expandedBlock->GetRangeOfRow(i)->from, lastBlock->GetRangeOfRow(i)->to);
  976.     expandedBlock->m_width = lastBlock->GetRangeOfRow(0)->to - expandedBlock->GetRangeOfRow(0)->from + 1;
  978.     Block *deletedBlock = NULL, *blockToDelete;
  979.     for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i) {
  980.         blockToDelete = m_blockMap[i].block;
  981.         if (blockToDelete == expandedBlock)
  982.             continue;
  983.         if (blockToDelete != deletedBlock) {
  984.             deletedBlock = blockToDelete;
  985.             RemoveBlock(blockToDelete);
  986.         }
  987.     }
  989.     if (!CheckAlignedBlock(expandedBlock))
  990.         ERROR_MESSAGE("BlockMultipleAlignment::MergeBlocks() - merge failed to create valid block");
  992.     UpdateBlockMap();
  993.     return true;
  994. }
  996. bool BlockMultipleAlignment::CreateBlock(unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex,
  997.     eUnalignedJustification justification)
  998. {
  999.     const BlockInfo& info = m_blockMap[fromAlignmentIndex];
  1000.     UnalignedBlock *prevUABlock = dynamic_cast<UnalignedBlock*>(info.block);
  1001.     if (!prevUABlock || info.block != m_blockMap[toAlignmentIndex].block)
  1002.         return false;
  1003.     unsigned int row, seqIndexFrom, seqIndexTo,
  1004.         newBlockWidth = toAlignmentIndex - fromAlignmentIndex + 1,
  1005.         origWidth = prevUABlock->m_width;
  1006.     vector < unsigned int > seqIndexesFrom(NRows()), seqIndexesTo(NRows());
  1007.     const Sequence *seq;
  1008. bool ignored;
  1009.     for (row=0; row<NRows(); ++row) {
  1010.         if (!GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, &seq, &seqIndexFrom, &ignored) ||
  1011.                 !GetSequenceAndIndexAt(toAlignmentIndex, row, justification, &seq, &seqIndexTo, &ignored) ||
  1012.                 seqIndexFrom < 0 || seqIndexTo < 0 ||
  1013.                 seqIndexTo - seqIndexFrom + 1 != newBlockWidth)
  1014.             return false;
  1015.         seqIndexesFrom[row] = seqIndexFrom;
  1016.         seqIndexesTo[row] = seqIndexTo;
  1017.     }
  1019.     TRACE_MESSAGE("creating new aligned and unaligned blocks");
  1021.     UnalignedBlock *nextUABlock = new UnalignedBlock(this);
  1022.     UngappedAlignedBlock *ABlock = new UngappedAlignedBlock(this);
  1023.     prevUABlock->m_width = nextUABlock->m_width = 0;
  1025.     bool deletePrevUABlock = true, deleteNextUABlock = true;
  1026.     const Block::Range *prevRange;
  1027.     int rangeWidth;
  1028.     for (row=0; row<NRows(); ++row) {
  1029.         prevRange = prevUABlock->GetRangeOfRow(row);
  1031.         nextUABlock->SetRangeOfRow(row, seqIndexesTo[row] + 1, prevRange->to);
  1032.         rangeWidth = prevRange->to - seqIndexesTo[row];
  1033.         if (rangeWidth < 0)
  1034.             ERROR_MESSAGE("BlockMultipleAlignment::CreateBlock() - negative nextRange width");
  1035.         else if (rangeWidth > 0) {
  1036.             if ((unsigned int) rangeWidth > nextUABlock->m_width)
  1037.                 nextUABlock->m_width = rangeWidth;
  1038.             deleteNextUABlock = false;
  1039.         }
  1041.         prevUABlock->SetRangeOfRow(row, prevRange->from, seqIndexesFrom[row] - 1);
  1042.         rangeWidth = seqIndexesFrom[row] - prevRange->from;
  1043.         if (rangeWidth < 0)
  1044.             ERROR_MESSAGE("BlockMultipleAlignment::CreateBlock() - negative prevRange width");
  1045.         else if (rangeWidth > 0) {
  1046.             if ((unsigned int) rangeWidth > prevUABlock->m_width)
  1047.                 prevUABlock->m_width = rangeWidth;
  1048.             deletePrevUABlock = false;
  1049.         }
  1051.         ABlock->SetRangeOfRow(row, seqIndexesFrom[row], seqIndexesTo[row]);
  1052.     }
  1054.     ABlock->m_width = newBlockWidth;
  1055.     if (prevUABlock->m_width + ABlock->m_width + nextUABlock->m_width != origWidth)
  1056.         ERROR_MESSAGE("BlockMultipleAlignment::CreateBlock() - bad block m_widths sum");
  1058.     InsertBlockAfter(prevUABlock, ABlock);
  1059.     InsertBlockAfter(ABlock, nextUABlock);
  1060.     if (deletePrevUABlock) {
  1061.         TRACE_MESSAGE("deleting zero-width unaligned block on left");
  1062.         RemoveBlock(prevUABlock);
  1063.     }
  1064.     if (deleteNextUABlock) {
  1065.         TRACE_MESSAGE("deleting zero-width unaligned block on right");
  1066.         RemoveBlock(nextUABlock);
  1067.     }
  1069.     if (!CheckAlignedBlock(ABlock))
  1070.         ERROR_MESSAGE("BlockMultipleAlignment::CreateBlock() - failed to create valid block");
  1072.     UpdateBlockMap();
  1073.     return true;
  1074. }
  1076. bool BlockMultipleAlignment::DeleteBlock(unsigned int alignmentIndex)
  1077. {
  1078.     Block *block = m_blockMap[alignmentIndex].block;
  1079.     if (!block || !block->IsAligned())
  1080.         return false;
  1082.     TRACE_MESSAGE("deleting block");
  1083.     Block
  1084.         *prevBlock = GetBlockBefore(block),
  1085.         *nextBlock = GetBlockAfter(block);
  1087.     // unaligned m_blocks on both sides - note that total alignment m_width can change!
  1088.     if (prevBlock && !prevBlock->IsAligned() && nextBlock && !nextBlock->IsAligned()) {
  1089.         const Block::Range *prevRange, *nextRange;
  1090.         unsigned int maxWidth = 0, width;
  1091.         for (unsigned int row=0; row<NRows(); ++row) {
  1092.             prevRange = prevBlock->GetRangeOfRow(row);
  1093.             nextRange = nextBlock->GetRangeOfRow(row);
  1094.             width = nextRange->to - prevRange->from + 1;
  1095.             prevBlock->SetRangeOfRow(row, prevRange->from, nextRange->to);
  1096.             if (width > maxWidth)
  1097.                 maxWidth = width;
  1098.         }
  1099.         prevBlock->m_width = maxWidth;
  1100.         TRACE_MESSAGE("removing extra unaligned block");
  1101.         RemoveBlock(nextBlock);
  1102.     }
  1104.     // unaligned block on left only
  1105.     else if (prevBlock && !prevBlock->IsAligned()) {
  1106.         const Block::Range *prevRange, *range;
  1107.         for (unsigned int row=0; row<NRows(); ++row) {
  1108.             prevRange = prevBlock->GetRangeOfRow(row);
  1109.             range = block->GetRangeOfRow(row);
  1110.             prevBlock->SetRangeOfRow(row, prevRange->from, range->to);
  1111.         }
  1112.         prevBlock->m_width += block->m_width;
  1113.     }
  1115.     // unaligned block on right only
  1116.     else if (nextBlock && !nextBlock->IsAligned()) {
  1117.         const Block::Range *range, *nextRange;
  1118.         for (unsigned int row=0; row<NRows(); ++row) {
  1119.             range = block->GetRangeOfRow(row);
  1120.             nextRange = nextBlock->GetRangeOfRow(row);
  1121.             nextBlock->SetRangeOfRow(row, range->from, nextRange->to);
  1122.         }
  1123.         nextBlock->m_width += block->m_width;
  1124.     }
  1126.     // no adjacent unaligned m_blocks
  1127.     else {
  1128.         TRACE_MESSAGE("creating new unaligned block");
  1129.         Block *newBlock = CreateNewUnalignedBlockBetween(prevBlock, nextBlock);
  1130.         if (prevBlock)
  1131.             InsertBlockAfter(prevBlock, newBlock);
  1132.         else if (nextBlock)
  1133.             InsertBlockBefore(newBlock, nextBlock);
  1134.         else
  1135.             m_blocks.push_back(CRef<Block>(newBlock));
  1136.     }
  1138.     RemoveBlock(block);
  1139.     UpdateBlockMap();
  1140.     return true;
  1141. }
  1143. bool BlockMultipleAlignment::DeleteAllBlocks(void)
  1144. {
  1145.     if (m_blocks.size() == 0)
  1146.         return false;
  1148.     m_blocks.clear();
  1149.     InitCache();
  1150.     AddUnalignedBlocks();   // one single unaligned block for whole alignment
  1151.     UpdateBlockMap();
  1152.     return true;
  1153. }
  1155. bool BlockMultipleAlignment::DeleteRow(unsigned int row)
  1156. {
  1157.     if (row >= NRows()) {
  1158.         ERROR_MESSAGE("BlockMultipleAlignment::DeleteRow() - row out of range");
  1159.         return false;
  1160.     }
  1162.     // remove sequence from list
  1163.     SequenceList::iterator s = m_sequences.begin();
  1164.     for (unsigned int i=0; i<row; ++i)
  1165.         ++s;
  1166.     m_sequences.erase(s);
  1168.     // delete row from all m_blocks, removing any zero-m_width m_blocks
  1169.     BlockList::iterator b = m_blocks.begin(), br, be = m_blocks.end();
  1170.     while (b != be) {
  1171.         (*b)->DeleteRow(row);
  1172.         if ((*b)->m_width == 0) {
  1173.             br = b;
  1174.             ++b;
  1175.             TRACE_MESSAGE("deleting block resized to zero width");
  1176.             RemoveBlock(*br);
  1177.         } else
  1178.             ++b;
  1179.     }
  1181.     // update total alignment m_width
  1182.     UpdateBlockMap();
  1183.     InitCache();
  1185.     return true;
  1186. }
  1188. void BlockMultipleAlignment::GetUngappedAlignedBlocks(UngappedAlignedBlockList *uabs) const
  1189. {
  1190.     uabs->clear();
  1191.     uabs->reserve(m_blocks.size());
  1192.     BlockList::const_iterator b, be = m_blocks.end();
  1193.     for (b=m_blocks.begin(); b!=be; ++b) {
  1194.         const UngappedAlignedBlock *uab = dynamic_cast<const UngappedAlignedBlock*>(b->GetPointer());
  1195.         if (uab)
  1196.             uabs->push_back(uab);
  1197.     }
  1198.     uabs->resize(uabs->size());
  1199. }
  1201. void BlockMultipleAlignment::GetModifiableUngappedAlignedBlocks(ModifiableUngappedAlignedBlockList *uabs)
  1202. {
  1203.     uabs->clear();
  1204.     uabs->reserve(m_blocks.size());
  1205.     BlockList::iterator b, be = m_blocks.end();
  1206.     for (b=m_blocks.begin(); b!=be; ++b) {
  1207.         UngappedAlignedBlock *uab = dynamic_cast<UngappedAlignedBlock*>(b->GetPointer());
  1208.         if (uab)
  1209.             uabs->push_back(uab);
  1210.     }
  1211.     uabs->resize(uabs->size());
  1212. }
  1214. bool BlockMultipleAlignment::ExtractRows(
  1215.     const vector < unsigned int >& slavesToRemove, AlignmentList *pairwiseAlignments)
  1216. {
  1217.     if (slavesToRemove.size() == 0) return false;
  1219.     // make a bool list of rows to remove, also checking to make sure slave list items are in range
  1220.     unsigned int i;
  1221.     vector < bool > removeRows(NRows(), false);
  1222.     for (i=0; i<slavesToRemove.size(); ++i) {
  1223.         if (slavesToRemove[i] > 0 && slavesToRemove[i] < NRows()) {
  1224.             removeRows[slavesToRemove[i]] = true;
  1225.         } else {
  1226.             ERROR_MESSAGE("BlockMultipleAlignment::ExtractRows() - can't extract row "
  1227.                 << slavesToRemove[i]);
  1228.             return false;
  1229.         }
  1230.     }
  1232.     if (pairwiseAlignments) {
  1233.         TRACE_MESSAGE("creating new pairwise alignments");
  1234.         SetDiagPostLevel(eDiag_Warning);    // otherwise, info messages take a long time if lots of rows
  1236.         UngappedAlignedBlockList uaBlocks;
  1237.         GetUngappedAlignedBlocks(&uaBlocks);
  1238.         UngappedAlignedBlockList::const_iterator u, ue = uaBlocks.end();
  1240.         for (i=0; i<slavesToRemove.size(); ++i) {
  1242.             // create new pairwise alignment from each removed row
  1243.             SequenceList newSeqs(2);
  1244.             newSeqs[0] = m_sequences[0];
  1245.             newSeqs[1] = m_sequences[slavesToRemove[i]];
  1246.             BlockMultipleAlignment *newAlignment = new BlockMultipleAlignment(newSeqs);
  1247.             for (u=uaBlocks.begin(); u!=ue; ++u) {
  1248.                 UngappedAlignedBlock *newABlock = new UngappedAlignedBlock(newAlignment);
  1249.                 const Block::Range *range = (*u)->GetRangeOfRow(0);
  1250.                 newABlock->SetRangeOfRow(0, range->from, range->to);
  1251.                 range = (*u)->GetRangeOfRow(slavesToRemove[i]);
  1252.                 newABlock->SetRangeOfRow(1, range->from, range->to);
  1253.                 newABlock->m_width = range->to - range->from + 1;
  1254.                 newAlignment->AddAlignedBlockAtEnd(newABlock);
  1255.             }
  1256.             if (!newAlignment->AddUnalignedBlocks() ||
  1257.                 !newAlignment->UpdateBlockMap()) {
  1258.                 ERROR_MESSAGE("BlockMultipleAlignment::ExtractRows() - error creating new alignment");
  1259.                 return false;
  1260.             }
  1262.             pairwiseAlignments->push_back(newAlignment);
  1263.         }
  1264.         SetDiagPostLevel(eDiag_Info);
  1265.     }
  1267.     // remove sequences
  1268.     TRACE_MESSAGE("deleting sequences");
  1269.     VectorRemoveElements(m_sequences, removeRows, slavesToRemove.size());
  1270.     VectorRemoveElements(m_rowDoubles, removeRows, slavesToRemove.size());
  1271.     VectorRemoveElements(m_rowStrings, removeRows, slavesToRemove.size());
  1273.     // delete row from all blocks, removing any zero-width blocks
  1274.     TRACE_MESSAGE("deleting alignment rows from blocks");
  1275.     BlockList::iterator b = m_blocks.begin(), br, be = m_blocks.end();
  1276.     while (b != be) {
  1277.         (*b)->DeleteRows(removeRows, slavesToRemove.size());
  1278.         if ((*b)->m_width == 0) {
  1279.             br = b;
  1280.             ++b;
  1281.             TRACE_MESSAGE("deleting block resized to zero width");
  1282.             RemoveBlock(*br);
  1283.         } else
  1284.             ++b;
  1285.     }
  1287.     // update total alignment m_width
  1288.     UpdateBlockMap();
  1289.     InitCache();
  1290.     return true;
  1291. }
  1293. bool BlockMultipleAlignment::MergeAlignment(const BlockMultipleAlignment *newAlignment)
  1294. {
  1295.     // check to see if the alignment is compatible - must have same master sequence
  1296.     // and blocks of new alignment must contain m_blocks of this alignment; at same time,
  1297.     // build up map of aligned blocks in new alignment that correspond to aligned blocks
  1298.     // of this object, for convenient lookup later
  1299.     if (newAlignment->GetMaster() != GetMaster())
  1300.         return false;
  1302.     const Block::Range *newRange, *thisRange;
  1303.     BlockList::const_iterator nb, nbe = newAlignment->m_blocks.end();
  1304.     BlockList::iterator b, be = m_blocks.end();
  1305.     typedef map < UngappedAlignedBlock *, const UngappedAlignedBlock * > AlignedBlockMap;
  1306.     AlignedBlockMap correspondingNewBlocks;
  1308.     for (b=m_blocks.begin(); b!=be; ++b) {
  1309.         if (!(*b)->IsAligned())
  1310.             continue;
  1311.         for (nb=newAlignment->m_blocks.begin(); nb!=nbe; ++nb) {
  1312.             if (!(*nb)->IsAligned())
  1313.                 continue;
  1315.             newRange = (*nb)->GetRangeOfRow(0);
  1316.             thisRange = (*b)->GetRangeOfRow(0);
  1317.             if (newRange->from <= thisRange->from && newRange->to >= thisRange->to) {
  1318.                 correspondingNewBlocks[dynamic_cast<UngappedAlignedBlock*>(b->GetPointer())] =
  1319.                     dynamic_cast<const UngappedAlignedBlock*>(nb->GetPointer());
  1320.                 break;
  1321.             }
  1322.         }
  1323.         if (nb == nbe) return false;    // no corresponding block found
  1324.     }
  1326.     // add slave sequences from new alignment; also copy scores/status
  1327.     unsigned int i, nNewRows = newAlignment->m_sequences.size() - 1;
  1328.     m_sequences.resize(m_sequences.size() + nNewRows);
  1329.     m_rowDoubles.resize(m_rowDoubles.size() + nNewRows);
  1330.     m_rowStrings.resize(m_rowStrings.size() + nNewRows);
  1331.     for (i=0; i<nNewRows; ++i) {
  1332.         m_sequences[m_sequences.size() + i - nNewRows] = newAlignment->m_sequences[i + 1];
  1333.         SetRowDouble(NRows() + i - nNewRows, newAlignment->GetRowDouble(i + 1));
  1334.         SetRowStatusLine(NRows() + i - nNewRows, newAlignment->GetRowStatusLine(i + 1));
  1335.     }
  1337.     // now that we know blocks are compatible, add new rows at end of this alignment, containing
  1338.     // all rows (except master) from new alignment; only that part of the aligned blocks from
  1339.     // the new alignment that intersect with the aligned blocks from this object are added, so
  1340.     // that this object's block structure is unchanged
  1342.     // resize all aligned blocks
  1343.     for (b=m_blocks.begin(); b!=be; ++b)
  1344.         (*b)->AddRows(nNewRows);
  1346.     // set ranges of aligned blocks, and add them to the list
  1347.     AlignedBlockMap::const_iterator ab, abe = correspondingNewBlocks.end();
  1348.     const Block::Range *thisMaster, *newMaster;
  1349.     for (ab=correspondingNewBlocks.begin(); ab!=abe; ++ab) {
  1350.         thisMaster = ab->first->GetRangeOfRow(0);
  1351.         newMaster = ab->second->GetRangeOfRow(0);
  1352.         for (i=0; i<nNewRows; ++i) {
  1353.             newRange = ab->second->GetRangeOfRow(i + 1);
  1354.             ab->first->SetRangeOfRow(NRows() + i - nNewRows,
  1355.                 newRange->from + thisMaster->from - newMaster->from,
  1356.                 newRange->to + thisMaster->to - newMaster->to);
  1357.         }
  1358.     }
  1360.     // delete then recreate the unaligned blocks, in case a merge requires the
  1361.     // creation of a new unaligned block
  1362.     for (b=m_blocks.begin(); b!=be; ) {
  1363.         if (!(*b)->IsAligned()) {
  1364.             BlockList::iterator bb(b);
  1365.             ++bb;
  1366.             m_blocks.erase(b);
  1367.             b = bb;
  1368.         } else
  1369.             ++b;
  1370.     }
  1371.     InitCache();
  1373.     // update this alignment, but leave row scores/status alone
  1374.     if (!AddUnalignedBlocks() || !UpdateBlockMap(false)) {
  1375.         ERROR_MESSAGE("BlockMultipleAlignment::MergeAlignment() - internal update after merge failed");
  1376.         return false;
  1377.     }
  1378.     return true;
  1379. }
  1381. template < class T >
  1382. bool ReorderVector(T& v, const std::vector < unsigned int >& newOrder)
  1383. {
  1384.     // check validity of new ordering
  1385.     if (newOrder.size() != v.size()) {
  1386.         ERROR_MESSAGE("ReorderVector() - wrong size newOrder");
  1387.         return false;
  1388.     }
  1389.     vector < bool > isPresent(v.size(), false);
  1390.     unsigned int r;
  1391.     for (r=0; r<v.size(); r++) {
  1392.         if (isPresent[newOrder[r]]) {
  1393.             ERROR_MESSAGE("ReorderVector() - invalid newOrder: repeated/missing row");
  1394.             return false;
  1395.         }
  1396.         isPresent[newOrder[r]] = true;
  1397.     }
  1399.     // not terribly efficient - makes a whole new copy with the new order, then re-copies back
  1400.     T reordered(v.size());
  1401.     for (r=0; r<v.size(); r++)
  1402.         reordered[r] = v[newOrder[r]];
  1403.     v = reordered;
  1405.     return true;
  1406. }
  1408. bool BlockMultipleAlignment::ReorderRows(const std::vector < unsigned int >& newOrder)
  1409. {
  1410.     // can't reorder master
  1411.     if (newOrder[0] != 0) {
  1412.         ERROR_MESSAGE("ReorderRows() - can't move master row");
  1413.         return false;
  1414.     }
  1415.     bool okay =
  1416.         (ReorderVector(m_sequences, newOrder) &&
  1417.          ReorderVector(m_rowDoubles, newOrder) &&
  1418.          ReorderVector(m_rowStrings, newOrder));
  1419.     if (!okay) {
  1420.         ERROR_MESSAGE("reordering of sequences and status info failed");
  1421.         return false;
  1422.     }
  1423.     BlockList::iterator b, be = m_blocks.end();
  1424.     for (b=m_blocks.begin(); b!=be; ++b)
  1425.         okay = (okay && (*b)->ReorderRows(newOrder));
  1426.     if (!okay)
  1427.         ERROR_MESSAGE("reordering of block ranges failed");
  1428.     return okay;
  1429. }
  1431. CSeq_align * CreatePairwiseSeqAlignFromMultipleRow(const BlockMultipleAlignment *multiple,
  1432.     const BlockMultipleAlignment::UngappedAlignedBlockList& m_blocks, unsigned int slaveRow)
  1433. {
  1434.     if (!multiple || slaveRow >= multiple->NRows()) {
  1435.         ERROR_MESSAGE("CreatePairwiseSeqAlignFromMultipleRow() - bad parameters");
  1436.         return NULL;
  1437.     }
  1439.     CSeq_align *seqAlign = new CSeq_align();
  1440.     seqAlign->SetType(CSeq_align::eType_partial);
  1441.     seqAlign->SetDim(2);
  1443.     CSeq_align::C_Segs::TDendiag& denDiags = seqAlign->SetSegs().SetDendiag();
  1444.     denDiags.resize((m_blocks.size() > 0) ? m_blocks.size() : 1);
  1446.     CSeq_align::C_Segs::TDendiag::iterator d, de = denDiags.end();
  1447.     BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b = m_blocks.begin();
  1448.     const Block::Range *range;
  1449.     for (d=denDiags.begin(); d!=de; ++d, ++b) {
  1451.         CDense_diag *denDiag = new CDense_diag();
  1452.         d->Reset(denDiag);
  1453.         denDiag->SetDim(2);
  1454.         denDiag->SetIds().resize(2);
  1456.         // master row
  1457.         CRef < CSeq_id > id(new CSeq_id());
  1458.         id->Assign(multiple->GetSequenceOfRow(0)->GetPreferredIdentifier());
  1459.         denDiag->SetIds().front() = id;
  1460.         if (m_blocks.size() > 0) {
  1461.             range = (*b)->GetRangeOfRow(0);
  1462.             denDiag->SetStarts().push_back(range->from);
  1463.         } else
  1464.             denDiag->SetStarts().push_back(0);
  1466.         // slave row
  1467.         id.Reset(new CSeq_id());
  1468.         id->Assign(multiple->GetSequenceOfRow(slaveRow)->GetPreferredIdentifier());
  1469.         denDiag->SetIds().back() = id;
  1470.         if (m_blocks.size() > 0) {
  1471.             range = (*b)->GetRangeOfRow(slaveRow);
  1472.             denDiag->SetStarts().push_back(range->from);
  1473.         } else
  1474.             denDiag->SetStarts().push_back(0);
  1476.         // block m_width
  1477.         denDiag->SetLen((m_blocks.size() > 0) ? (*b)->m_width : 0);
  1478.     }
  1480.     return seqAlign;
  1481. }
  1483. unsigned int BlockMultipleAlignment::NAlignedBlocks(void) const
  1484. {
  1485.     unsigned int n = 0;
  1486.     BlockList::const_iterator b, be = m_blocks.end();
  1487.     for (b=m_blocks.begin(); b!=be; ++b)
  1488.         if ((*b)->IsAligned())
  1489.             ++n;
  1490.     return n;
  1491. }
  1493. unsigned int BlockMultipleAlignment::GetAlignmentIndex(unsigned int row, unsigned int seqIndex, eUnalignedJustification justification)
  1494. {
  1495.     if (row >= NRows() || seqIndex >= GetSequenceOfRow(row)->Length()) {
  1496.         ERROR_MESSAGE("BlockMultipleAlignment::GetAlignmentIndex() - coordinate out of range");
  1497.         return eUndefined;
  1498.     }
  1500.     unsigned int alignmentIndex, blockColumn;
  1501.     const Block *block = NULL;
  1502.     const Block::Range *range;
  1504.     for (alignmentIndex=0; alignmentIndex<m_totalWidth; ++alignmentIndex) {
  1506.         // check each block to see if index is in range
  1507.         if (block != m_blockMap[alignmentIndex].block) {
  1508.             block = m_blockMap[alignmentIndex].block;
  1510.             range = block->GetRangeOfRow(row);
  1511.             if (seqIndex >= range->from && seqIndex <= range->to) {
  1513.                 // override requested justification for end blocks
  1514.                 if (block == m_blocks.back()) // also true if there's a single aligned block
  1515.                     justification = eLeft;
  1516.                 else if (block == m_blocks.front())
  1517.                     justification = eRight;
  1519.                 // search block columns to find index (inefficient, but avoids having to write
  1520.                 // inverse functions of Block::GetIndexAt()
  1521.                 for (blockColumn=0; blockColumn<block->m_width; ++blockColumn) {
  1522.                     if (seqIndex == block->GetIndexAt(blockColumn, row, justification))
  1523.                         return alignmentIndex + blockColumn;
  1524.                 }
  1525.                 ERROR_MESSAGE("BlockMultipleAlignment::GetAlignmentIndex() - can't find index in block");
  1526.                 return eUndefined;
  1527.             }
  1528.         }
  1529.     }
  1531.     // should never get here...
  1532.     ERROR_MESSAGE("BlockMultipleAlignment::GetAlignmentIndex() - confused");
  1533.     return eUndefined;
  1534. }
  1537. bool Block::ReorderRows(const std::vector < unsigned int >& newOrder)
  1538. {
  1539.     return ReorderVector(m_ranges, newOrder);
  1540. }
  1543. ///// UngappedAlignedBlock methods /////
  1545. char UngappedAlignedBlock::GetCharacterAt(unsigned int blockColumn, unsigned int row) const
  1546. {
  1547.     return m_parentAlignment->GetSequenceOfRow(row)->m_sequenceString[GetIndexAt(blockColumn, row)];
  1548. }
  1550. Block * UngappedAlignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
  1551. {
  1552.     UngappedAlignedBlock *copy = new UngappedAlignedBlock(newMultiple);
  1553.     const Block::Range *range;
  1554.     for (unsigned int row=0; row<NSequences(); ++row) {
  1555.         range = GetRangeOfRow(row);
  1556.         copy->SetRangeOfRow(row, range->from, range->to);
  1557.     }
  1558.     copy->m_width = m_width;
  1559.     return copy;
  1560. }
  1562. void UngappedAlignedBlock::DeleteRow(unsigned int row)
  1563. {
  1564.     RangeList::iterator r = m_ranges.begin();
  1565.     for (unsigned int i=0; i<row; ++i)
  1566.         ++r;
  1567.     m_ranges.erase(r);
  1568. }
  1570. void UngappedAlignedBlock::DeleteRows(vector < bool >& removeRows, unsigned int nToRemove)
  1571. {
  1572.     VectorRemoveElements(m_ranges, removeRows, nToRemove);
  1573. }
  1576. ///// UnalignedBlock methods /////
  1578. unsigned int UnalignedBlock::GetIndexAt(unsigned int blockColumn, unsigned int row,
  1579.         BlockMultipleAlignment::eUnalignedJustification justification) const
  1580. {
  1581.     const Block::Range *range = GetRangeOfRow(row);
  1582.     unsigned int seqIndex = BlockMultipleAlignment::eUndefined, rangeWidth, rangeMiddle, extraSpace;
  1584.     switch (justification) {
  1585.         case BlockMultipleAlignment::eLeft:
  1586.             seqIndex = range->from + blockColumn;
  1587.             break;
  1588.         case BlockMultipleAlignment::eRight:
  1589.             seqIndex = range->to - m_width + blockColumn + 1;
  1590.             break;
  1591.         case BlockMultipleAlignment::eCenter:
  1592.             rangeWidth = (range->to - range->from + 1);
  1593.             extraSpace = (m_width - rangeWidth) / 2;
  1594.             if (blockColumn < extraSpace || blockColumn >= extraSpace + rangeWidth)
  1595.                 seqIndex = BlockMultipleAlignment::eUndefined;
  1596.             else
  1597.                 seqIndex = range->from + blockColumn - extraSpace;
  1598.             break;
  1599.         case BlockMultipleAlignment::eSplit:
  1600.             rangeWidth = (range->to - range->from + 1);
  1601.             rangeMiddle = (rangeWidth / 2) + (rangeWidth % 2);
  1602.             extraSpace = m_width - rangeWidth;
  1603.             if (blockColumn < rangeMiddle)
  1604.                 seqIndex = range->from + blockColumn;
  1605.             else if (blockColumn >= extraSpace + rangeMiddle)
  1606.                 seqIndex = range->to - m_width + blockColumn + 1;
  1607.             else
  1608.                 seqIndex = BlockMultipleAlignment::eUndefined;
  1609.             break;
  1610.     }
  1611.     if (seqIndex < range->from || seqIndex > range->to)
  1612.         seqIndex = BlockMultipleAlignment::eUndefined;
  1614.     return seqIndex;
  1615. }
  1617. void UnalignedBlock::Resize(void)
  1618. {
  1619.     m_width = 0;
  1620.     for (unsigned int i=0; i<NSequences(); ++i) {
  1621.         unsigned int blockWidth = m_ranges[i].to - m_ranges[i].from + 1;
  1622.         if (blockWidth > m_width)
  1623.             m_width = blockWidth;
  1624.     }
  1625. }
  1627. void UnalignedBlock::DeleteRow(unsigned int row)
  1628. {
  1629.     RangeList::iterator r = m_ranges.begin();
  1630.     for (unsigned int i=0; i<row; ++i)
  1631.         ++r;
  1632.     m_ranges.erase(r);
  1633.     Resize();
  1634. }
  1636. void UnalignedBlock::DeleteRows(vector < bool >& removeRows, unsigned int nToRemove)
  1637. {
  1638.     VectorRemoveElements(m_ranges, removeRows, nToRemove);
  1639.     Resize();
  1640. }
  1642. Block * UnalignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
  1643. {
  1644.     UnalignedBlock *copy = new UnalignedBlock(newMultiple);
  1645.     const Block::Range *range;
  1646.     for (unsigned int row=0; row<NSequences(); ++row) {
  1647.         range = GetRangeOfRow(row);
  1648.         copy->SetRangeOfRow(row, range->from, range->to);
  1649.     }
  1650.     copy->m_width = m_width;
  1651.     return copy;
  1652. }
  1654. END_SCOPE(struct_util)
  1656. /*
  1657. * ---------------------------------------------------------------------------
  1658. * $Log: su_block_multiple_alignment.cpp,v $
  1659. * Revision 1000.0  2004/06/01 18:14:28  gouriano
  1660. * PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.9
  1661. *
  1662. * Revision 1.9  2004/05/28 10:07:39  thiessen
  1663. * fix GCC warning
  1664. *
  1665. * Revision 1.8  2004/05/28 09:46:57  thiessen
  1666. * restructure C-toolkit header usage ; move C Bioseq storage into su_sequence_set
  1667. *
  1668. * Revision 1.7  2004/05/27 21:34:08  thiessen
  1669. * add PSSM calculation (requires C-toolkit)
  1670. *
  1671. * Revision 1.6  2004/05/26 14:49:59  thiessen
  1672. * UNDEFINED -> eUndefined
  1673. *
  1674. * Revision 1.5  2004/05/26 14:30:16  thiessen
  1675. * adjust handling of alingment data ; add row ordering
  1676. *
  1677. * Revision 1.4  2004/05/26 02:40:24  thiessen
  1678. * progress towards LOO - all but PSSM and row ordering
  1679. *
  1680. * Revision 1.3  2004/05/25 16:24:50  thiessen
  1681. * remove WorkShop warnings
  1682. *
  1683. * Revision 1.2  2004/05/25 16:12:30  thiessen
  1684. * fix GCC warnings
  1685. *
  1686. * Revision 1.1  2004/05/25 15:52:17  thiessen
  1687. * add BlockMultipleAlignment, IBM algorithm
  1688. *
  1689. */