1 /* $Id: block_multiple_alignment.cpp 103491 2007-05-04 17:18:18Z kazimird $
2 * ===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Authors: Paul Thiessen
27 *
28 * File Description:
29 * Classes to hold alignment data
30 *
31 * ===========================================================================
32 */
33
34 #include <ncbi_pch.hpp>
35 #include <corelib/ncbistd.hpp>
36
37 #include <objects/seqalign/Dense_diag.hpp>
38
39 #include "remove_header_conflicts.hpp"
40
41 #include "block_multiple_alignment.hpp"
42 #include "sequence_set.hpp"
43 #include "molecule.hpp"
44 #include "conservation_colorer.hpp"
45 #include "style_manager.hpp"
46 #include "structure_set.hpp"
47 #include "messenger.hpp"
48 #include "cn3d_colors.hpp"
49 #include "alignment_manager.hpp"
50 #include "cn3d_tools.hpp"
51 #include "molecule_identifier.hpp"
52 #include "cn3d_threader.hpp"
53 #include "cn3d_pssm.hpp"
54
55 USING_NCBI_SCOPE;
56 USING_SCOPE(objects);
57
58
BEGIN_SCOPE(Cn3D)59 BEGIN_SCOPE(Cn3D)
60
61 BlockMultipleAlignment::BlockMultipleAlignment(SequenceList *sequenceList, AlignmentManager *alnMgr) :
62 alignmentManager(alnMgr), sequences(sequenceList), conservationColorer(NULL),
63 pssm(NULL), showGeometryViolations(false),
64 alignMasterFrom(-1), alignMasterTo(-1), alignDependentFrom(-1), alignDependentTo(-1)
65 {
66 InitCache();
67 rowDoubles.resize(sequenceList->size(), 0.0);
68 rowStrings.resize(sequenceList->size());
69 geometryViolations.resize(sequenceList->size());
70
71 // create conservation colorer
72 conservationColorer = new ConservationColorer(this);
73 }
74
InitCache(void)75 void BlockMultipleAlignment::InitCache(void)
76 {
77 cachePrevRow = kMax_UInt;
78 cachePrevBlock = NULL;
79 cacheBlockIterator = blocks.begin();
80 }
81
~BlockMultipleAlignment(void)82 BlockMultipleAlignment::~BlockMultipleAlignment(void)
83 {
84 BlockList::iterator i, ie = blocks.end();
85 for (i=blocks.begin(); i!=ie; ++i) delete *i;
86 delete sequences;
87 delete conservationColorer;
88 RemovePSSM();
89 }
90
Clone(void) const91 BlockMultipleAlignment * BlockMultipleAlignment::Clone(void) const
92 {
93 // must actually copy the list
94 BlockMultipleAlignment *copy = new BlockMultipleAlignment(new SequenceList(*sequences), alignmentManager);
95 BlockList::const_iterator b, be = blocks.end();
96 for (b=blocks.begin(); b!=be; ++b) {
97 Block *newBlock = (*b)->Clone(copy);
98 copy->blocks.push_back(newBlock);
99 if ((*b)->IsAligned()) {
100 MarkBlockMap::const_iterator r = markBlocks.find(*b);
101 if (r != markBlocks.end())
102 copy->markBlocks[newBlock] = true;
103 }
104 }
105 copy->UpdateBlockMapAndColors();
106 copy->rowDoubles = rowDoubles;
107 copy->rowStrings = rowStrings;
108 copy->geometryViolations = geometryViolations;
109 copy->showGeometryViolations = showGeometryViolations;
110 copy->updateOrigin = updateOrigin;
111 copy->alignMasterFrom = alignMasterFrom;
112 copy->alignMasterTo = alignMasterTo;
113 copy->alignDependentFrom = alignDependentFrom;
114 copy->alignDependentTo = alignDependentTo;
115 return copy;
116 }
117
GetPSSM(void) const118 const PSSMWrapper& BlockMultipleAlignment::GetPSSM(void) const
119 {
120 if (!pssm)
121 pssm = new PSSMWrapper(this);
122 return *pssm;
123 }
124
RemovePSSM(void) const125 void BlockMultipleAlignment::RemovePSSM(void) const
126 {
127 if (pssm) {
128 delete pssm;
129 pssm = NULL;
130 }
131 }
132
FreeColors(void)133 void BlockMultipleAlignment::FreeColors(void)
134 {
135 conservationColorer->FreeColors();
136 RemovePSSM();
137 }
138
CheckAlignedBlock(const Block * block) const139 bool BlockMultipleAlignment::CheckAlignedBlock(const Block *block) const
140 {
141 if (!block || !block->IsAligned()) {
142 ERRORMSG("BlockMultipleAlignment::CheckAlignedBlock() - checks aligned blocks only");
143 return false;
144 }
145 if (block->NSequences() != sequences->size()) {
146 ERRORMSG("BlockMultipleAlignment::CheckAlignedBlock() - block size mismatch");
147 return false;
148 }
149
150 // make sure ranges are reasonable for each sequence
151 unsigned int row;
152 const Block
153 *prevBlock = GetBlockBefore(block),
154 *nextBlock = GetBlockAfter(block);
155 const Block::Range *range, *prevRange = NULL, *nextRange = NULL;
156 SequenceList::const_iterator sequence = sequences->begin();
157 for (row=0; row<block->NSequences(); ++row, ++sequence) {
158 range = block->GetRangeOfRow(row);
159 if (prevBlock) prevRange = prevBlock->GetRangeOfRow(row);
160 if (nextBlock) nextRange = nextBlock->GetRangeOfRow(row);
161 if (range->to - range->from + 1 != (int)block->width || // check block width
162 (prevRange && range->from <= prevRange->to) || // check for range overlap
163 (nextRange && range->to >= nextRange->from) || // check for range overlap
164 range->from > range->to || // check range values
165 range->to >= (int)(*sequence)->Length()) { // check bounds of end
166 ERRORMSG("BlockMultipleAlignment::CheckAlignedBlock() - range error");
167 return false;
168 }
169 }
170
171 return true;
172 }
173
AddAlignedBlockAtEnd(UngappedAlignedBlock * newBlock)174 bool BlockMultipleAlignment::AddAlignedBlockAtEnd(UngappedAlignedBlock *newBlock)
175 {
176 bool okay = CheckAlignedBlock(newBlock);
177 if (okay)
178 blocks.push_back(newBlock);
179 return okay;
180 }
181
182 UnalignedBlock * BlockMultipleAlignment::
CreateNewUnalignedBlockBetween(const Block * leftBlock,const Block * rightBlock)183 CreateNewUnalignedBlockBetween(const Block *leftBlock, const Block *rightBlock)
184 {
185 if ((leftBlock && !leftBlock->IsAligned()) || (rightBlock && !rightBlock->IsAligned())) {
186 ERRORMSG("CreateNewUnalignedBlockBetween() - passed an unaligned block");
187 return NULL;
188 }
189
190 unsigned int row, length;
191 int from, to;
192 SequenceList::const_iterator s, se = sequences->end();
193
194 UnalignedBlock *newBlock = new UnalignedBlock(this);
195 newBlock->width = 0;
196
197 for (row=0, s=sequences->begin(); s!=se; ++row, ++s) {
198
199 if (leftBlock)
200 from = leftBlock->GetRangeOfRow(row)->to + 1;
201 else
202 from = 0;
203
204 if (rightBlock)
205 to = rightBlock->GetRangeOfRow(row)->from - 1;
206 else
207 to = (*s)->Length() - 1;
208
209 newBlock->SetRangeOfRow(row, from, to);
210
211 if (to + 1 < from) { // just to make sure...
212 ERRORMSG("CreateNewUnalignedBlockBetween() - unaligned length < 0");
213 return NULL;
214 }
215 length = to - from + 1;
216 if (length > newBlock->width)
217 newBlock->width = length;
218 }
219
220 if (newBlock->width == 0) {
221 delete newBlock;
222 return NULL;
223 } else
224 return newBlock;
225 }
226
AddUnalignedBlocks(void)227 bool BlockMultipleAlignment::AddUnalignedBlocks(void)
228 {
229 BlockList::iterator a, ae = blocks.end();
230 const Block *alignedBlock = NULL, *prevAlignedBlock = NULL;
231 Block *newUnalignedBlock;
232
233 // unaligned blocks to the left of each aligned block
234 for (a=blocks.begin(); a!=ae; ++a) {
235 alignedBlock = *a;
236 newUnalignedBlock = CreateNewUnalignedBlockBetween(prevAlignedBlock, alignedBlock);
237 if (newUnalignedBlock) blocks.insert(a, newUnalignedBlock);
238 prevAlignedBlock = alignedBlock;
239 }
240
241 // right tail
242 newUnalignedBlock = CreateNewUnalignedBlockBetween(alignedBlock, NULL);
243 if (newUnalignedBlock)
244 blocks.insert(a, newUnalignedBlock);
245
246 return true;
247 }
248
UpdateBlockMapAndColors(bool clearRowInfo)249 bool BlockMultipleAlignment::UpdateBlockMapAndColors(bool clearRowInfo)
250 {
251 unsigned int i = 0, j;
252 int n = 0;
253 BlockList::iterator b, be = blocks.end();
254
255 // reset old stuff, recalculate width
256 totalWidth = 0;
257 for (b=blocks.begin(); b!=be; ++b) totalWidth += (*b)->width;
258 // TESTMSG("alignment display size: " << totalWidth << " x " << NRows());
259
260 // fill out the block map
261 conservationColorer->Clear();
262 blockMap.resize(totalWidth);
263 UngappedAlignedBlock *aBlock;
264 for (b=blocks.begin(); b!=be; ++b) {
265 aBlock = dynamic_cast<UngappedAlignedBlock*>(*b);
266 if (aBlock) {
267 conservationColorer->AddBlock(aBlock);
268 ++n;
269 }
270 for (j=0; j<(*b)->width; ++j, ++i) {
271 blockMap[i].block = *b;
272 blockMap[i].blockColumn = j;
273 blockMap[i].alignedBlockNum = (aBlock ? n : -1);
274 }
275 }
276
277 // if alignment changes, any pssm/scores/status/special colors become invalid
278 RemovePSSM();
279 if (clearRowInfo) ClearRowInfo();
280 ShowGeometryViolations(showGeometryViolations); // recalculate GV's
281
282 return true;
283 }
284
GetCharacterTraitsAt(unsigned int alignmentColumn,unsigned int row,eUnalignedJustification justification,char * character,Vector * color,bool * isHighlighted,bool * drawBackground,Vector * cellBackgroundColor) const285 bool BlockMultipleAlignment::GetCharacterTraitsAt(
286 unsigned int alignmentColumn, unsigned int row, eUnalignedJustification justification,
287 char *character, Vector *color, bool *isHighlighted,
288 bool *drawBackground, Vector *cellBackgroundColor) const
289 {
290 const Sequence *sequence;
291 int seqIndex;
292 bool isAligned;
293
294 if (!GetSequenceAndIndexAt(alignmentColumn, row, justification, &sequence, &seqIndex, &isAligned))
295 return false;
296
297 *character = (seqIndex >= 0) ? sequence->sequenceString[seqIndex] : '~';
298 if (isAligned)
299 *character = toupper((unsigned char)(*character));
300 else
301 *character = tolower((unsigned char)(*character));
302
303 // try to color by molecule first
304 if (sequence->molecule) {
305 *color = (seqIndex >= 0) ?
306 sequence->molecule->GetResidueColor(seqIndex) :
307 GlobalColors()->Get(Colors::eNoCoordinates);;
308 }
309
310 // otherwise (for unstructured sequence):
311 else {
312 StyleSettings::eColorScheme colorScheme =
313 sequence->parentSet->styleManager->GetGlobalStyle().proteinBackbone.colorScheme;
314
315 // color by hydrophobicity
316 if (sequence->isProtein && colorScheme == StyleSettings::eHydrophobicity) {
317 double hydrophobicity = GetHydrophobicity(toupper((unsigned char)(*character)));
318 *color = (hydrophobicity != UNKNOWN_HYDROPHOBICITY) ?
319 GlobalColors()->Get(Colors::eHydrophobicityMap, hydrophobicity) :
320 GlobalColors()->Get(Colors::eNoHydrophobicity);
321 }
322 // or color by charge
323 else if (sequence->isProtein && colorScheme == StyleSettings::eCharge) {
324 int charge = GetCharge(toupper((unsigned char)(*character)));
325 *color = GlobalColors()->Get(
326 (charge > 0) ? Colors::ePositive : ((charge < 0) ? Colors::eNegative : Colors::eNeutral));;
327 }
328 // else, color by alignment color
329 else {
330 const Vector *aColor;
331 if (isAligned && (aColor = GetAlignmentColor(row, seqIndex, colorScheme)) != NULL) {
332 *color = *aColor;
333 } else {
334 *color = GlobalColors()->Get(Colors::eUnaligned);
335 }
336 }
337 }
338
339 if (seqIndex >= 0)
340 *isHighlighted = GlobalMessenger()->IsHighlighted(sequence, seqIndex);
341 else
342 *isHighlighted = false;
343
344 // add special alignment coloring (but don't override highlight)
345 *drawBackground = false;
346 if (!(*isHighlighted)) {
347
348 // check for block flagged for realignment
349 if (markBlocks.find(blockMap[alignmentColumn].block) != markBlocks.end()) {
350 *drawBackground = true;
351 *cellBackgroundColor = GlobalColors()->Get(Colors::eMarkBlock);
352 }
353
354 // optionally show geometry violations
355 if (showGeometryViolations && seqIndex >= 0 && geometryViolations[row][seqIndex]) {
356 *drawBackground = true;
357 *cellBackgroundColor = GlobalColors()->Get(Colors::eGeometryViolation);
358 }
359
360 // check for unmergeable alignment
361 const BlockMultipleAlignment *referenceAlignment = alignmentManager->GetCurrentMultipleAlignment();
362 if (referenceAlignment && referenceAlignment != this &&
363 seqIndex >= 0 && GetMaster() == referenceAlignment->GetMaster()) {
364
365 bool unmergeable = false;
366 const Block *block = blockMap[alignmentColumn].block;
367
368 // case where master residues are aligned in multiple but not in this one
369 if (row == 0 && !isAligned && referenceAlignment->IsAligned(row, seqIndex))
370 unmergeable = true;
371
372 // block boundaries in this inside aligned block of multiple
373 else if (row > 0 && isAligned) {
374 const Block::Range *range = block->GetRangeOfRow(row);
375 bool
376 isLeftEdge = (seqIndex == range->from),
377 isRightEdge = (seqIndex == range->to);
378 if (isLeftEdge || isRightEdge) {
379
380 // get corresponding block of multiple
381 const Block::Range *masterRange = block->GetRangeOfRow(0);
382 unsigned int masterIndex = masterRange->from + seqIndex - range->from;
383 const Block *multipleBlock = referenceAlignment->GetBlock(0, masterIndex);
384 masterRange = multipleBlock->GetRangeOfRow(0);
385
386 if (multipleBlock->IsAligned() &&
387 ((isLeftEdge && (int)masterIndex > masterRange->from) ||
388 (isRightEdge && (int)masterIndex < masterRange->to)))
389 unmergeable = true;
390 }
391 }
392
393 // check for inserts relative to the multiple
394 else if (row > 0 && !isAligned) {
395 const Block
396 *blockBefore = GetBlockBefore(block),
397 *blockAfter = GetBlockAfter(block),
398 *refBlock;
399 if (blockBefore && blockBefore->IsAligned() && blockAfter && blockAfter->IsAligned() &&
400 referenceAlignment->GetBlock(0, blockBefore->GetRangeOfRow(0)->to) ==
401 (refBlock=referenceAlignment->GetBlock(0, blockAfter->GetRangeOfRow(0)->from)) &&
402 refBlock->IsAligned()
403 )
404 unmergeable = true;
405 }
406
407 if (unmergeable) {
408 *drawBackground = true;
409 *cellBackgroundColor = GlobalColors()->Get(Colors::eMergeFail);
410 }
411 }
412 }
413
414 return true;
415 }
416
GetSequenceAndIndexAt(unsigned int alignmentColumn,unsigned int row,eUnalignedJustification requestedJustification,const Sequence ** sequence,int * index,bool * isAligned) const417 bool BlockMultipleAlignment::GetSequenceAndIndexAt(
418 unsigned int alignmentColumn, unsigned int row, eUnalignedJustification requestedJustification,
419 const Sequence **sequence, int *index, bool *isAligned) const
420 {
421 if (sequence) *sequence = (*sequences)[row];
422
423 const BlockInfo& blockInfo = blockMap[alignmentColumn];
424
425 if (!blockInfo.block->IsAligned()) {
426 if (isAligned) *isAligned = false;
427 // override requested justification for end blocks
428 if (blockInfo.block == blocks.back()) // also true if there's a single aligned block
429 requestedJustification = eLeft;
430 else if (blockInfo.block == blocks.front())
431 requestedJustification = eRight;
432 } else
433 if (isAligned) *isAligned = true;
434
435 if (index)
436 *index = blockInfo.block->GetIndexAt(blockInfo.blockColumn, row, requestedJustification);
437
438 return true;
439 }
440
GetRowForSequence(const Sequence * sequence) const441 int BlockMultipleAlignment::GetRowForSequence(const Sequence *sequence) const
442 {
443 // this only works for structured sequences, since non-structure sequences can
444 // be repeated any number of times in the alignment; assumes repeated structures
445 // will each have a unique Sequence object
446 if (!sequence || !sequence->molecule) {
447 ERRORMSG("BlockMultipleAlignment::GetRowForSequence() - Sequence must have associated structure");
448 return -1;
449 }
450
451 if (cachePrevRow >= NRows() || sequence != (*sequences)[cachePrevRow]) {
452 unsigned int row;
453 for (row=0; row<NRows(); ++row) if ((*sequences)[row] == sequence) break;
454 if (row == NRows()) {
455 // ERRORMSG("BlockMultipleAlignment::GetRowForSequence() - can't find given Sequence");
456 return -1;
457 }
458 cachePrevRow = row;
459 }
460 return cachePrevRow;
461 }
462
GetAlignmentColor(unsigned int row,unsigned int seqIndex,StyleSettings::eColorScheme colorScheme) const463 const Vector * BlockMultipleAlignment::GetAlignmentColor(
464 unsigned int row, unsigned int seqIndex, StyleSettings::eColorScheme colorScheme) const
465 {
466 const UngappedAlignedBlock *block = dynamic_cast<const UngappedAlignedBlock*>(GetBlock(row, seqIndex));
467 if (!block || !block->IsAligned()) {
468 WARNINGMSG("BlockMultipleAlignment::GetAlignmentColor() called on unaligned residue");
469 return NULL;
470 }
471
472 const Vector *alignedColor;
473
474 switch (colorScheme) {
475 case StyleSettings::eAligned:
476 alignedColor = GlobalColors()->Get(Colors::eConservationMap, Colors::nConservationMap - 1);
477 break;
478 case StyleSettings::eIdentity:
479 alignedColor = conservationColorer->
480 GetIdentityColor(block, seqIndex - block->GetRangeOfRow(row)->from);
481 break;
482 case StyleSettings::eVariety:
483 alignedColor = conservationColorer->
484 GetVarietyColor(block, seqIndex - block->GetRangeOfRow(row)->from);
485 break;
486 case StyleSettings::eWeightedVariety:
487 alignedColor = conservationColorer->
488 GetWeightedVarietyColor(block, seqIndex - block->GetRangeOfRow(row)->from);
489 break;
490 case StyleSettings::eInformationContent:
491 alignedColor = conservationColorer->
492 GetInformationContentColor(block, seqIndex - block->GetRangeOfRow(row)->from);
493 break;
494 case StyleSettings::eFit:
495 alignedColor = conservationColorer->
496 GetFitColor(block, seqIndex - block->GetRangeOfRow(row)->from, row);
497 break;
498 case StyleSettings::eBlockFit:
499 alignedColor = conservationColorer->GetBlockFitColor(block, row);
500 break;
501 case StyleSettings::eBlockZFit:
502 alignedColor = conservationColorer->GetBlockZFitColor(block, row);
503 break;
504 case StyleSettings::eBlockRowFit:
505 alignedColor = conservationColorer->GetBlockRowFitColor(block, row);
506 break;
507 default:
508 alignedColor = NULL;
509 }
510
511 return alignedColor;
512 }
513
GetAlignmentColor(const Sequence * sequence,unsigned int seqIndex,StyleSettings::eColorScheme colorScheme) const514 const Vector * BlockMultipleAlignment::GetAlignmentColor(
515 const Sequence *sequence, unsigned int seqIndex, StyleSettings::eColorScheme colorScheme) const
516 {
517 int row = GetRowForSequence(sequence);
518 if (row < 0) return NULL;
519 return GetAlignmentColor(row, seqIndex, colorScheme);
520 }
521
IsAligned(unsigned int row,unsigned int seqIndex) const522 bool BlockMultipleAlignment::IsAligned(unsigned int row, unsigned int seqIndex) const
523 {
524 const Block *block = GetBlock(row, seqIndex);
525 return (block && block->IsAligned());
526 }
527
GetBlock(unsigned int row,unsigned int seqIndex) const528 const Block * BlockMultipleAlignment::GetBlock(unsigned int row, unsigned int seqIndex) const
529 {
530 // make sure we're in range for this sequence
531 if (row >= NRows() || seqIndex >= (*sequences)[row]->Length()) {
532 ERRORMSG("BlockMultipleAlignment::GetBlock() - coordinate out of range");
533 return NULL;
534 }
535
536 const Block::Range *range;
537
538 // first check to see if it's in the same block as last time.
539 if (cachePrevBlock) {
540 range = cachePrevBlock->GetRangeOfRow(row);
541 if ((int)seqIndex >= range->from && (int)seqIndex <= range->to) return cachePrevBlock;
542 ++cacheBlockIterator; // start search at next block
543 } else {
544 cacheBlockIterator = blocks.begin();
545 }
546
547 // otherwise, perform block search. This search is most efficient when queries
548 // happen in order from left to right along a given row.
549 do {
550 if (cacheBlockIterator == blocks.end()) cacheBlockIterator = blocks.begin();
551 range = (*cacheBlockIterator)->GetRangeOfRow(row);
552 if ((int)seqIndex >= range->from && (int)seqIndex <= range->to) {
553 cachePrevBlock = *cacheBlockIterator; // cache this block
554 return cachePrevBlock;
555 }
556 ++cacheBlockIterator;
557 } while (1);
558 }
559
GetFirstAlignedBlockPosition(void) const560 int BlockMultipleAlignment::GetFirstAlignedBlockPosition(void) const
561 {
562 BlockList::const_iterator b = blocks.begin();
563 if (blocks.size() > 0 && (*b)->IsAligned()) // first block is aligned
564 return 0;
565 else if (blocks.size() >= 2 && (*(++b))->IsAligned()) // second block is aligned
566 return blocks.front()->width;
567 else
568 return -1;
569 }
570
GetAlignedDependentIndex(unsigned int masterSeqIndex,unsigned int dependentRow) const571 int BlockMultipleAlignment::GetAlignedDependentIndex(unsigned int masterSeqIndex, unsigned int dependentRow) const
572 {
573 const UngappedAlignedBlock
574 *aBlock = dynamic_cast<const UngappedAlignedBlock*>(GetBlock(0, masterSeqIndex));
575 if (!aBlock) return -1;
576
577 const Block::Range
578 *masterRange = aBlock->GetRangeOfRow(0),
579 *dependentRange = aBlock->GetRangeOfRow(dependentRow);
580 return (dependentRange->from + masterSeqIndex - masterRange->from);
581 }
582
SelectedRange(unsigned int row,unsigned int alnIndexFrom,unsigned int alnIndexTo,eUnalignedJustification justification,bool toggle) const583 void BlockMultipleAlignment::SelectedRange(unsigned int row, unsigned int alnIndexFrom, unsigned int alnIndexTo,
584 eUnalignedJustification justification, bool toggle) const
585 {
586 // translate from,to (alignment columns) into sequence indexes
587 const Sequence *sequence;
588 int fromSeqIndex, toSeqIndex;
589 bool ignored;
590
591 // trim selection area to size of this alignment
592 if (alnIndexTo >= totalWidth) alnIndexTo = totalWidth - 1;
593
594 // find first residue within range
595 while (alnIndexFrom <= alnIndexTo) {
596 GetSequenceAndIndexAt(alnIndexFrom, row, justification, &sequence, &fromSeqIndex, &ignored);
597 if (fromSeqIndex >= 0) break;
598 ++alnIndexFrom;
599 }
600 if (alnIndexFrom > alnIndexTo) return;
601
602 // find last residue within range
603 while (alnIndexTo >= alnIndexFrom) {
604 GetSequenceAndIndexAt(alnIndexTo, row, justification, &sequence, &toSeqIndex, &ignored);
605 if (toSeqIndex >= 0) break;
606 --alnIndexTo;
607 }
608
609 if (toggle)
610 GlobalMessenger()->ToggleHighlights(sequence, fromSeqIndex, toSeqIndex);
611 else
612 GlobalMessenger()->AddHighlights(sequence, fromSeqIndex, toSeqIndex);
613 }
614
SelectBlocks(unsigned int alnIndexFrom,unsigned int alnIndexTo,bool toggle) const615 void BlockMultipleAlignment::SelectBlocks(unsigned int alnIndexFrom, unsigned int alnIndexTo, bool toggle) const
616 {
617 typedef vector < pair < unsigned int, unsigned int > > VecPair;
618 VecPair highlightRanges;
619
620 BlockList::const_iterator b, be = blocks.end();
621 unsigned int start = 0, from, to;
622 for (b=blocks.begin(); b!=be; ++b) {
623 if ((*b)->IsAligned()) {
624 from = start;
625 to = start + (*b)->width - 1;
626 // highlight if the block intersects with the selection
627 if ((alnIndexFrom >= from && alnIndexFrom <= to) || (alnIndexTo >= from && alnIndexTo <= to)
628 || (alnIndexFrom < from && alnIndexTo > to))
629 highlightRanges.push_back(make_pair(from, to));
630 }
631 start += (*b)->width;
632 }
633
634 if (highlightRanges.size() == 0)
635 return;
636
637 // select all ranges in each row; justification is irrelevant since we're in aligned blocks only
638 for (unsigned int row=0; row<NRows(); ++row) {
639 VecPair::const_iterator p, pe = highlightRanges.end();
640 for (p=highlightRanges.begin(); p!=pe; ++p)
641 SelectedRange(row, p->first, p->second, eLeft, toggle);
642 }
643 }
644
GetAlignedBlockPosition(unsigned int alignmentIndex,int * blockColumn,int * blockWidth) const645 void BlockMultipleAlignment::GetAlignedBlockPosition(unsigned int alignmentIndex,
646 int *blockColumn, int *blockWidth) const
647 {
648 *blockColumn = *blockWidth = -1;
649 const BlockInfo& info = blockMap[alignmentIndex];
650 if (info.block->IsAligned()) {
651 *blockColumn = info.blockColumn;
652 *blockWidth = info.block->width;
653 }
654 }
655
GetBlockBefore(const Block * block) const656 Block * BlockMultipleAlignment::GetBlockBefore(const Block *block) const
657 {
658 Block *prevBlock = NULL;
659 BlockList::const_iterator b, be = blocks.end();
660 for (b=blocks.begin(); b!=be; ++b) {
661 if (*b == block) break;
662 prevBlock = *b;
663 }
664 return prevBlock;
665 }
666
GetUnalignedBlockBefore(const UngappedAlignedBlock * aBlock) const667 const UnalignedBlock * BlockMultipleAlignment::GetUnalignedBlockBefore(
668 const UngappedAlignedBlock *aBlock) const
669 {
670 const Block *prevBlock;
671 if (aBlock)
672 prevBlock = GetBlockBefore(aBlock);
673 else
674 prevBlock = blocks.back();
675 return dynamic_cast<const UnalignedBlock*>(prevBlock);
676 }
677
GetUnalignedBlockAfter(const UngappedAlignedBlock * aBlock) const678 const UnalignedBlock * BlockMultipleAlignment::GetUnalignedBlockAfter(
679 const UngappedAlignedBlock *aBlock) const
680 {
681 const Block *nextBlock;
682 if (aBlock)
683 nextBlock = GetBlockAfter(aBlock);
684 else
685 nextBlock = blocks.front();
686 return dynamic_cast<const UnalignedBlock*>(nextBlock);
687 }
688
GetBlockAfter(const Block * block) const689 Block * BlockMultipleAlignment::GetBlockAfter(const Block *block) const
690 {
691 BlockList::const_iterator b, be = blocks.end();
692 for (b=blocks.begin(); b!=be; ++b) {
693 if (*b == block) {
694 ++b;
695 if (b == be) break;
696 return *b;
697 }
698 }
699 return NULL;
700 }
701
InsertBlockBefore(Block * newBlock,const Block * insertAt)702 void BlockMultipleAlignment::InsertBlockBefore(Block *newBlock, const Block *insertAt)
703 {
704 BlockList::iterator b, be = blocks.end();
705 for (b=blocks.begin(); b!=be; ++b) {
706 if (*b == insertAt) {
707 blocks.insert(b, newBlock);
708 return;
709 }
710 }
711 WARNINGMSG("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
712 }
713
InsertBlockAfter(const Block * insertAt,Block * newBlock)714 void BlockMultipleAlignment::InsertBlockAfter(const Block *insertAt, Block *newBlock)
715 {
716 BlockList::iterator b, be = blocks.end();
717 for (b=blocks.begin(); b!=be; ++b) {
718 if (*b == insertAt) {
719 ++b;
720 blocks.insert(b, newBlock);
721 return;
722 }
723 }
724 WARNINGMSG("BlockMultipleAlignment::InsertBlockBefore() - couldn't find insertAt block");
725 }
726
RemoveBlock(Block * block)727 void BlockMultipleAlignment::RemoveBlock(Block *block)
728 {
729 BlockList::iterator b, be = blocks.end();
730 for (b=blocks.begin(); b!=be; ++b) {
731 if (*b == block) {
732 delete *b;
733 blocks.erase(b);
734 InitCache();
735 return;
736 }
737 }
738 WARNINGMSG("BlockMultipleAlignment::RemoveBlock() - couldn't find block");
739 }
740
MoveBlockBoundary(unsigned int columnFrom,unsigned int columnTo)741 bool BlockMultipleAlignment::MoveBlockBoundary(unsigned int columnFrom, unsigned int columnTo)
742 {
743 int blockColumn, blockWidth;
744 GetAlignedBlockPosition(columnFrom, &blockColumn, &blockWidth);
745 if (blockColumn < 0 || blockWidth <= 0) return false;
746
747 TRACEMSG("trying to move block boundary from " << columnFrom << " to " << columnTo);
748
749 const BlockInfo& info = blockMap[columnFrom];
750 unsigned int row;
751 int requestedShift = columnTo - columnFrom, actualShift = 0;
752 const Block::Range *range;
753
754 // shrink block from left
755 if (blockColumn == 0 && requestedShift > 0 && requestedShift < (int)info.block->width) {
756 actualShift = requestedShift;
757 TRACEMSG("shrinking block from left");
758 for (row=0; row<NRows(); ++row) {
759 range = info.block->GetRangeOfRow(row);
760 info.block->SetRangeOfRow(row, range->from + requestedShift, range->to);
761 }
762 info.block->width -= requestedShift;
763 Block *prevBlock = GetBlockBefore(info.block);
764 if (prevBlock && !prevBlock->IsAligned()) {
765 for (row=0; row<NRows(); ++row) {
766 range = prevBlock->GetRangeOfRow(row);
767 prevBlock->SetRangeOfRow(row, range->from, range->to + requestedShift);
768 }
769 prevBlock->width += requestedShift;
770 } else {
771 Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(prevBlock, info.block);
772 if (newUnalignedBlock) InsertBlockBefore(newUnalignedBlock, info.block);
773 TRACEMSG("added new unaligned block");
774 }
775 }
776
777 // shrink block from right (requestedShift < 0)
778 else if (blockColumn == (int)info.block->width - 1 &&
779 requestedShift < 0 && -requestedShift < (int)info.block->width) {
780 actualShift = requestedShift;
781 TRACEMSG("shrinking block from right");
782 for (row=0; row<NRows(); ++row) {
783 range = info.block->GetRangeOfRow(row);
784 info.block->SetRangeOfRow(row, range->from, range->to + requestedShift);
785 }
786 info.block->width += requestedShift;
787 Block *nextBlock = GetBlockAfter(info.block);
788 if (nextBlock && !nextBlock->IsAligned()) {
789 for (row=0; row<NRows(); ++row) {
790 range = nextBlock->GetRangeOfRow(row);
791 nextBlock->SetRangeOfRow(row, range->from + requestedShift, range->to);
792 }
793 nextBlock->width -= requestedShift;
794 } else {
795 Block *newUnalignedBlock = CreateNewUnalignedBlockBetween(info.block, nextBlock);
796 if (newUnalignedBlock) InsertBlockAfter(info.block, newUnalignedBlock);
797 TRACEMSG("added new unaligned block");
798 }
799 }
800
801 // grow block to right
802 else if (blockColumn == (int)info.block->width - 1 && requestedShift > 0) {
803 Block *nextBlock = GetBlockAfter(info.block);
804 if (nextBlock && !nextBlock->IsAligned()) {
805 int nRes;
806 actualShift = requestedShift;
807 for (row=0; row<NRows(); ++row) {
808 range = nextBlock->GetRangeOfRow(row);
809 nRes = range->to - range->from + 1;
810 if (nRes < actualShift) actualShift = nRes;
811 }
812 if (actualShift) {
813 TRACEMSG("growing block to right");
814 for (row=0; row<NRows(); ++row) {
815 range = info.block->GetRangeOfRow(row);
816 info.block->SetRangeOfRow(row, range->from, range->to + actualShift);
817 range = nextBlock->GetRangeOfRow(row);
818 nextBlock->SetRangeOfRow(row, range->from + actualShift, range->to);
819 }
820 info.block->width += actualShift;
821 nextBlock->width -= actualShift;
822 if (nextBlock->width == 0) {
823 RemoveBlock(nextBlock);
824 TRACEMSG("removed empty block");
825 }
826 }
827 }
828 }
829
830 // grow block to left (requestedShift < 0)
831 else if (blockColumn == 0 && requestedShift < 0) {
832 Block *prevBlock = GetBlockBefore(info.block);
833 if (prevBlock && !prevBlock->IsAligned()) {
834 int nRes;
835 actualShift = requestedShift;
836 for (row=0; row<NRows(); ++row) {
837 range = prevBlock->GetRangeOfRow(row);
838 nRes = range->to - range->from + 1;
839 if (nRes < -actualShift) actualShift = -nRes;
840 }
841 if (actualShift) {
842 TRACEMSG("growing block to left");
843 for (row=0; row<NRows(); ++row) {
844 range = info.block->GetRangeOfRow(row);
845 info.block->SetRangeOfRow(row, range->from + actualShift, range->to);
846 range = prevBlock->GetRangeOfRow(row);
847 prevBlock->SetRangeOfRow(row, range->from, range->to + actualShift);
848 }
849 info.block->width -= actualShift;
850 prevBlock->width += actualShift;
851 if (prevBlock->width == 0) {
852 RemoveBlock(prevBlock);
853 TRACEMSG("removed empty block");
854 }
855 }
856 }
857 }
858
859 if (actualShift != 0) {
860 UpdateBlockMapAndColors();
861 return true;
862 } else
863 return false;
864 }
865
ShiftRow(unsigned int row,unsigned int fromAlignmentIndex,unsigned int toAlignmentIndex,eUnalignedJustification justification)866 bool BlockMultipleAlignment::ShiftRow(unsigned int row, unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex,
867 eUnalignedJustification justification)
868 {
869 if (fromAlignmentIndex == toAlignmentIndex) return false;
870
871 Block
872 *blockFrom = blockMap[fromAlignmentIndex].block,
873 *blockTo = blockMap[toAlignmentIndex].block;
874
875 // at least one end of the drag must be in an aligned block
876 UngappedAlignedBlock *ABlock =
877 dynamic_cast<UngappedAlignedBlock*>(blockFrom);
878 if (ABlock) {
879 if (blockTo != blockFrom && blockTo->IsAligned()) return false;
880 } else {
881 ABlock = dynamic_cast<UngappedAlignedBlock*>(blockTo);
882 if (!ABlock) return false;
883 }
884
885 // and the other must be in the same aligned block or an adjacent unaligned block
886 UnalignedBlock
887 *prevUABlock = dynamic_cast<UnalignedBlock*>(GetBlockBefore(ABlock)),
888 *nextUABlock = dynamic_cast<UnalignedBlock*>(GetBlockAfter(ABlock));
889 if (blockFrom != blockTo &&
890 ((ABlock == blockFrom && prevUABlock != blockTo && nextUABlock != blockTo) ||
891 (ABlock == blockTo && prevUABlock != blockFrom && nextUABlock != blockFrom)))
892 return false;
893
894 int requestedShift, actualShift = 0, width = 0;
895
896 // slightly different behaviour when dragging from unaligned to aligned...
897 if (!blockFrom->IsAligned()) {
898 int fromSeqIndex, toSeqIndex;
899 GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, NULL, &fromSeqIndex, NULL);
900 GetSequenceAndIndexAt(toAlignmentIndex, row, justification, NULL, &toSeqIndex, NULL);
901 if (fromSeqIndex < 0 || toSeqIndex < 0) return false;
902 requestedShift = toSeqIndex - fromSeqIndex;
903 }
904
905 // vs. dragging from aligned
906 else {
907 requestedShift = toAlignmentIndex - fromAlignmentIndex;
908 }
909
910 const Block::Range *prevRange = NULL, *nextRange = NULL,
911 *range = ABlock->GetRangeOfRow(row);
912 if (prevUABlock) prevRange = prevUABlock->GetRangeOfRow(row);
913 if (nextUABlock) nextRange = nextUABlock->GetRangeOfRow(row);
914 if (requestedShift > 0) {
915 if (prevUABlock) width = prevRange->to - prevRange->from + 1;
916 actualShift = (width > requestedShift) ? requestedShift : width;
917 } else {
918 if (nextUABlock) width = nextRange->to - nextRange->from + 1;
919 actualShift = (width > -requestedShift) ? requestedShift : -width;
920 }
921 if (actualShift == 0) return false;
922
923 TRACEMSG("shifting row " << row << " by " << actualShift);
924
925 ABlock->SetRangeOfRow(row, range->from - actualShift, range->to - actualShift);
926
927 if (prevUABlock) {
928 prevUABlock->SetRangeOfRow(row, prevRange->from, prevRange->to - actualShift);
929 prevUABlock->width = 0;
930 for (unsigned int i=0; i<NRows(); ++i) {
931 prevRange = prevUABlock->GetRangeOfRow(i);
932 width = prevRange->to - prevRange->from + 1;
933 if (width < 0) ERRORMSG("BlockMultipleAlignment::ShiftRow() - negative width on left");
934 if (width > (int)prevUABlock->width) prevUABlock->width = width;
935 }
936 if (prevUABlock->width == 0) {
937 TRACEMSG("removing zero-width unaligned block on left");
938 RemoveBlock(prevUABlock);
939 }
940 } else {
941 TRACEMSG("creating unaligned block on left");
942 prevUABlock = CreateNewUnalignedBlockBetween(GetBlockBefore(ABlock), ABlock);
943 InsertBlockBefore(prevUABlock, ABlock);
944 }
945
946 if (nextUABlock) {
947 nextUABlock->SetRangeOfRow(row, nextRange->from - actualShift, nextRange->to);
948 nextUABlock->width = 0;
949 for (unsigned int i=0; i<NRows(); ++i) {
950 nextRange = nextUABlock->GetRangeOfRow(i);
951 width = nextRange->to - nextRange->from + 1;
952 if (width < 0) ERRORMSG("BlockMultipleAlignment::ShiftRow() - negative width on right");
953 if (width > (int)nextUABlock->width) nextUABlock->width = width;
954 }
955 if (nextUABlock->width == 0) {
956 TRACEMSG("removing zero-width unaligned block on right");
957 RemoveBlock(nextUABlock);
958 }
959 } else {
960 TRACEMSG("creating unaligned block on right");
961 nextUABlock = CreateNewUnalignedBlockBetween(ABlock, GetBlockAfter(ABlock));
962 InsertBlockAfter(ABlock, nextUABlock);
963 }
964
965 if (!CheckAlignedBlock(ABlock))
966 ERRORMSG("BlockMultipleAlignment::ShiftRow() - shift failed to create valid aligned block");
967
968 UpdateBlockMapAndColors();
969 return true;
970 }
971
ZipAlignResidue(unsigned int row,unsigned int alignmentIndex,bool moveRight,eUnalignedJustification justification)972 bool BlockMultipleAlignment::ZipAlignResidue(unsigned int row, unsigned int alignmentIndex, bool moveRight, eUnalignedJustification justification)
973 {
974 if (blockMap[alignmentIndex].block->IsAligned()) {
975 TRACEMSG("residue is already aligned");
976 return false;
977 }
978
979 UngappedAlignedBlock *aBlock = dynamic_cast<UngappedAlignedBlock*>(
980 moveRight ? GetBlockAfter(blockMap[alignmentIndex].block) : GetBlockBefore(blockMap[alignmentIndex].block));
981 if (!aBlock) {
982 TRACEMSG("no aligned block to the " << (moveRight ? "right" : "left"));
983 return false;
984 }
985
986 return ShiftRow(row, alignmentIndex,
987 GetAlignmentIndex(row,
988 (moveRight ? aBlock->GetRangeOfRow(row)->from : aBlock->GetRangeOfRow(row)->to),
989 justification),
990 justification);
991 }
992
OptimizeBlock(unsigned int row,unsigned int alignmentIndex,eUnalignedJustification justification)993 bool BlockMultipleAlignment::OptimizeBlock(unsigned int row, unsigned int alignmentIndex, eUnalignedJustification justification)
994 {
995 // is this an aligned block?
996 UngappedAlignedBlock *block = dynamic_cast<UngappedAlignedBlock*>(blockMap[alignmentIndex].block);
997 if (!block) {
998 TRACEMSG("block is unaligned");
999 return false;
1000 }
1001
1002 // see if there's any room to shift
1003 const UnalignedBlock *prevBlock = GetUnalignedBlockBefore(block), *nextBlock = GetUnalignedBlockAfter(block);
1004 unsigned int maxShiftRight = 0, maxShiftLeft = 0;
1005 const Block::Range *range;
1006 if (prevBlock) {
1007 range = prevBlock->GetRangeOfRow(row);
1008 maxShiftRight = range->to - range->from + 1;
1009 }
1010 if (nextBlock) {
1011 range = nextBlock->GetRangeOfRow(row);
1012 maxShiftLeft = range->to - range->from + 1;
1013 }
1014
1015 // if this is first or last block, constrain by footprint
1016 int blockNum = GetAlignedBlockNumber(alignmentIndex);
1017 if (blockNum == 1 || blockNum == (int)NAlignedBlocks()) {
1018 int excess = 0;
1019 if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
1020 WARNINGMSG("Can't get footprint excess residues from registry");
1021 if (blockNum == 1) {
1022 if ((int)maxShiftRight > excess) {
1023 WARNINGMSG("maxShiftRight constrained to " << excess << " by footprint excess residues preference");
1024 maxShiftRight = excess;
1025 }
1026 } else { // last block
1027 if ((int)maxShiftLeft > excess) {
1028 WARNINGMSG("maxShiftLeft constrained to " << excess << " by footprint excess residues preference");
1029 maxShiftLeft = excess;
1030 }
1031 }
1032 }
1033
1034 TRACEMSG("maxShiftRight " << maxShiftRight << ", maxShiftLeft " << maxShiftLeft);
1035 if (maxShiftRight == 0 && maxShiftLeft == 0) {
1036 TRACEMSG("no room to shift this block");
1037 return false;
1038 }
1039
1040 // scan all possible block positions, find max information content
1041 range = block->GetRangeOfRow(row);
1042 unsigned int bestSeqIndexStart = range->from;
1043 float score, maxScore = kMin_Float;
1044 for (unsigned int seqIndexStart = range->from - maxShiftRight; seqIndexStart <= range->from + maxShiftLeft; ++seqIndexStart) {
1045
1046 // calculate block's info content given each position of this row
1047 score = 0.0f;
1048 for (unsigned int blockColumn=0; blockColumn<block->width; ++blockColumn) {
1049
1050 // create profile for this column
1051 typedef std::map < char, unsigned int > ColumnProfile;
1052 ColumnProfile profile;
1053 ColumnProfile::iterator p, pe;
1054 for (unsigned int r=0; r<NRows(); ++r) {
1055 unsigned int seqIndex = (r == row) ? (seqIndexStart + blockColumn) : (block->GetRangeOfRow(r)->from + blockColumn);
1056 char ch = ScreenResidueCharacter(GetSequenceOfRow(r)->sequenceString[seqIndex]);
1057 if ((p=profile.find(ch)) != profile.end())
1058 ++(p->second);
1059 else
1060 profile[ch] = 1;
1061 }
1062
1063 // information content (calculated in bits -> logs of base 2) for this column
1064 for (p=profile.begin(), pe=profile.end(); p!=pe; ++p) {
1065 static const float ln2 = log(2.0f), threshhold = 0.0001f;
1066 float expFreq = GetStandardProbability(p->first);
1067 if (expFreq > threshhold) {
1068 float obsFreq = 1.0f * p->second / NRows(),
1069 freqRatio = obsFreq / expFreq;
1070 if (freqRatio > threshhold)
1071 score += obsFreq * ((float) log(freqRatio)) / ln2; // information content
1072 }
1073 }
1074 }
1075
1076 // keep track of best position
1077 if (seqIndexStart == range->from - maxShiftRight || score > maxScore) {
1078 maxScore = score;
1079 bestSeqIndexStart = seqIndexStart;
1080 }
1081 }
1082
1083 // if the best position is other than the current, then shift the block accordingly
1084 bool moved = ((int)bestSeqIndexStart != range->from);
1085 if (moved) {
1086 if ((int)bestSeqIndexStart < range->from)
1087 TRACEMSG("shifting block right by " << (range->from - bestSeqIndexStart));
1088 else
1089 TRACEMSG("shifting block left by " << (bestSeqIndexStart - range->from));
1090 alignmentIndex = GetAlignmentIndex(row, range->from, justification);
1091 unsigned int alnIdx2 = GetAlignmentIndex(row, bestSeqIndexStart, justification);
1092 moved = ShiftRow(row, alnIdx2, alignmentIndex, justification);
1093 if (!moved)
1094 ERRORMSG("ShiftRow() failed!");
1095 } else
1096 TRACEMSG("block was not moved");
1097
1098 return moved;
1099 }
1100
SplitBlock(unsigned int alignmentIndex)1101 bool BlockMultipleAlignment::SplitBlock(unsigned int alignmentIndex)
1102 {
1103 const BlockInfo& info = blockMap[alignmentIndex];
1104 if (!info.block->IsAligned() || info.block->width < 2 || info.blockColumn == 0)
1105 return false;
1106
1107 TRACEMSG("splitting block");
1108 UngappedAlignedBlock *newAlignedBlock = new UngappedAlignedBlock(this);
1109 newAlignedBlock->width = info.block->width - info.blockColumn;
1110 info.block->width = info.blockColumn;
1111
1112 const Block::Range *range;
1113 unsigned int oldTo;
1114 for (unsigned int row=0; row<NRows(); ++row) {
1115 range = info.block->GetRangeOfRow(row);
1116 oldTo = range->to;
1117 info.block->SetRangeOfRow(row, range->from, range->from + info.block->width - 1);
1118 newAlignedBlock->SetRangeOfRow(row, oldTo - newAlignedBlock->width + 1, oldTo);
1119 }
1120
1121 InsertBlockAfter(info.block, newAlignedBlock);
1122 if (!CheckAlignedBlock(info.block) || !CheckAlignedBlock(newAlignedBlock))
1123 ERRORMSG("BlockMultipleAlignment::SplitBlock() - split failed to create valid blocks");
1124
1125 UpdateBlockMapAndColors();
1126 return true;
1127 }
1128
MergeBlocks(unsigned int fromAlignmentIndex,unsigned int toAlignmentIndex)1129 bool BlockMultipleAlignment::MergeBlocks(unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex)
1130 {
1131 Block
1132 *expandedBlock = blockMap[fromAlignmentIndex].block,
1133 *lastBlock = blockMap[toAlignmentIndex].block;
1134 if (expandedBlock == lastBlock) return false;
1135 unsigned int i;
1136 for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i)
1137 if (!blockMap[i].block->IsAligned()) return false;
1138
1139 TRACEMSG("merging block(s)");
1140 for (i=0; i<NRows(); ++i)
1141 expandedBlock->SetRangeOfRow(i,
1142 expandedBlock->GetRangeOfRow(i)->from, lastBlock->GetRangeOfRow(i)->to);
1143 expandedBlock->width = lastBlock->GetRangeOfRow(0)->to - expandedBlock->GetRangeOfRow(0)->from + 1;
1144
1145 Block *deletedBlock = NULL, *blockToDelete;
1146 for (i=fromAlignmentIndex; i<=toAlignmentIndex; ++i) {
1147 blockToDelete = blockMap[i].block;
1148 if (blockToDelete == expandedBlock) continue;
1149 if (blockToDelete != deletedBlock) {
1150 deletedBlock = blockToDelete;
1151 RemoveBlock(blockToDelete);
1152 }
1153 }
1154
1155 if (!CheckAlignedBlock(expandedBlock))
1156 ERRORMSG("BlockMultipleAlignment::MergeBlocks() - merge failed to create valid block");
1157
1158 UpdateBlockMapAndColors();
1159 return true;
1160 }
1161
CreateBlock(unsigned int fromAlignmentIndex,unsigned int toAlignmentIndex,eUnalignedJustification justification)1162 bool BlockMultipleAlignment::CreateBlock(unsigned int fromAlignmentIndex, unsigned int toAlignmentIndex,
1163 eUnalignedJustification justification)
1164 {
1165 const BlockInfo& info = blockMap[fromAlignmentIndex];
1166 UnalignedBlock *prevUABlock = dynamic_cast<UnalignedBlock*>(info.block);
1167 if (!prevUABlock || info.block != blockMap[toAlignmentIndex].block) return false;
1168 int seqIndexFrom, seqIndexTo;
1169 unsigned int row,
1170 newBlockWidth = toAlignmentIndex - fromAlignmentIndex + 1,
1171 origWidth = prevUABlock->width;
1172 vector < unsigned int > seqIndexesFrom(NRows()), seqIndexesTo(NRows());
1173 const Sequence *seq;
1174 bool ignored;
1175 for (row=0; row<NRows(); ++row) {
1176 if (!GetSequenceAndIndexAt(fromAlignmentIndex, row, justification, &seq, &seqIndexFrom, &ignored) ||
1177 !GetSequenceAndIndexAt(toAlignmentIndex, row, justification, &seq, &seqIndexTo, &ignored) ||
1178 seqIndexFrom < 0 || seqIndexTo < 0 ||
1179 seqIndexTo - seqIndexFrom + 1 != (int)newBlockWidth)
1180 return false;
1181 seqIndexesFrom[row] = seqIndexFrom;
1182 seqIndexesTo[row] = seqIndexTo;
1183 }
1184
1185 TRACEMSG("creating new aligned and unaligned blocks");
1186
1187 UnalignedBlock *nextUABlock = new UnalignedBlock(this);
1188 UngappedAlignedBlock *ABlock = new UngappedAlignedBlock(this);
1189 prevUABlock->width = nextUABlock->width = 0;
1190
1191 bool deletePrevUABlock = true, deleteNextUABlock = true;
1192 const Block::Range *prevRange;
1193 int rangeWidth;
1194 for (row=0; row<NRows(); ++row) {
1195 prevRange = prevUABlock->GetRangeOfRow(row);
1196
1197 nextUABlock->SetRangeOfRow(row, seqIndexesTo[row] + 1, prevRange->to);
1198 rangeWidth = prevRange->to - seqIndexesTo[row];
1199 if (rangeWidth < 0)
1200 ERRORMSG("BlockMultipleAlignment::CreateBlock() - negative nextRange width");
1201 else if (rangeWidth > 0) {
1202 if (rangeWidth > (int)nextUABlock->width) nextUABlock->width = rangeWidth;
1203 deleteNextUABlock = false;
1204 }
1205
1206 prevUABlock->SetRangeOfRow(row, prevRange->from, seqIndexesFrom[row] - 1);
1207 rangeWidth = seqIndexesFrom[row] - prevRange->from;
1208 if (rangeWidth < 0)
1209 ERRORMSG("BlockMultipleAlignment::CreateBlock() - negative prevRange width");
1210 else if (rangeWidth > 0) {
1211 if (rangeWidth > (int)prevUABlock->width) prevUABlock->width = rangeWidth;
1212 deletePrevUABlock = false;
1213 }
1214
1215 ABlock->SetRangeOfRow(row, seqIndexesFrom[row], seqIndexesTo[row]);
1216 }
1217
1218 ABlock->width = newBlockWidth;
1219 if (prevUABlock->width + ABlock->width + nextUABlock->width != origWidth)
1220 ERRORMSG("BlockMultipleAlignment::CreateBlock() - bad block widths sum");
1221
1222 InsertBlockAfter(prevUABlock, ABlock);
1223 InsertBlockAfter(ABlock, nextUABlock);
1224 if (deletePrevUABlock) {
1225 TRACEMSG("deleting zero-width unaligned block on left");
1226 RemoveBlock(prevUABlock);
1227 }
1228 if (deleteNextUABlock) {
1229 TRACEMSG("deleting zero-width unaligned block on right");
1230 RemoveBlock(nextUABlock);
1231 }
1232
1233 if (!CheckAlignedBlock(ABlock))
1234 ERRORMSG("BlockMultipleAlignment::CreateBlock() - failed to create valid block");
1235
1236 UpdateBlockMapAndColors();
1237 return true;
1238 }
1239
DeleteBlock(unsigned int alignmentIndex)1240 bool BlockMultipleAlignment::DeleteBlock(unsigned int alignmentIndex)
1241 {
1242 Block *block = blockMap[alignmentIndex].block;
1243 if (!block || !block->IsAligned()) return false;
1244
1245 TRACEMSG("deleting block");
1246 Block
1247 *prevBlock = GetBlockBefore(block),
1248 *nextBlock = GetBlockAfter(block);
1249
1250 // unaligned blocks on both sides - note that total alignment width can change!
1251 if (prevBlock && !prevBlock->IsAligned() && nextBlock && !nextBlock->IsAligned()) {
1252 const Block::Range *prevRange, *nextRange;
1253 unsigned int maxWidth = 0, width;
1254 for (unsigned int row=0; row<NRows(); ++row) {
1255 prevRange = prevBlock->GetRangeOfRow(row);
1256 nextRange = nextBlock->GetRangeOfRow(row);
1257 width = nextRange->to - prevRange->from + 1;
1258 prevBlock->SetRangeOfRow(row, prevRange->from, nextRange->to);
1259 if (width > maxWidth) maxWidth = width;
1260 }
1261 prevBlock->width = maxWidth;
1262 TRACEMSG("removing extra unaligned block");
1263 RemoveBlock(nextBlock);
1264 }
1265
1266 // unaligned block on left only
1267 else if (prevBlock && !prevBlock->IsAligned()) {
1268 const Block::Range *prevRange, *range;
1269 for (unsigned int row=0; row<NRows(); ++row) {
1270 prevRange = prevBlock->GetRangeOfRow(row);
1271 range = block->GetRangeOfRow(row);
1272 prevBlock->SetRangeOfRow(row, prevRange->from, range->to);
1273 }
1274 prevBlock->width += block->width;
1275 }
1276
1277 // unaligned block on right only
1278 else if (nextBlock && !nextBlock->IsAligned()) {
1279 const Block::Range *range, *nextRange;
1280 for (unsigned int row=0; row<NRows(); ++row) {
1281 range = block->GetRangeOfRow(row);
1282 nextRange = nextBlock->GetRangeOfRow(row);
1283 nextBlock->SetRangeOfRow(row, range->from, nextRange->to);
1284 }
1285 nextBlock->width += block->width;
1286 }
1287
1288 // no adjacent unaligned blocks
1289 else {
1290 TRACEMSG("creating new unaligned block");
1291 Block *newBlock = CreateNewUnalignedBlockBetween(prevBlock, nextBlock);
1292 if (prevBlock)
1293 InsertBlockAfter(prevBlock, newBlock);
1294 else if (nextBlock)
1295 InsertBlockBefore(newBlock, nextBlock);
1296 else
1297 blocks.push_back(newBlock);
1298 }
1299
1300 RemoveBlock(block);
1301 UpdateBlockMapAndColors();
1302 return true;
1303 }
1304
DeleteAllBlocks(void)1305 bool BlockMultipleAlignment::DeleteAllBlocks(void)
1306 {
1307 if (blocks.size() == 0) return false;
1308
1309 DELETE_ALL_AND_CLEAR(blocks, BlockList);
1310 InitCache();
1311 AddUnalignedBlocks(); // one single unaligned block for whole alignment
1312 UpdateBlockMapAndColors();
1313 return true;
1314 }
1315
DeleteRow(unsigned int row)1316 bool BlockMultipleAlignment::DeleteRow(unsigned int row)
1317 {
1318 if (row >= NRows()) {
1319 ERRORMSG("BlockMultipleAlignment::DeleteRow() - row out of range");
1320 return false;
1321 }
1322
1323 // remove sequence from list
1324 vector < bool > removeRows(NRows(), false);
1325 removeRows[row] = true;
1326 VectorRemoveElements(*sequences, removeRows, 1);
1327 VectorRemoveElements(rowDoubles, removeRows, 1);
1328 VectorRemoveElements(rowStrings, removeRows, 1);
1329 VectorRemoveElements(geometryViolations, removeRows, 1);
1330
1331 // delete row from all blocks, removing any zero-width blocks
1332 BlockList::iterator b = blocks.begin(), br, be = blocks.end();
1333 while (b != be) {
1334 (*b)->DeleteRow(row);
1335 if ((*b)->width == 0) {
1336 br = b;
1337 ++b;
1338 TRACEMSG("deleting block resized to zero width");
1339 RemoveBlock(*br);
1340 } else
1341 ++b;
1342 }
1343
1344 // update total alignment width
1345 UpdateBlockMapAndColors();
1346 InitCache();
1347
1348 return true;
1349 }
1350
GetBlocks(ConstBlockList * cbl) const1351 void BlockMultipleAlignment::GetBlocks(ConstBlockList *cbl) const
1352 {
1353 cbl->clear();
1354 cbl->reserve(blocks.size());
1355 BlockList::const_iterator b, be = blocks.end();
1356 for (b=blocks.begin(); b!=be; ++b)
1357 cbl->push_back(*b);
1358 }
1359
GetUngappedAlignedBlocks(UngappedAlignedBlockList * uabs) const1360 unsigned int BlockMultipleAlignment::GetUngappedAlignedBlocks(UngappedAlignedBlockList *uabs) const
1361 {
1362 uabs->clear();
1363 uabs->reserve(blocks.size());
1364 BlockList::const_iterator b, be = blocks.end();
1365 for (b=blocks.begin(); b!=be; ++b) {
1366 UngappedAlignedBlock *uab = dynamic_cast<UngappedAlignedBlock*>(*b);
1367 if (uab) uabs->push_back(uab);
1368 }
1369 uabs->resize(uabs->size());
1370 return uabs->size();
1371 }
1372
ExtractRows(const vector<unsigned int> & dependentsToRemove,AlignmentList * pairwiseAlignments)1373 bool BlockMultipleAlignment::ExtractRows(
1374 const vector < unsigned int >& dependentsToRemove, AlignmentList *pairwiseAlignments)
1375 {
1376 if (dependentsToRemove.size() == 0) return false;
1377
1378 // make a bool list of rows to remove, also checking to make sure dependent list items are in range
1379 unsigned int i;
1380 vector < bool > removeRows(NRows(), false);
1381 for (i=0; i<dependentsToRemove.size(); ++i) {
1382 if (dependentsToRemove[i] > 0 && dependentsToRemove[i] < NRows()) {
1383 removeRows[dependentsToRemove[i]] = true;
1384 } else {
1385 ERRORMSG("BlockMultipleAlignment::ExtractRows() - can't extract row "
1386 << dependentsToRemove[i]);
1387 return false;
1388 }
1389 }
1390
1391 if (pairwiseAlignments) {
1392 TRACEMSG("creating new pairwise alignments");
1393 SetDiagPostLevel(eDiag_Warning); // otherwise, info messages take a long time if lots of rows
1394
1395 UngappedAlignedBlockList uaBlocks;
1396 GetUngappedAlignedBlocks(&uaBlocks);
1397 UngappedAlignedBlockList::const_iterator u, ue = uaBlocks.end();
1398
1399 for (i=0; i<dependentsToRemove.size(); ++i) {
1400
1401 // redraw molecule associated with removed row
1402 const Molecule *molecule = GetSequenceOfRow(dependentsToRemove[i])->molecule;
1403 if (molecule) GlobalMessenger()->PostRedrawMolecule(molecule);
1404
1405 // create new pairwise alignment from each removed row
1406 SequenceList *newSeqs = new SequenceList(2);
1407 (*newSeqs)[0] = (*sequences)[0];
1408 (*newSeqs)[1] = (*sequences)[dependentsToRemove[i]];
1409 BlockMultipleAlignment *newAlignment = new BlockMultipleAlignment(newSeqs, alignmentManager);
1410 for (u=uaBlocks.begin(); u!=ue; ++u) {
1411 // only copy blocks that aren't flagged to be realigned
1412 if (markBlocks.find(*u) == markBlocks.end()) {
1413 UngappedAlignedBlock *newABlock = new UngappedAlignedBlock(newAlignment);
1414 const Block::Range *range = (*u)->GetRangeOfRow(0);
1415 newABlock->SetRangeOfRow(0, range->from, range->to);
1416 range = (*u)->GetRangeOfRow(dependentsToRemove[i]);
1417 newABlock->SetRangeOfRow(1, range->from, range->to);
1418 newABlock->width = range->to - range->from + 1;
1419 newAlignment->AddAlignedBlockAtEnd(newABlock);
1420 }
1421 }
1422 if (!newAlignment->AddUnalignedBlocks() ||
1423 !newAlignment->UpdateBlockMapAndColors()) {
1424 ERRORMSG("BlockMultipleAlignment::ExtractRows() - error creating new alignment");
1425 return false;
1426 }
1427
1428 // add aligned region info (for threader to use later on)
1429 if (uaBlocks.size() > 0) {
1430 int excess = 0;
1431 if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
1432 WARNINGMSG("Can't get footprint excess residues from registry");
1433 newAlignment->alignDependentFrom =
1434 uaBlocks.front()->GetRangeOfRow(dependentsToRemove[i])->from - excess;
1435 if (newAlignment->alignDependentFrom < 0)
1436 newAlignment->alignDependentFrom = 0;
1437 newAlignment->alignDependentTo =
1438 uaBlocks.back()->GetRangeOfRow(dependentsToRemove[i])->to + excess;
1439 if (newAlignment->alignDependentTo >= (int)((*newSeqs)[1]->Length()))
1440 newAlignment->alignDependentTo = (*newSeqs)[1]->Length() - 1;
1441 TRACEMSG((*newSeqs)[1]->identifier->ToString() << " aligned from "
1442 << newAlignment->alignDependentFrom << " to " << newAlignment->alignDependentTo);
1443 }
1444
1445 pairwiseAlignments->push_back(newAlignment);
1446 }
1447 SetDiagPostLevel(eDiag_Info);
1448 }
1449
1450 // remove sequences
1451 TRACEMSG("deleting sequences");
1452 VectorRemoveElements(*sequences, removeRows, dependentsToRemove.size());
1453 VectorRemoveElements(rowDoubles, removeRows, dependentsToRemove.size());
1454 VectorRemoveElements(rowStrings, removeRows, dependentsToRemove.size());
1455 VectorRemoveElements(geometryViolations, removeRows, dependentsToRemove.size());
1456
1457 // delete row from all blocks, removing any zero-width blocks
1458 TRACEMSG("deleting alignment rows from blocks");
1459 BlockList::const_iterator b = blocks.begin(), br, be = blocks.end();
1460 while (b != be) {
1461 (*b)->DeleteRows(removeRows, dependentsToRemove.size());
1462 if ((*b)->width == 0) {
1463 br = b;
1464 ++b;
1465 TRACEMSG("deleting block resized to zero width");
1466 RemoveBlock(*br);
1467 } else
1468 ++b;
1469 }
1470
1471 // update total alignment width
1472 UpdateBlockMapAndColors();
1473 InitCache();
1474 return true;
1475 }
1476
MergeAlignment(const BlockMultipleAlignment * newAlignment)1477 bool BlockMultipleAlignment::MergeAlignment(const BlockMultipleAlignment *newAlignment)
1478 {
1479 // check to see if the alignment is compatible - must have same master sequence
1480 // and blocks of new alignment must contain blocks of this alignment; at same time,
1481 // build up map of aligned blocks in new alignment that correspond to aligned blocks
1482 // of this object, for convenient lookup later
1483 if (newAlignment->GetMaster() != GetMaster()) return false;
1484
1485 const Block::Range *newRange, *thisRange;
1486 BlockList::const_iterator nb, nbe = newAlignment->blocks.end();
1487 BlockList::iterator b, be = blocks.end();
1488 typedef map < UngappedAlignedBlock *, const UngappedAlignedBlock * > AlignedBlockMap;
1489 AlignedBlockMap correspondingNewBlocks;
1490
1491 for (b=blocks.begin(); b!=be; ++b) {
1492 if (!(*b)->IsAligned()) continue;
1493 for (nb=newAlignment->blocks.begin(); nb!=nbe; ++nb) {
1494 if (!(*nb)->IsAligned()) continue;
1495
1496 newRange = (*nb)->GetRangeOfRow(0);
1497 thisRange = (*b)->GetRangeOfRow(0);
1498 if (newRange->from <= thisRange->from && newRange->to >= thisRange->to) {
1499 correspondingNewBlocks[dynamic_cast<UngappedAlignedBlock*>(*b)] =
1500 dynamic_cast<const UngappedAlignedBlock*>(*nb);
1501 break;
1502 }
1503 }
1504 if (nb == nbe) return false; // no corresponding block found
1505 }
1506
1507 // add dependent sequences from new alignment; also copy other row-associated info
1508 unsigned int i, nNewRows = newAlignment->sequences->size() - 1;
1509 sequences->resize(sequences->size() + nNewRows);
1510 rowDoubles.resize(rowDoubles.size() + nNewRows);
1511 rowStrings.resize(rowStrings.size() + nNewRows);
1512 geometryViolations.resize(geometryViolations.size() + nNewRows);
1513 for (i=0; i<nNewRows; ++i) {
1514 unsigned int j = NRows() + i - nNewRows;
1515 (*sequences)[j] = (*(newAlignment->sequences))[i + 1];
1516 SetRowDouble(j, newAlignment->GetRowDouble(i + 1));
1517 SetRowStatusLine(j, newAlignment->GetRowStatusLine(i + 1));
1518 geometryViolations[j] = newAlignment->geometryViolations[i + 1];
1519 }
1520
1521 // now that we know blocks are compatible, add new rows at end of this alignment, containing
1522 // all rows (except master) from new alignment; only that part of the aligned blocks from
1523 // the new alignment that intersect with the aligned blocks from this object are added, so
1524 // that this object's block structure is unchanged
1525
1526 // resize all aligned blocks
1527 for (b=blocks.begin(); b!=be; ++b) (*b)->AddRows(nNewRows);
1528
1529 // set ranges of aligned blocks, and add them to the list
1530 AlignedBlockMap::const_iterator ab, abe = correspondingNewBlocks.end();
1531 const Block::Range *thisMaster, *newMaster;
1532 for (ab=correspondingNewBlocks.begin(); ab!=abe; ++ab) {
1533 thisMaster = ab->first->GetRangeOfRow(0);
1534 newMaster = ab->second->GetRangeOfRow(0);
1535 for (i=0; i<nNewRows; ++i) {
1536 newRange = ab->second->GetRangeOfRow(i + 1);
1537 ab->first->SetRangeOfRow(NRows() + i - nNewRows,
1538 newRange->from + thisMaster->from - newMaster->from,
1539 newRange->to + thisMaster->to - newMaster->to);
1540 }
1541 }
1542
1543 // delete then recreate the unaligned blocks, in case a merge requires the
1544 // creation of a new unaligned block
1545 for (b=blocks.begin(); b!=be; ) {
1546 if (!(*b)->IsAligned()) {
1547 BlockList::iterator bb(b);
1548 ++bb;
1549 delete *b;
1550 blocks.erase(b);
1551 b = bb;
1552 } else
1553 ++b;
1554 }
1555 InitCache();
1556
1557 // update this alignment, but leave row scores/status alone
1558 if (!AddUnalignedBlocks() || !UpdateBlockMapAndColors(false)) {
1559 ERRORMSG("BlockMultipleAlignment::MergeAlignment() - internal update after merge failed");
1560 return false;
1561 }
1562 return true;
1563 }
1564
ShowGeometryViolations(bool showGV)1565 unsigned int BlockMultipleAlignment::ShowGeometryViolations(bool showGV)
1566 {
1567 geometryViolations.clear();
1568 geometryViolations.resize(NRows());
1569
1570 if (!showGV || !GetMaster()->molecule || GetMaster()->molecule->parentSet->isAlphaOnly) {
1571 showGeometryViolations = false;
1572 return 0;
1573 }
1574
1575 Threader::GeometryViolationsForRow violations;
1576 unsigned int nViolations = alignmentManager->threader->GetGeometryViolations(this, &violations);
1577 if (violations.size() != NRows()) {
1578 ERRORMSG("BlockMultipleAlignment::ShowGeometryViolations() - wrong size list");
1579 showGeometryViolations = false;
1580 return 0;
1581 }
1582
1583 for (unsigned int row=0; row<NRows(); ++row) {
1584 geometryViolations[row].resize(GetSequenceOfRow(row)->Length(), false);
1585 Threader::IntervalList::const_iterator i, ie = violations[row].end();
1586 for (i=violations[row].begin(); i!=ie; ++i)
1587 for (unsigned int l=i->first; l<=i->second; ++l)
1588 geometryViolations[row][l] = true;
1589 }
1590
1591 showGeometryViolations = true;
1592 return nViolations;
1593 }
1594
CreatePairwiseSeqAlignFromMultipleRow(const BlockMultipleAlignment * multiple,const BlockMultipleAlignment::UngappedAlignedBlockList & blocks,unsigned int dependentRow)1595 CSeq_align * CreatePairwiseSeqAlignFromMultipleRow(const BlockMultipleAlignment *multiple,
1596 const BlockMultipleAlignment::UngappedAlignedBlockList& blocks, unsigned int dependentRow)
1597 {
1598 if (!multiple || dependentRow >= multiple->NRows()) {
1599 ERRORMSG("CreatePairwiseSeqAlignFromMultipleRow() - bad parameters");
1600 return NULL;
1601 }
1602
1603 CSeq_align *seqAlign = new CSeq_align();
1604 seqAlign->SetType(CSeq_align::eType_partial);
1605 seqAlign->SetDim(2);
1606
1607 CSeq_align::C_Segs::TDendiag& denDiags = seqAlign->SetSegs().SetDendiag();
1608 denDiags.resize((blocks.size() > 0) ? blocks.size() : 1);
1609
1610 CSeq_align::C_Segs::TDendiag::iterator d, de = denDiags.end();
1611 BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b = blocks.begin();
1612 const Block::Range *range;
1613 for (d=denDiags.begin(); d!=de; ++d, ++b) {
1614
1615 CDense_diag *denDiag = new CDense_diag();
1616 d->Reset(denDiag);
1617 denDiag->SetDim(2);
1618 denDiag->SetIds().resize(2);
1619
1620 // master row
1621 denDiag->SetIds().front().Reset(multiple->GetSequenceOfRow(0)->CreateSeqId());
1622 if (blocks.size() > 0) {
1623 range = (*b)->GetRangeOfRow(0);
1624 denDiag->SetStarts().push_back(range->from);
1625 } else
1626 denDiag->SetStarts().push_back(0);
1627
1628 // dependent row
1629 denDiag->SetIds().back().Reset(multiple->GetSequenceOfRow(dependentRow)->CreateSeqId());
1630 if (blocks.size() > 0) {
1631 range = (*b)->GetRangeOfRow(dependentRow);
1632 denDiag->SetStarts().push_back(range->from);
1633 } else
1634 denDiag->SetStarts().push_back(0);
1635
1636 // block width
1637 denDiag->SetLen((blocks.size() > 0) ? (*b)->width : 0);
1638 }
1639
1640 return seqAlign;
1641 }
1642
MarkBlock(unsigned int column)1643 bool BlockMultipleAlignment::MarkBlock(unsigned int column)
1644 {
1645 if (column < blockMap.size() && blockMap[column].block->IsAligned()) {
1646 TRACEMSG("marked block for realignment");
1647 markBlocks[blockMap[column].block] = true;
1648 return true;
1649 }
1650 return false;
1651 }
1652
ClearMarks(void)1653 bool BlockMultipleAlignment::ClearMarks(void)
1654 {
1655 if (markBlocks.size() == 0) return false;
1656 markBlocks.clear();
1657 return true;
1658 }
1659
GetMarkedBlockNumbers(std::vector<unsigned int> * alignedBlockNumbers) const1660 void BlockMultipleAlignment::GetMarkedBlockNumbers(std::vector < unsigned int > *alignedBlockNumbers) const
1661 {
1662 alignedBlockNumbers->resize(markBlocks.size());
1663 unsigned int a = 0, i = 0;
1664 BlockList::const_iterator b, be = blocks.end();
1665 for (b=blocks.begin(); b!=be; ++b) {
1666 if ((*b)->IsAligned()) {
1667 if (markBlocks.find(*b) != markBlocks.end())
1668 (*alignedBlockNumbers)[i++] = a;
1669 ++a;
1670 }
1671 }
1672 }
1673
HighlightAlignedColumnsOfMasterRange(unsigned int from,unsigned int to) const1674 bool BlockMultipleAlignment::HighlightAlignedColumnsOfMasterRange(unsigned int from, unsigned int to) const
1675 {
1676 const Sequence *master = GetMaster();
1677
1678 // must do one column at a time, rather than range, in case there are inserts wrt the master
1679 bool anyError = false;
1680 for (unsigned int i=from; i<=to; ++i) {
1681
1682 // sanity check
1683 if (i >= master->Length() || !IsAligned(0U, i)) {
1684 WARNINGMSG("Can't highlight alignment at master residue " << (i+1));
1685 anyError = true;
1686 // highlight unaligned residues, but master only
1687 if (i < master->Length())
1688 GlobalMessenger()->AddHighlights(GetSequenceOfRow(0), i, i);
1689 continue;
1690 }
1691
1692 // get block and offset
1693 const Block *block = GetBlock(0, i);
1694 unsigned int blockOffset = i - block->GetRangeOfRow(0)->from;
1695
1696 // highlight aligned residue in each row
1697 for (unsigned int row=0; row<NRows(); ++row) {
1698 unsigned int dependentIndex = block->GetRangeOfRow(row)->from + blockOffset;
1699 GlobalMessenger()->AddHighlights(GetSequenceOfRow(row), dependentIndex, dependentIndex);
1700 }
1701 }
1702
1703 return !anyError;
1704 }
1705
NAlignedBlocks(void) const1706 unsigned int BlockMultipleAlignment::NAlignedBlocks(void) const
1707 {
1708 unsigned int n = 0;
1709 BlockList::const_iterator b, be = blocks.end();
1710 for (b=blocks.begin(); b!=be; ++b) if ((*b)->IsAligned()) ++n;
1711 return n;
1712 }
1713
SumOfAlignedBlockWidths(void) const1714 unsigned int BlockMultipleAlignment::SumOfAlignedBlockWidths(void) const
1715 {
1716 unsigned int w = 0;
1717 BlockList::const_iterator b, be = blocks.end();
1718 for (b=blocks.begin(); b!=be; ++b)
1719 if ((*b)->IsAligned())
1720 w += (*b)->width;
1721 return w;
1722 }
1723
GetRelativeAlignmentFraction(unsigned int alignmentIndex) const1724 double BlockMultipleAlignment::GetRelativeAlignmentFraction(unsigned int alignmentIndex) const
1725 {
1726 if (!blockMap[alignmentIndex].block->IsAligned())
1727 return -1.0;
1728
1729 unsigned int f = 0;
1730 BlockList::const_iterator b, be = blocks.end();
1731 for (b=blocks.begin(); b!=be; ++b) {
1732 if ((*b)->IsAligned()) {
1733 if (*b == blockMap[alignmentIndex].block) {
1734 f += blockMap[alignmentIndex].blockColumn;
1735 return (((double) f) / (SumOfAlignedBlockWidths() - 1));
1736 } else {
1737 f += (*b)->width;
1738 }
1739 }
1740 }
1741 return -1.0; // should never get here
1742 }
1743
HasNoAlignedBlocks(void) const1744 bool BlockMultipleAlignment::HasNoAlignedBlocks(void) const
1745 {
1746 return (blocks.size() == 0 || (blocks.size() == 1 && !blocks.front()->IsAligned()));
1747 }
1748
GetAlignmentIndex(unsigned int row,unsigned int seqIndex,eUnalignedJustification justification) const1749 int BlockMultipleAlignment::GetAlignmentIndex(unsigned int row, unsigned int seqIndex, eUnalignedJustification justification) const
1750 {
1751 if (row >= NRows() || seqIndex >= GetSequenceOfRow(row)->Length()) {
1752 ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - coordinate out of range");
1753 return -1;
1754 }
1755
1756 unsigned int alignmentIndex, blockColumn;
1757 const Block *block = NULL;
1758 const Block::Range *range;
1759
1760 for (alignmentIndex=0; alignmentIndex<totalWidth; ++alignmentIndex) {
1761
1762 // check each block to see if index is in range
1763 if (block != blockMap[alignmentIndex].block) {
1764 block = blockMap[alignmentIndex].block;
1765
1766 range = block->GetRangeOfRow(row);
1767 if ((int)seqIndex >= range->from && (int)seqIndex <= range->to) {
1768
1769 // override requested justification for end blocks
1770 if (block == blocks.back()) // also true if there's a single aligned block
1771 justification = eLeft;
1772 else if (block == blocks.front())
1773 justification = eRight;
1774
1775 // search block columns to find index (inefficient, but avoids having to write
1776 // inverse functions of Block::GetIndexAt()
1777 for (blockColumn=0; blockColumn<block->width; ++blockColumn) {
1778 if ((int)seqIndex == block->GetIndexAt(blockColumn, row, justification))
1779 return alignmentIndex + blockColumn;
1780 }
1781 ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - can't find index in block");
1782 return -1;
1783 }
1784 }
1785 }
1786
1787 // should never get here...
1788 ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - confused");
1789 return -1;
1790 }
1791
1792 // returns the loop length of the unaligned region at the given row and alignment index; -1 if that position is aligned
GetLoopLength(unsigned int row,unsigned int alignmentIndex)1793 int BlockMultipleAlignment::GetLoopLength(unsigned int row, unsigned int alignmentIndex)
1794 {
1795 UnalignedBlock *block = dynamic_cast<UnalignedBlock*>(blockMap[alignmentIndex].block);
1796 if (!block)
1797 return -1;
1798
1799 const Block::Range *range = block->GetRangeOfRow(row);
1800 return (range->to - range->from + 1);
1801 }
1802
1803
1804 ///// UngappedAlignedBlock methods /////
1805
GetCharacterAt(unsigned int blockColumn,unsigned int row) const1806 char UngappedAlignedBlock::GetCharacterAt(unsigned int blockColumn, unsigned int row) const
1807 {
1808 return (*(parentAlignment->GetSequences()))[row]->sequenceString[GetIndexAt(blockColumn, row)];
1809 }
1810
Clone(const BlockMultipleAlignment * newMultiple) const1811 Block * UngappedAlignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
1812 {
1813 UngappedAlignedBlock *copy = new UngappedAlignedBlock(newMultiple);
1814 const Block::Range *range;
1815 for (unsigned int row=0; row<NSequences(); ++row) {
1816 range = GetRangeOfRow(row);
1817 copy->SetRangeOfRow(row, range->from, range->to);
1818 }
1819 copy->width = width;
1820 return copy;
1821 }
1822
DeleteRow(unsigned int row)1823 void UngappedAlignedBlock::DeleteRow(unsigned int row)
1824 {
1825 RangeList::iterator r = ranges.begin();
1826 for (unsigned int i=0; i<row; ++i) ++r;
1827 ranges.erase(r);
1828 }
1829
DeleteRows(vector<bool> & removeRows,unsigned int nToRemove)1830 void UngappedAlignedBlock::DeleteRows(vector < bool >& removeRows, unsigned int nToRemove)
1831 {
1832 VectorRemoveElements(ranges, removeRows, nToRemove);
1833 }
1834
1835
1836 ///// UnalignedBlock methods /////
1837
GetIndexAt(unsigned int blockColumn,unsigned int row,BlockMultipleAlignment::eUnalignedJustification justification) const1838 int UnalignedBlock::GetIndexAt(unsigned int blockColumn, unsigned int row,
1839 BlockMultipleAlignment::eUnalignedJustification justification) const
1840 {
1841 const Block::Range *range = GetRangeOfRow(row);
1842 int seqIndex = -1;
1843 unsigned int rangeWidth, rangeMiddle, extraSpace;
1844
1845 switch (justification) {
1846 case BlockMultipleAlignment::eLeft:
1847 seqIndex = range->from + blockColumn;
1848 break;
1849 case BlockMultipleAlignment::eRight:
1850 seqIndex = range->to - width + blockColumn + 1;
1851 break;
1852 case BlockMultipleAlignment::eCenter:
1853 rangeWidth = (range->to - range->from + 1);
1854 extraSpace = (width - rangeWidth) / 2;
1855 if (blockColumn < extraSpace || blockColumn >= extraSpace + rangeWidth)
1856 seqIndex = -1;
1857 else
1858 seqIndex = range->from + blockColumn - extraSpace;
1859 break;
1860 case BlockMultipleAlignment::eSplit:
1861 rangeWidth = (range->to - range->from + 1);
1862 rangeMiddle = (rangeWidth / 2) + (rangeWidth % 2);
1863 extraSpace = width - rangeWidth;
1864 if (blockColumn < rangeMiddle)
1865 seqIndex = range->from + blockColumn;
1866 else if (blockColumn >= extraSpace + rangeMiddle)
1867 seqIndex = range->to - width + blockColumn + 1;
1868 else
1869 seqIndex = -1;
1870 break;
1871 }
1872 if (seqIndex < range->from || seqIndex > range->to) seqIndex = -1;
1873
1874 return seqIndex;
1875 }
1876
Resize(void)1877 void UnalignedBlock::Resize(void)
1878 {
1879 width = 0;
1880 for (unsigned int i=0; i<NSequences(); ++i) {
1881 unsigned int blockWidth = ranges[i].to - ranges[i].from + 1;
1882 if (blockWidth > width) width = blockWidth;
1883 }
1884 }
1885
DeleteRow(unsigned int row)1886 void UnalignedBlock::DeleteRow(unsigned int row)
1887 {
1888 RangeList::iterator r = ranges.begin();
1889 for (unsigned int i=0; i<row; ++i) ++r;
1890 ranges.erase(r);
1891 Resize();
1892 }
1893
DeleteRows(vector<bool> & removeRows,unsigned int nToRemove)1894 void UnalignedBlock::DeleteRows(vector < bool >& removeRows, unsigned int nToRemove)
1895 {
1896 VectorRemoveElements(ranges, removeRows, nToRemove);
1897 Resize();
1898 }
1899
Clone(const BlockMultipleAlignment * newMultiple) const1900 Block * UnalignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
1901 {
1902 UnalignedBlock *copy = new UnalignedBlock(newMultiple);
1903 const Block::Range *range;
1904 for (unsigned int row=0; row<NSequences(); ++row) {
1905 range = GetRangeOfRow(row);
1906 copy->SetRangeOfRow(row, range->from, range->to);
1907 }
1908 copy->width = width;
1909 return copy;
1910 }
1911
MinResidues(void) const1912 unsigned int UnalignedBlock::MinResidues(void) const
1913 {
1914 int min = 0, m;
1915 const Block::Range *range;
1916 for (unsigned int row=0; row<NSequences(); ++row) {
1917 range = GetRangeOfRow(row);
1918 m = range->to - range->from + 1;
1919 if (row == 0 || m < min)
1920 min = m;
1921 }
1922 return min;
1923 }
1924
1925 END_SCOPE(Cn3D)
1926