1 /* $Id: alignment_manager.cpp 411368 2013-08-28 11:25:58Z thiessen $
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 manipulate alignments
30 *
31 * ===========================================================================
32 */
33
34 #include <ncbi_pch.hpp>
35 #include <corelib/ncbistd.hpp>
36
37 #include <objects/seqalign/Dense_diag.hpp>
38 #include <objects/seqalign/Dense_seg.hpp>
39 #include <objects/seqset/Bioseq_set.hpp>
40 #include <objects/scoremat/scoremat__.hpp>
41 #include <objects/cdd/Cdd_id.hpp>
42 #include <objects/cdd/Cdd_id_set.hpp>
43
44 #include "remove_header_conflicts.hpp"
45
46 #include "alignment_manager.hpp"
47 #include "sequence_set.hpp"
48 #include "alignment_set.hpp"
49 #include "block_multiple_alignment.hpp"
50 #include "messenger.hpp"
51 #include "structure_set.hpp"
52 #include "sequence_viewer.hpp"
53 #include "molecule.hpp"
54 #include "show_hide_manager.hpp"
55 #include "cn3d_threader.hpp"
56 #include "update_viewer.hpp"
57 #include "sequence_display.hpp"
58 #include "cn3d_tools.hpp"
59 #include "molecule_identifier.hpp"
60 #include "cn3d_blast.hpp"
61 #include "style_manager.hpp"
62 #include "cn3d_ba_interface.hpp"
63 #include "cn3d_pssm.hpp"
64 #include "progress_meter.hpp"
65
66 #include <algo/structure/bma_refine/Interface.hpp>
67 #include <algo/structure/bma_refine/InterfaceGUI.hpp>
68
69 USING_NCBI_SCOPE;
70 USING_SCOPE(objects);
71
72
BEGIN_SCOPE(Cn3D)73 BEGIN_SCOPE(Cn3D)
74
75 ///// AlignmentManager methods /////
76
77 void AlignmentManager::Init(void)
78 {
79 sequenceViewer = new SequenceViewer(this);
80 GlobalMessenger()->AddSequenceViewer(sequenceViewer);
81
82 updateViewer = new UpdateViewer(this);
83 GlobalMessenger()->AddSequenceViewer(updateViewer);
84
85 threader = new Threader(this);
86 blaster = new BLASTer();
87 blockAligner = new BlockAligner();
88 // bmaRefiner = new BMARefiner();
89
90 originalMultiple = NULL;
91 }
92
AlignmentManager(const SequenceSet * sSet,const AlignmentSet * aSet)93 AlignmentManager::AlignmentManager(const SequenceSet *sSet, const AlignmentSet *aSet)
94 {
95 Init();
96 NewAlignments(sSet, aSet);
97 }
98
AlignmentManager(const SequenceSet * sSet,const AlignmentSet * aSet,const UpdateAlignList & updates)99 AlignmentManager::AlignmentManager(const SequenceSet *sSet,
100 const AlignmentSet *aSet, const UpdateAlignList& updates)
101 {
102 Init();
103 NewAlignments(sSet, aSet);
104
105 // create BlockMultipleAlignments from updates; add to update viewer
106 PairwiseAlignmentList pairwise(1);
107 UpdateViewer::AlignmentList updateAlignments;
108 UpdateAlignList::const_iterator u, ue = updates.end();
109 for (u=updates.begin(); u!=ue; ++u) {
110 if (u->GetObject().IsSetSeqannot() && u->GetObject().GetSeqannot().GetData().IsAlign()) {
111 CSeq_annot::C_Data::TAlign::const_iterator
112 s, se = u->GetObject().GetSeqannot().GetData().GetAlign().end();
113 for (s=u->GetObject().GetSeqannot().GetData().GetAlign().begin(); s!=se; ++s) {
114
115 // determine master sequence
116 const Sequence *master = NULL;
117 if (aSet) {
118 master = aSet->master;
119 } else if (updateAlignments.size() > 0) {
120 master = updateAlignments.front()->GetMaster();
121 } else {
122 SequenceSet::SequenceList::const_iterator q, qe = sSet->sequences.end();
123 for (q=sSet->sequences.begin(); q!=qe; ++q) {
124 if ((*q)->identifier->MatchesSeqId(
125 (*s)->GetSegs().IsDendiag() ?
126 (*s)->GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
127 (*s)->GetSegs().GetDenseg().GetIds().front().GetObject()
128 )) {
129 master = *q;
130 break;
131 }
132 }
133 }
134 if (!master) {
135 ERRORMSG("AlignmentManager::AlignmentManager() - "
136 << "can't determine master sequence for updates");
137 return;
138 }
139
140 const MasterDependentAlignment *alignment =
141 new MasterDependentAlignment(NULL, master, s->GetObject());
142 pairwise.front() = alignment;
143 BlockMultipleAlignment *multiple = CreateMultipleFromPairwiseWithIBM(pairwise);
144 multiple->updateOrigin = *u; // to keep track of which Update-align this came from
145 updateAlignments.push_back(multiple);
146 delete alignment;
147 }
148 }
149 }
150 updateViewer->AddAlignments(updateAlignments);
151
152 // set this set of updates as the initial state of the editor's undo stack
153 updateViewer->SetInitialState();
154 }
155
~AlignmentManager(void)156 AlignmentManager::~AlignmentManager(void)
157 {
158 GlobalMessenger()->RemoveSequenceViewer(sequenceViewer);
159 GlobalMessenger()->RemoveSequenceViewer(updateViewer);
160 delete sequenceViewer;
161 delete updateViewer;
162 delete threader;
163 delete blaster;
164 delete blockAligner;
165 // delete bmaRefiner;
166 if (originalMultiple) delete originalMultiple;
167 }
168
NewAlignments(const SequenceSet * sSet,const AlignmentSet * aSet)169 void AlignmentManager::NewAlignments(const SequenceSet *sSet, const AlignmentSet *aSet)
170 {
171 sequenceSet = sSet;
172 alignmentSet = aSet;
173
174 if (!alignmentSet) {
175 sequenceViewer->DisplaySequences(&(sequenceSet->sequences));
176 return;
177 }
178
179 // all dependents start out visible
180 dependentsVisible.resize(alignmentSet->alignments.size());
181 for (unsigned int i=0; i<dependentsVisible.size(); ++i)
182 dependentsVisible[i] = true;
183
184 NewMultipleWithRows(dependentsVisible);
185 originalMultiple = GetCurrentMultipleAlignment()->Clone();
186 originalRowOrder.resize(originalMultiple->NRows());
187 for (unsigned int r=0; r<originalMultiple->NRows(); ++r)
188 originalRowOrder[r] = r;
189 }
190
SavePairwiseFromMultiple(const BlockMultipleAlignment * multiple,const vector<unsigned int> & rowOrder)191 void AlignmentManager::SavePairwiseFromMultiple(const BlockMultipleAlignment *multiple,
192 const vector < unsigned int >& rowOrder)
193 {
194 // create new AlignmentSet based on this multiple alignment, feed back into StructureSet
195 AlignmentSet *newAlignmentSet =
196 AlignmentSet::CreateFromMultiple(multiple->GetMaster()->parentSet, multiple, rowOrder);
197 if (newAlignmentSet) {
198 multiple->GetMaster()->parentSet->ReplaceAlignmentSet(newAlignmentSet);
199 alignmentSet = newAlignmentSet;
200 } else {
201 ERRORMSG("Couldn't create pairwise alignments from the current multiple!\n"
202 << "Alignment data in output file will be left unchanged.");
203 return;
204 }
205
206 // see whether PSSM and/or row order have changed
207 if (!originalMultiple) {
208 originalMultiple = multiple->Clone();
209 originalRowOrder = rowOrder;
210 alignmentSet->parentSet->SetDataChanged(StructureSet::ePSSMData);
211 alignmentSet->parentSet->SetDataChanged(StructureSet::eRowOrderData);
212 } else {
213
214 // check for row order change
215 TRACEMSG("checking for row order changes..." << originalMultiple << ' ' << multiple);
216 if (originalMultiple->NRows() != multiple->NRows()) {
217 TRACEMSG("row order changed");
218 alignmentSet->parentSet->SetDataChanged(StructureSet::eRowOrderData);
219 } else {
220 for (unsigned int row=0; row<originalMultiple->NRows(); ++row) {
221 if (originalMultiple->GetSequenceOfRow(originalRowOrder[row]) !=
222 multiple->GetSequenceOfRow(rowOrder[row]) ||
223 originalRowOrder[row] != rowOrder[row])
224 {
225 TRACEMSG("row order changed");
226 alignmentSet->parentSet->SetDataChanged(StructureSet::eRowOrderData);
227 break;
228 }
229 }
230 }
231
232 // check for PSSM change
233 if ((multiple->HasNoAlignedBlocks() && !originalMultiple->HasNoAlignedBlocks()) ||
234 (originalMultiple->HasNoAlignedBlocks() && !multiple->HasNoAlignedBlocks())) {
235 TRACEMSG("PSSM changed - zero blocks before or after");
236 alignmentSet->parentSet->SetDataChanged(StructureSet::ePSSMData);
237 } else if (!multiple->HasNoAlignedBlocks() && !originalMultiple->HasNoAlignedBlocks()) {
238 const CPssmWithParameters
239 &originalPSSM = originalMultiple->GetPSSM().GetPSSM(),
240 ¤tPSSM = multiple->GetPSSM().GetPSSM();
241 TRACEMSG("checking for PSSM changes... ");
242 if (originalPSSM.GetPssm().GetNumRows() != currentPSSM.GetPssm().GetNumRows() ||
243 originalPSSM.GetPssm().GetNumColumns() != currentPSSM.GetPssm().GetNumColumns() ||
244 originalPSSM.GetPssm().GetByRow() != currentPSSM.GetPssm().GetByRow() ||
245 !originalPSSM.GetPssm().IsSetFinalData() || !currentPSSM.GetPssm().IsSetFinalData() ||
246 originalPSSM.GetPssm().GetFinalData().GetLambda() != currentPSSM.GetPssm().GetFinalData().GetLambda() ||
247 originalPSSM.GetPssm().GetFinalData().GetKappa() != currentPSSM.GetPssm().GetFinalData().GetKappa() ||
248 originalPSSM.GetPssm().GetFinalData().GetH() != currentPSSM.GetPssm().GetFinalData().GetH() ||
249 originalPSSM.GetPssm().GetFinalData().GetScalingFactor() != currentPSSM.GetPssm().GetFinalData().GetScalingFactor()
250 ) {
251 TRACEMSG("PSSM changed");
252 alignmentSet->parentSet->SetDataChanged(StructureSet::ePSSMData);
253 } else {
254 CPssmFinalData::TScores::const_iterator
255 o = originalPSSM.GetPssm().GetFinalData().GetScores().begin(),
256 oe = originalPSSM.GetPssm().GetFinalData().GetScores().end(),
257 c = currentPSSM.GetPssm().GetFinalData().GetScores().begin();
258 for (; o!=oe; ++o, ++c) {
259 if ((*o) != (*c)) {
260 TRACEMSG("PSSM changed");
261 alignmentSet->parentSet->SetDataChanged(StructureSet::ePSSMData);
262 break;
263 }
264 }
265 }
266 }
267
268 // keep for comparison on next save
269 delete originalMultiple;
270 originalMultiple = multiple->Clone();
271 originalMultiple->RemovePSSM();
272 originalRowOrder = rowOrder;
273 }
274 }
275
GetCurrentMultipleAlignment(void) const276 const BlockMultipleAlignment * AlignmentManager::GetCurrentMultipleAlignment(void) const
277 {
278 const ViewerBase::AlignmentList& currentAlignments = sequenceViewer->GetCurrentAlignments();
279 return ((currentAlignments.size() > 0) ? currentAlignments.front() : NULL);
280 }
281
AlignedToAllDependents(int masterResidue,const AlignmentManager::PairwiseAlignmentList & alignments)282 static bool AlignedToAllDependents(int masterResidue,
283 const AlignmentManager::PairwiseAlignmentList& alignments)
284 {
285 AlignmentManager::PairwiseAlignmentList::const_iterator a, ae = alignments.end();
286 for (a=alignments.begin(); a!=ae; ++a) {
287 if ((*a)->masterToDependent[masterResidue] == -1) return false;
288 }
289 return true;
290 }
291
NoDependentInsertionsBetween(int masterFrom,int masterTo,const AlignmentManager::PairwiseAlignmentList & alignments)292 static bool NoDependentInsertionsBetween(int masterFrom, int masterTo,
293 const AlignmentManager::PairwiseAlignmentList& alignments)
294 {
295 AlignmentManager::PairwiseAlignmentList::const_iterator a, ae = alignments.end();
296 for (a=alignments.begin(); a!=ae; ++a) {
297 if (((*a)->masterToDependent[masterTo] - (*a)->masterToDependent[masterFrom]) != (masterTo - masterFrom))
298 return false;
299 }
300 return true;
301 }
302
NoBlockBoundariesBetween(int masterFrom,int masterTo,const AlignmentManager::PairwiseAlignmentList & alignments)303 static bool NoBlockBoundariesBetween(int masterFrom, int masterTo,
304 const AlignmentManager::PairwiseAlignmentList& alignments)
305 {
306 AlignmentManager::PairwiseAlignmentList::const_iterator a, ae = alignments.end();
307 for (a=alignments.begin(); a!=ae; ++a) {
308 if ((*a)->blockStructure[masterTo] != (*a)->blockStructure[masterFrom])
309 return false;
310 }
311 return true;
312 }
313
314 BlockMultipleAlignment *
CreateMultipleFromPairwiseWithIBM(const PairwiseAlignmentList & alignments)315 AlignmentManager::CreateMultipleFromPairwiseWithIBM(const PairwiseAlignmentList& alignments)
316 {
317 PairwiseAlignmentList::const_iterator a, ae = alignments.end();
318
319 // create sequence list; fill with sequences of master + dependents
320 BlockMultipleAlignment::SequenceList
321 *sequenceList = new BlockMultipleAlignment::SequenceList(alignments.size() + 1);
322 BlockMultipleAlignment::SequenceList::iterator s = sequenceList->begin();
323 *(s++) = alignments.front()->master;
324 for (a=alignments.begin(); a!=ae; ++a) {
325 *(s++) = (*a)->dependent;
326 if ((*a)->master != sequenceList->front()) {
327 ERRORMSG("AlignmentManager::CreateMultipleFromPairwiseWithIBM() -\n"
328 << "all pairwise alignments must have the same master sequence");
329 return NULL;
330 }
331 }
332 BlockMultipleAlignment *multipleAlignment = new BlockMultipleAlignment(sequenceList, this);
333
334 // each block is a continuous region on the master, over which each master
335 // residue is aligned to a residue of each dependent, and where there are no
336 // insertions relative to the master in any of the dependents
337 int masterFrom = 0, masterTo;
338 unsigned int row;
339 UngappedAlignedBlock *newBlock;
340
341 while (masterFrom < (int)multipleAlignment->GetMaster()->Length()) {
342
343 // look for first all-aligned residue
344 if (!AlignedToAllDependents(masterFrom, alignments)) {
345 ++masterFrom;
346 continue;
347 }
348
349 // find all next continuous all-aligned residues, but checking for
350 // block boundaries from the original master-dependent pairs, so that
351 // blocks don't get merged
352 for (masterTo=masterFrom+1;
353 masterTo < (int)multipleAlignment->GetMaster()->Length() &&
354 AlignedToAllDependents(masterTo, alignments) &&
355 NoDependentInsertionsBetween(masterFrom, masterTo, alignments) &&
356 NoBlockBoundariesBetween(masterFrom, masterTo, alignments);
357 ++masterTo) ;
358 --masterTo; // after loop, masterTo = first residue past block
359
360 // create new block with ranges from master and all dependents
361 newBlock = new UngappedAlignedBlock(multipleAlignment);
362 newBlock->SetRangeOfRow(0, masterFrom, masterTo);
363 newBlock->width = masterTo - masterFrom + 1;
364
365 //TESTMSG("masterFrom " << masterFrom+1 << ", masterTo " << masterTo+1);
366 for (a=alignments.begin(), row=1; a!=ae; ++a, ++row) {
367 newBlock->SetRangeOfRow(row,
368 (*a)->masterToDependent[masterFrom],
369 (*a)->masterToDependent[masterTo]);
370 //TESTMSG("dependent->from " << b->from+1 << ", dependent->to " << b->to+1);
371 }
372
373 // copy new block into alignment
374 multipleAlignment->AddAlignedBlockAtEnd(newBlock);
375
376 // start looking for next block
377 masterFrom = masterTo + 1;
378 }
379
380 if (!multipleAlignment->AddUnalignedBlocks() ||
381 !multipleAlignment->UpdateBlockMapAndColors()) {
382 ERRORMSG("AlignmentManager::CreateMultipleFromPairwiseWithIBM() - error finalizing alignment");
383 return NULL;
384 }
385
386 return multipleAlignment;
387 }
388
GetAlignedResidueIndexes(BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator & b,BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator & be,int row,int * seqIndexes,bool countHighlights=false,const BlockMultipleAlignment * multiple=NULL)389 static int GetAlignedResidueIndexes(
390 BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator& b,
391 BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator& be,
392 int row, int *seqIndexes, bool countHighlights = false, const BlockMultipleAlignment *multiple = NULL)
393 {
394 unsigned int i = 0, c, highlighted = 0;
395 const Block::Range *range;
396 for (; b!=be; ++b) {
397 range = (*b)->GetRangeOfRow(row);
398 for (c=0; c<(*b)->width; ++c, ++i) {
399 seqIndexes[i] = range->from + c;
400 if (countHighlights) {
401 if (GlobalMessenger()->IsHighlighted(multiple->GetSequenceOfRow(row), seqIndexes[i]))
402 ++highlighted;
403 else
404 seqIndexes[i] = -1;
405 }
406 }
407 }
408 return highlighted;
409 }
410
RealignAllDependentStructures(bool highlightedOnly) const411 void AlignmentManager::RealignAllDependentStructures(bool highlightedOnly) const
412 {
413 const BlockMultipleAlignment *multiple = GetCurrentMultipleAlignment();
414 if (!multiple) return;
415 BlockMultipleAlignment::UngappedAlignedBlockList blocks;
416 multiple->GetUngappedAlignedBlocks(&blocks);
417 BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
418 int nResidues = 0;
419 for (b=blocks.begin(); b!=be; ++b)
420 nResidues += (*b)->width;
421 if (nResidues == 0) {
422 WARNINGMSG("Can't realign dependents with no aligned residues!");
423 return;
424 }
425
426 const Sequence *masterSeq = multiple->GetSequenceOfRow(0), *dependentSeq;
427 const Molecule *masterMol, *dependentMol;
428 if (!masterSeq || !(masterMol = masterSeq->molecule)) {
429 WARNINGMSG("Can't realign dependents to non-structured master!");
430 return;
431 }
432
433 int *masterSeqIndexes = new int[nResidues], *dependentSeqIndexes = new int[nResidues];
434 b = blocks.begin();
435 int nHighlightedAligned = GetAlignedResidueIndexes(b, be, 0, masterSeqIndexes, highlightedOnly, multiple);
436 if ((highlightedOnly ? nHighlightedAligned : nResidues) < 3) {
437 WARNINGMSG("Can't realign dependents using < 3 residues!");
438 delete[] masterSeqIndexes;
439 delete[] dependentSeqIndexes;
440 return;
441 }
442
443 double *weights = new double[nResidues];
444 const StructureObject *dependentObj;
445
446 typedef const Vector * CVP;
447 CVP *masterCoords = new CVP[nResidues], *dependentCoords = new CVP[nResidues];
448 if (!masterMol->GetAlphaCoords(nResidues, masterSeqIndexes, masterCoords)) {
449 WARNINGMSG("Can't get master alpha coords");
450 } else if (masterSeq->GetOrSetMMDBLink() == MoleculeIdentifier::VALUE_NOT_SET) {
451 WARNINGMSG("Don't know master MMDB ID");
452 } else {
453
454 masterMol->parentSet->InitStructureAlignments(masterSeq->identifier->mmdbID);
455
456 unsigned int nStructureAlignments = 0;
457 for (unsigned int i=1; i<multiple->NRows(); ++i) {
458 dependentSeq = multiple->GetSequenceOfRow(i);
459 if (!dependentSeq || !(dependentMol = dependentSeq->molecule)) continue;
460
461 b = blocks.begin();
462 GetAlignedResidueIndexes(b, be, i, dependentSeqIndexes);
463 if (dependentMol->GetAlphaCoords(nResidues, dependentSeqIndexes, dependentCoords) < 3) {
464 ERRORMSG("can't realign dependent " << dependentSeq->identifier->pdbID << ", not enough coordinates in aligned region");
465 continue;
466 }
467
468 if (!dependentMol->GetParentOfType(&dependentObj)) continue;
469
470 // if any Vector* is NULL, make sure that weight is 0 so the pointer won't be accessed
471 int nWeighted = 0;
472 for (int j=0; j<nResidues; ++j) {
473 if (!masterCoords[j] || !dependentCoords[j] || (highlightedOnly && masterSeqIndexes[j] < 0)) {
474 weights[j] = 0.0;
475 } else {
476 weights[j] = 1.0; // for now, just use flat weighting
477 ++nWeighted;
478 }
479 }
480 if (nWeighted < 3) {
481 WARNINGMSG("Can't realign dependent #" << (i+1) << " using < 3 residues!");
482 continue;
483 }
484
485 INFOMSG("realigning dependent " << dependentSeq->identifier->pdbID << " against master " << masterSeq->identifier->pdbID
486 << " using coordinates of " << nWeighted << " residues");
487 (const_cast<StructureObject*>(dependentObj))->RealignStructure(nResidues, masterCoords, dependentCoords, weights, i);
488 ++nStructureAlignments;
489 }
490
491 // if no structure alignments, remove the list entirely
492 if (nStructureAlignments == 0) masterMol->parentSet->RemoveStructureAlignments();
493 }
494
495 delete[] masterSeqIndexes;
496 delete[] dependentSeqIndexes;
497 delete[] masterCoords;
498 delete[] dependentCoords;
499 delete[] weights;
500 return;
501 }
502
GetAlignmentSetDependentSequences(vector<const Sequence * > * sequences) const503 void AlignmentManager::GetAlignmentSetDependentSequences(vector < const Sequence * > *sequences) const
504 {
505 sequences->resize(alignmentSet->alignments.size());
506
507 AlignmentSet::AlignmentList::const_iterator a, ae = alignmentSet->alignments.end();
508 int i = 0;
509 for (a=alignmentSet->alignments.begin(); a!=ae; ++a, ++i) {
510 (*sequences)[i] = (*a)->dependent;
511 }
512 }
513
GetAlignmentSetDependentVisibilities(vector<bool> * visibilities) const514 void AlignmentManager::GetAlignmentSetDependentVisibilities(vector < bool > *visibilities) const
515 {
516 if (dependentsVisible.size() != alignmentSet->alignments.size()) // can happen if row is added/deleted
517 dependentsVisible.resize(alignmentSet->alignments.size(), true);
518
519 // copy visibility list
520 *visibilities = dependentsVisible;
521 }
522
ShowHideCallbackFunction(const vector<bool> & itemsEnabled)523 void AlignmentManager::ShowHideCallbackFunction(const vector < bool >& itemsEnabled)
524 {
525 if (itemsEnabled.size() != dependentsVisible.size() ||
526 itemsEnabled.size() != alignmentSet->alignments.size()) {
527 ERRORMSG("AlignmentManager::ShowHideCallbackFunction() - wrong size list");
528 return;
529 }
530
531 dependentsVisible = itemsEnabled;
532 NewMultipleWithRows(dependentsVisible);
533
534 AlignmentSet::AlignmentList::const_iterator
535 a = alignmentSet->alignments.begin(), ae = alignmentSet->alignments.end();
536 const StructureObject *object;
537
538 if ((*a)->master->molecule) {
539 // Show() redraws whole StructureObject only if necessary
540 if ((*a)->master->molecule->GetParentOfType(&object))
541 object->parentSet->showHideManager->Show(object, true);
542 // always redraw aligned molecule, in case alignment colors change
543 GlobalMessenger()->PostRedrawMolecule((*a)->master->molecule);
544 }
545 for (int i=0; a!=ae; ++a, ++i) {
546 if ((*a)->dependent->molecule) {
547 if ((*a)->dependent->molecule->GetParentOfType(&object))
548 object->parentSet->showHideManager->Show(object, dependentsVisible[i]);
549 GlobalMessenger()->PostRedrawMolecule((*a)->dependent->molecule);
550 }
551 }
552
553 // do necessary redraws + show/hides: sequences + chains in the alignment
554 sequenceViewer->Refresh();
555 GlobalMessenger()->PostRedrawAllSequenceViewers();
556 GlobalMessenger()->UnPostRedrawSequenceViewer(sequenceViewer); // Refresh() does this already
557 }
558
NewMultipleWithRows(const vector<bool> & visibilities)559 void AlignmentManager::NewMultipleWithRows(const vector < bool >& visibilities)
560 {
561 if (visibilities.size() != alignmentSet->alignments.size()) {
562 ERRORMSG("AlignmentManager::NewMultipleWithRows() - wrong size visibility vector");
563 return;
564 }
565
566 // make a multiple from all visible rows
567 PairwiseAlignmentList alignments;
568 AlignmentSet::AlignmentList::const_iterator a, ae=alignmentSet->alignments.end();
569 int i = 0;
570 for (a=alignmentSet->alignments.begin(); a!=ae; ++a, ++i) {
571 if (visibilities[i])
572 alignments.push_back(*a);
573 }
574
575 // sequenceViewer will own the resulting alignment
576 sequenceViewer->DisplayAlignment(CreateMultipleFromPairwiseWithIBM(alignments));
577 }
578
IsAligned(const Sequence * sequence,unsigned int seqIndex) const579 bool AlignmentManager::IsAligned(const Sequence *sequence, unsigned int seqIndex) const
580 {
581 if (!sequence) return false;
582 const BlockMultipleAlignment *currentAlignment = GetCurrentMultipleAlignment();
583 if (currentAlignment)
584 return currentAlignment->IsAligned(sequence, seqIndex);
585 else
586 return false;
587 }
588
IsInAlignment(const Sequence * sequence) const589 bool AlignmentManager::IsInAlignment(const Sequence *sequence) const
590 {
591 if (!sequence) return false;
592 const BlockMultipleAlignment *currentAlignment = GetCurrentMultipleAlignment();
593 if (currentAlignment) {
594 for (unsigned int i=0; i<currentAlignment->GetSequences()->size(); ++i) {
595 if ((*(currentAlignment->GetSequences()))[i] == sequence)
596 return true;
597 }
598 }
599 return false;
600 }
601
GetAlignmentColor(const Sequence * sequence,unsigned int seqIndex,StyleSettings::eColorScheme colorScheme) const602 const Vector * AlignmentManager::GetAlignmentColor(const Sequence *sequence, unsigned int seqIndex,
603 StyleSettings::eColorScheme colorScheme) const
604 {
605 const BlockMultipleAlignment *currentAlignment = GetCurrentMultipleAlignment();
606 if (currentAlignment)
607 return currentAlignment->GetAlignmentColor(sequence, seqIndex, colorScheme);
608 else
609 return NULL;
610 }
611
ShowSequenceViewer(bool showNow) const612 void AlignmentManager::ShowSequenceViewer(bool showNow) const
613 {
614 sequenceViewer->CreateSequenceWindow(showNow);
615 }
616
ShowUpdateWindow(void) const617 void AlignmentManager::ShowUpdateWindow(void) const
618 {
619 updateViewer->CreateUpdateWindow();
620 }
621
RealignDependentSequences(BlockMultipleAlignment * multiple,const vector<unsigned int> & dependentsToRealign)622 void AlignmentManager::RealignDependentSequences(
623 BlockMultipleAlignment *multiple, const vector < unsigned int >& dependentsToRealign)
624 {
625 if (!multiple || sequenceViewer->GetCurrentAlignments().size() == 0 ||
626 multiple != sequenceViewer->GetCurrentAlignments().front()) {
627 ERRORMSG("AlignmentManager::RealignDependentSequences() - wrong multiple alignment");
628 return;
629 }
630 if (dependentsToRealign.size() == 0) return;
631
632 // create alignments for each master/dependent pair, then update displays
633 UpdateViewer::AlignmentList alignments;
634 TRACEMSG("extracting rows");
635 if (multiple->ExtractRows(dependentsToRealign, &alignments)) {
636 TRACEMSG("recreating display");
637 sequenceViewer->GetCurrentDisplay()->RowsRemoved(dependentsToRealign, multiple);
638 TRACEMSG("adding to update window");
639 SetDiagPostLevel(eDiag_Warning); // otherwise, info messages take a long time if lots of rows
640 updateViewer->AddAlignments(alignments);
641 SetDiagPostLevel(eDiag_Info);
642 TRACEMSG("done");
643 updateViewer->CreateUpdateWindow();
644 }
645 }
646
ThreadUpdate(const ThreaderOptions & options,BlockMultipleAlignment * single)647 void AlignmentManager::ThreadUpdate(const ThreaderOptions& options, BlockMultipleAlignment *single)
648 {
649 const ViewerBase::AlignmentList& currentAlignments = sequenceViewer->GetCurrentAlignments();
650 if (currentAlignments.size() == 0) return;
651
652 // run the threader on the given alignment
653 UpdateViewer::AlignmentList singleList, replacedList;
654 Threader::AlignmentList newAlignments;
655 unsigned int nRowsAddedToMultiple;
656 bool foundSingle = false; // sanity check
657
658 singleList.push_back(single);
659 if (threader->Realign(
660 options, currentAlignments.front(), &singleList, &newAlignments, &nRowsAddedToMultiple, sequenceViewer)) {
661
662 // replace threaded alignment with new one(s) (or leftover where threader/merge failed)
663 UpdateViewer::AlignmentList::const_iterator a, ae = updateViewer->GetCurrentAlignments().end();
664 for (a=updateViewer->GetCurrentAlignments().begin(); a!=ae; ++a) {
665 if (*a == single) {
666 Threader::AlignmentList::const_iterator n, ne = newAlignments.end();
667 for (n=newAlignments.begin(); n!=ne; ++n)
668 replacedList.push_back(*n);
669 foundSingle = true;
670 } else
671 replacedList.push_back((*a)->Clone());
672 }
673 if (!foundSingle) ERRORMSG(
674 "AlignmentManager::ThreadUpdate() - threaded alignment not found in update viewer!");
675 updateViewer->ReplaceAlignments(replacedList);
676
677 // tell the sequenceViewer that rows have been merged into the multiple
678 if (nRowsAddedToMultiple > 0)
679 sequenceViewer->GetCurrentDisplay()->RowsAdded(nRowsAddedToMultiple, currentAlignments.front());
680 }
681 }
682
ThreadAllUpdates(const ThreaderOptions & options)683 void AlignmentManager::ThreadAllUpdates(const ThreaderOptions& options)
684 {
685 const ViewerBase::AlignmentList& currentAlignments = sequenceViewer->GetCurrentAlignments();
686 if (currentAlignments.size() == 0) return;
687
688 // run the threader on update pairwise alignments
689 Threader::AlignmentList newAlignments;
690 unsigned int nRowsAddedToMultiple;
691 if (threader->Realign(
692 options, currentAlignments.front(), &(updateViewer->GetCurrentAlignments()),
693 &newAlignments, &nRowsAddedToMultiple, sequenceViewer)) {
694
695 // replace update alignments with new ones (or leftovers where threader/merge failed)
696 updateViewer->ReplaceAlignments(newAlignments);
697
698 // tell the sequenceViewer that rows have been merged into the multiple
699 if (nRowsAddedToMultiple > 0)
700 sequenceViewer->GetCurrentDisplay()->
701 RowsAdded(nRowsAddedToMultiple, currentAlignments.front());
702 }
703 }
704
MergeUpdates(const AlignmentManager::UpdateMap & updatesToMerge,bool mergeToNeighbor)705 void AlignmentManager::MergeUpdates(const AlignmentManager::UpdateMap& updatesToMerge, bool mergeToNeighbor)
706 {
707 if (updatesToMerge.size() == 0) return;
708 const ViewerBase::AlignmentList& currentUpdates = updateViewer->GetCurrentAlignments();
709 if (currentUpdates.size() == 0) return;
710
711 // transform this structure view into an alignment view, and turn on the editor.
712 ViewerBase::AlignmentList::const_iterator u, ue = currentUpdates.end();
713 const BlockMultipleAlignment *newMultiple = NULL;
714 if (sequenceViewer->GetCurrentAlignments().size() == 0) {
715
716 for (u=currentUpdates.begin(); u!=ue; ++u) { // find first update alignment
717 if (updatesToMerge.find(*u) != updatesToMerge.end()) {
718 newMultiple = *u;
719
720 // create new alignment, then call SavePairwiseFromMultiple to create
721 // an AlignmentSet and the initial ASN data
722 sequenceViewer->DisplayAlignment(newMultiple->Clone());
723 vector < unsigned int > rowOrder(newMultiple->NRows());
724 for (unsigned int i=0; i<newMultiple->NRows(); ++i) rowOrder[i] = i;
725 SavePairwiseFromMultiple(newMultiple, rowOrder);
726
727 // editor needs to be on if >1 update is to be merged in
728 sequenceViewer->TurnOnEditor();
729
730 // set default alignment-type style
731 newMultiple->GetMaster()->parentSet->styleManager->
732 SetGlobalRenderingStyle(StyleSettings::eTubeShortcut);
733 newMultiple->GetMaster()->parentSet->styleManager->
734 SetGlobalColorScheme(StyleSettings::eAlignedShortcut);
735 GlobalMessenger()->PostRedrawAllStructures();
736 break;
737 }
738 }
739 }
740
741 BlockMultipleAlignment *multiple =
742 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
743 sequenceViewer->GetCurrentAlignments().front() : NULL;
744 if (!multiple) {
745 ERRORMSG("Must have an alignment in the sequence viewer to merge with");
746 return;
747 }
748
749 // make sure the editor is on in the sequenceViewer
750 if (!sequenceViewer->EditorIsOn())
751 sequenceViewer->TurnOnEditor();
752
753 int nSuccessfulMerges = 0;
754 ViewerBase::AlignmentList updatesToKeep;
755 for (u=currentUpdates.begin(); u!=ue; ++u) {
756 if (*u == newMultiple) continue;
757 bool merged = false;
758 if (updatesToMerge.find(*u) != updatesToMerge.end()) {
759 merged = multiple->MergeAlignment(*u);
760 if (merged) {
761 nSuccessfulMerges += (*u)->NRows() - 1;
762 } else {
763 for (unsigned int i=0; i<(*u)->NRows(); ++i) {
764 string status = (*u)->GetRowStatusLine(i);
765 if (status.size() > 0)
766 status += "; merge failed!";
767 else
768 status = "Merge failed!";
769 (*u)->SetRowStatusLine(i, status);
770 }
771 }
772 }
773 if (!merged) {
774 BlockMultipleAlignment *keep = (*u)->Clone();
775 updatesToKeep.push_back(keep);
776 }
777 }
778
779 updateViewer->ReplaceAlignments(updatesToKeep);
780 if (nSuccessfulMerges > 0) {
781 int where = -1;
782
783 // if necessary, find nearest neighbor to merged sequence, and add new row after it
784 if (mergeToNeighbor && nSuccessfulMerges == 1) {
785 BlockMultipleAlignment::UngappedAlignedBlockList blocks;
786 multiple->GetUngappedAlignedBlocks(&blocks);
787 BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = blocks.end();
788 int rowScore, bestScore = 0;
789 unsigned int col, row, lastRow = multiple->NRows() - 1;
790 const Sequence *mergeSeq = multiple->GetSequenceOfRow(lastRow);
791 for (row=0; row<lastRow; ++row) {
792 const Sequence *otherSeq = multiple->GetSequenceOfRow(row);
793 rowScore = 0;
794 for (b=blocks.begin(); b!=be; ++b) {
795 for (col=0; col<(*b)->width; ++col) {
796 rowScore += GetBLOSUM62Score(
797 mergeSeq->sequenceString[(*b)->GetRangeOfRow(lastRow)->from + col],
798 otherSeq->sequenceString[(*b)->GetRangeOfRow(row)->from + col]);
799 }
800 }
801 TRACEMSG("Merge score with row " << (row + 1) << " (" << multiple->GetSequenceOfRow(row)->identifier->ToString() << ") is " << rowScore);
802 if (row == 0 || rowScore > bestScore) {
803 where = row;
804 bestScore = rowScore;
805 }
806 }
807 INFOMSG("Closest row is #" << (where+1) << ", "
808 << multiple->GetSequenceOfRow(where)->identifier->ToString());
809 }
810
811 sequenceViewer->GetCurrentDisplay()->RowsAdded(nSuccessfulMerges, multiple, where);
812 }
813
814 // add pending imported structures to asn data
815 updateViewer->SavePendingStructures();
816 }
817
CalculateRowScoresWithThreader(double weightPSSM)818 void AlignmentManager::CalculateRowScoresWithThreader(double weightPSSM)
819 {
820 threader->CalculateScores(GetCurrentMultipleAlignment(), weightPSSM);
821 }
822
NUpdates(void) const823 unsigned int AlignmentManager::NUpdates(void) const
824 {
825 return updateViewer->GetCurrentAlignments().size();
826 }
827
GetUpdateSequences(list<const Sequence * > * updateSequences) const828 void AlignmentManager::GetUpdateSequences(list < const Sequence * > *updateSequences) const
829 {
830 updateSequences->clear();
831 const ViewerBase::AlignmentList& currentUpdates = updateViewer->GetCurrentAlignments();
832 if (currentUpdates.size() == 0) return;
833 ViewerBase::AlignmentList::const_iterator u, ue = currentUpdates.end();
834 for (u=currentUpdates.begin(); u!=ue; ++u)
835 updateSequences->push_back((*u)->GetSequenceOfRow(1)); // assume update aln has just one dependent...
836 }
837
GetStructureProteins(vector<const Sequence * > * chains) const838 bool AlignmentManager::GetStructureProteins(vector < const Sequence * > *chains) const
839 {
840 if (!chains || GetCurrentMultipleAlignment() != NULL ||
841 !sequenceViewer || !sequenceViewer->GetCurrentDisplay()) return false;
842
843 sequenceViewer->GetCurrentDisplay()->GetProteinSequences(chains);
844 return (chains->size() > 0);
845 }
846
ReplaceUpdatesInASN(ncbi::objects::CCdd::TPending & newUpdates) const847 void AlignmentManager::ReplaceUpdatesInASN(ncbi::objects::CCdd::TPending& newUpdates) const
848 {
849 if (sequenceSet)
850 sequenceSet->parentSet->ReplaceUpdates(newUpdates);
851 else
852 ERRORMSG("AlignmentManager::ReplaceUpdatesInASN() - can't get StructureSet");
853 }
854
PurgeSequence(const MoleculeIdentifier * identifier)855 void AlignmentManager::PurgeSequence(const MoleculeIdentifier *identifier)
856 {
857 BlockMultipleAlignment *multiple =
858 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
859 sequenceViewer->GetCurrentAlignments().front() : NULL;
860 if (!multiple) return;
861
862 // remove matching rows from multiple alignment
863 vector < unsigned int > rowsToRemove;
864 for (unsigned int i=1; i<multiple->NRows(); ++i)
865 if (multiple->GetSequenceOfRow(i)->identifier == identifier)
866 rowsToRemove.push_back(i);
867
868 if (rowsToRemove.size() > 0) {
869
870 // turn on editor, and update multiple pointer
871 if (!sequenceViewer->EditorIsOn()) {
872 sequenceViewer->TurnOnEditor();
873 multiple = sequenceViewer->GetCurrentAlignments().front();
874 }
875
876 if (!multiple->ExtractRows(rowsToRemove, NULL)) {
877 ERRORMSG("AlignmentManager::PurgeSequence() - ExtractRows failed!");
878 return;
879 }
880
881 // remove rows from SequenceDisplay
882 SequenceDisplay *display = sequenceViewer->GetCurrentDisplay();
883 if (!display) {
884 ERRORMSG("AlignmentManager::PurgeSequence() - can't get SequenceDisplay!");
885 return;
886 }
887 display->RowsRemoved(rowsToRemove, multiple);
888 }
889
890 // remove matching alignments from Update window
891 const ViewerBase::AlignmentList& currentUpdates = updateViewer->GetCurrentAlignments();
892 if (currentUpdates.size() == 0) return;
893 ViewerBase::AlignmentList::const_iterator u, ue = currentUpdates.end();
894
895 for (u=currentUpdates.begin(); u!=ue; ++u) // quick check if any match found
896 if ((*u)->GetSequenceOfRow(1)->identifier == identifier) break;
897
898 if (u != ue) {
899 ViewerBase::AlignmentList updatesToKeep;
900 for (u=currentUpdates.begin(); u!=ue; ++u) {
901 if ((*u)->GetSequenceOfRow(1)->identifier != identifier) {
902 BlockMultipleAlignment *keep = (*u)->Clone();
903 updatesToKeep.push_back(keep);
904 }
905 }
906 updateViewer->ReplaceAlignments(updatesToKeep);
907 }
908 }
909
BlockAlignAllUpdates(void)910 void AlignmentManager::BlockAlignAllUpdates(void)
911 {
912 BlockMultipleAlignment *currentMultiple =
913 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
914 sequenceViewer->GetCurrentAlignments().front() : NULL;
915 if (!currentMultiple) return;
916
917 const UpdateViewer::AlignmentList& currentUpdates = updateViewer->GetCurrentAlignments();
918 if (currentUpdates.size() == 0)
919 return;
920
921 // run the block aligner on update pairwise alignments
922 BlockAligner::AlignmentList newAlignmentsList;
923 int nRowsAddedToMultiple;
924
925 if (blockAligner->CreateNewPairwiseAlignmentsByBlockAlignment(currentMultiple,
926 currentUpdates, &newAlignmentsList, &nRowsAddedToMultiple, sequenceViewer)) {
927
928 // replace update alignments with new ones (or leftovers where algorithm failed)
929 updateViewer->ReplaceAlignments(newAlignmentsList);
930
931 // tell the sequenceViewer that rows have been merged into the multiple
932 if (nRowsAddedToMultiple > 0)
933 sequenceViewer->GetCurrentDisplay()->RowsAdded(nRowsAddedToMultiple, currentMultiple);
934 }
935 }
936
BlockAlignUpdate(BlockMultipleAlignment * single)937 void AlignmentManager::BlockAlignUpdate(BlockMultipleAlignment *single)
938 {
939 BlockMultipleAlignment *currentMultiple =
940 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
941 sequenceViewer->GetCurrentAlignments().front() : NULL;
942 if (!currentMultiple) return;
943
944 // run the block aligner on the given alignment
945 UpdateViewer::AlignmentList singleList, replacedList;
946 BlockAligner::AlignmentList newAlignments;
947 int nRowsAddedToMultiple;
948 bool foundSingle = false; // sanity check
949
950 singleList.push_back(single);
951 if (blockAligner->CreateNewPairwiseAlignmentsByBlockAlignment(
952 currentMultiple, singleList, &newAlignments, &nRowsAddedToMultiple, sequenceViewer)) {
953
954 // replace threaded alignment with new one
955 UpdateViewer::AlignmentList::const_iterator a, ae = updateViewer->GetCurrentAlignments().end();
956 for (a=updateViewer->GetCurrentAlignments().begin(); a!=ae; ++a) {
957 if (*a == single) {
958 BlockAligner::AlignmentList::const_iterator n, ne = newAlignments.end();
959 for (n=newAlignments.begin(); n!=ne; ++n)
960 replacedList.push_back(*n);
961 foundSingle = true;
962 } else
963 replacedList.push_back((*a)->Clone());
964 }
965 if (!foundSingle) ERRORMSG(
966 "AlignmentManager::BlockAlignUpdate() - changed alignment not found in update viewer!");
967 updateViewer->ReplaceAlignments(replacedList);
968
969 // tell the sequenceViewer that rows have been merged into the multiple
970 if (nRowsAddedToMultiple > 0)
971 sequenceViewer->GetCurrentDisplay()->RowsAdded(nRowsAddedToMultiple, currentMultiple);
972 }
973 }
974
975 typedef struct {
976 unsigned int multBlock, pairBlock;
977 bool extendPairBlockLeft;
978 int nResidues;
979 } ExtendInfo;
980
CreateNewPairwiseAlignmentsByBlockExtension(const BlockMultipleAlignment & multiple,const UpdateViewer::AlignmentList & toExtend,UpdateViewer::AlignmentList * newAlignments)981 static bool CreateNewPairwiseAlignmentsByBlockExtension(const BlockMultipleAlignment& multiple,
982 const UpdateViewer::AlignmentList& toExtend, UpdateViewer::AlignmentList *newAlignments)
983 {
984 if (multiple.HasNoAlignedBlocks())
985 return false;
986
987 bool anyChanges = false;
988
989 BlockMultipleAlignment::UngappedAlignedBlockList multBlocks, pairBlocks;
990 multiple.GetUngappedAlignedBlocks(&multBlocks);
991
992 newAlignments->clear();
993 UpdateViewer::AlignmentList::const_iterator t, te = toExtend.end();
994 for (t=toExtend.begin(); t!=te; ++t) {
995 BlockMultipleAlignment *p = (*t)->Clone();
996 newAlignments->push_back(p);
997
998 pairBlocks.clear();
999 p->GetUngappedAlignedBlocks(&pairBlocks);
1000 if (pairBlocks.size() == 0)
1001 continue;
1002
1003 typedef list < ExtendInfo > ExtendList;
1004 typedef map < const UnalignedBlock *, ExtendList > ExtendMap;
1005 ExtendMap extendMap;
1006
1007 // find cases where a pairwise block ends inside a multiple block
1008 const UnalignedBlock *ub;
1009 for (unsigned int mb=0; mb<multBlocks.size(); ++mb) {
1010 const Block::Range *multRange = multBlocks[mb]->GetRangeOfRow(0);
1011 for (unsigned int pb=0; pb<pairBlocks.size(); ++pb) {
1012 const Block::Range *pairRange = pairBlocks[pb]->GetRangeOfRow(0);
1013
1014 if (pairRange->from > multRange->from && pairRange->from <= multRange->to &&
1015 (ub = p->GetUnalignedBlockBefore(pairBlocks[pb])) != NULL)
1016 {
1017 ExtendList& el = extendMap[ub];
1018 el.resize(el.size() + 1);
1019 el.back().multBlock = mb;
1020 el.back().pairBlock = pb;
1021 el.back().extendPairBlockLeft = true;
1022 el.back().nResidues = pairRange->from - multRange->from;
1023 TRACEMSG("block " << (pb+1) << " wants to be extended to the left");
1024 }
1025
1026 if (pairRange->to < multRange->to && pairRange->to >= multRange->from &&
1027 (ub = p->GetUnalignedBlockAfter(pairBlocks[pb])) != NULL)
1028 {
1029 ExtendList& el = extendMap[ub];
1030 el.resize(el.size() + 1);
1031 el.back().multBlock = mb;
1032 el.back().pairBlock = pb;
1033 el.back().extendPairBlockLeft = false;
1034 el.back().nResidues = multRange->to - pairRange->to;
1035 TRACEMSG("block " << (pb+1) << " wants to be extended to the right");
1036 }
1037 }
1038 }
1039
1040 ExtendMap::const_iterator u, ue = extendMap.end();
1041 for (u=extendMap.begin(); u!=ue; ++u) {
1042
1043 // don't slurp residues into non-adjacent aligned blocks from the multiple
1044 if (u->second.size() == 2 && (u->second.back().multBlock - u->second.front().multBlock > 1)) {
1045 TRACEMSG("can't extend with intervening block(s) in multiple between blocks "
1046 << (u->second.front().pairBlock+1) << " and " << (u->second.back().pairBlock+1));
1047 continue;
1048 }
1049
1050 ExtendList::const_iterator e, ee;
1051
1052 // extend from one side only if an unaligned block is contained entirely within a block in the multiple,
1053 // and this unaligned block is "complete" - that is, MinResidues == width.
1054 if (u->first->MinResidues() == u->first->width && u->second.size() == 2 &&
1055 u->second.front().multBlock == u->second.back().multBlock)
1056 {
1057 ExtendList& modList = const_cast<ExtendList&>(u->second);
1058 modList.front().nResidues = u->first->width;
1059 modList.erase(++(modList.begin()));
1060 }
1061
1062 else {
1063 // check to see if there are sufficient unaligned residues to extend from both sides
1064 int available = u->first->MinResidues(), totalShifts = 0;
1065 for (e=u->second.begin(), ee=u->second.end(); e!=ee; ++e)
1066 totalShifts += e->nResidues;
1067 if (totalShifts > available) {
1068 TRACEMSG("inadequate residues to the "
1069 << (u->second.front().extendPairBlockLeft ? "left" : "right")
1070 << " of block " << (u->second.front().pairBlock+1) << "; no extension performed");
1071 continue;
1072 }
1073 }
1074
1075 // perform any allowed extensions
1076 for (e=u->second.begin(), ee=u->second.end(); e!=ee; ++e) {
1077 int alnIdx = p->GetAlignmentIndex(0,
1078 (e->extendPairBlockLeft ?
1079 pairBlocks[e->pairBlock]->GetRangeOfRow(0)->from :
1080 pairBlocks[e->pairBlock]->GetRangeOfRow(0)->to),
1081 BlockMultipleAlignment::eLeft); // justification is irrelevant since this is an aligned block
1082 TRACEMSG("extending " << (e->extendPairBlockLeft ? "left" : "right")
1083 << " side of block " << (e->pairBlock+1) << " by " << e->nResidues << " residues");
1084 if (p->MoveBlockBoundary(alnIdx, alnIdx + (e->extendPairBlockLeft ? -e->nResidues : e->nResidues)))
1085 anyChanges = true;
1086 else
1087 ERRORMSG("MoveBlockBoundary() failed!");
1088 }
1089 }
1090
1091 // now, merge any adjacent blocks that end within a block of the multiple
1092 const Block::Range *prevRange = NULL;
1093 vector < int > mergeIndexes;
1094 for (unsigned int pb=0; pb<pairBlocks.size(); ++pb) {
1095 const Block::Range *pairRange = pairBlocks[pb]->GetRangeOfRow(0);
1096 if (prevRange && prevRange->to == pairRange->from - 1 && multiple.IsAligned(0U, prevRange->to)) {
1097 // justification is irrelevant since this is an aligned block
1098 int mAlnIdx1 = multiple.GetAlignmentIndex(0, prevRange->to, BlockMultipleAlignment::eLeft),
1099 mAlnIdx2 = multiple.GetAlignmentIndex(0, pairRange->from, BlockMultipleAlignment::eLeft),
1100 pAlnIdx1 = p->GetAlignmentIndex(0, prevRange->to, BlockMultipleAlignment::eLeft),
1101 pAlnIdx2 = p->GetAlignmentIndex(0, pairRange->from, BlockMultipleAlignment::eLeft);
1102 if (pAlnIdx1 == pAlnIdx2 - 1 &&
1103 multiple.GetAlignedBlockNumber(mAlnIdx1) == multiple.GetAlignedBlockNumber(mAlnIdx2))
1104 {
1105 TRACEMSG("merging blocks " << pb << " and " << (pb+1));
1106 // just flag to merge later, since actually merging would mess up pb block indexes
1107 mergeIndexes.push_back(pAlnIdx1);
1108 }
1109 }
1110 prevRange = pairRange;
1111 }
1112 for (unsigned int i=0; i<mergeIndexes.size(); ++i) {
1113 if (p->MergeBlocks(mergeIndexes[i], mergeIndexes[i] + 1))
1114 anyChanges = true;
1115 else
1116 ERRORMSG("MergeBlocks() failed!");
1117 }
1118 }
1119
1120 return anyChanges;
1121 }
1122
ExtendAllUpdates(void)1123 void AlignmentManager::ExtendAllUpdates(void)
1124 {
1125 BlockMultipleAlignment *currentMultiple =
1126 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
1127 sequenceViewer->GetCurrentAlignments().front() : NULL;
1128 if (!currentMultiple) return;
1129
1130 const UpdateViewer::AlignmentList& currentUpdates = updateViewer->GetCurrentAlignments();
1131 if (currentUpdates.size() == 0)
1132 return;
1133
1134 // run the block extender on update pairwise alignments
1135 UpdateViewer::AlignmentList newAlignmentsList;
1136 if (CreateNewPairwiseAlignmentsByBlockExtension(*currentMultiple, currentUpdates, &newAlignmentsList))
1137 updateViewer->ReplaceAlignments(newAlignmentsList);
1138 }
1139
ExtendUpdate(BlockMultipleAlignment * single)1140 void AlignmentManager::ExtendUpdate(BlockMultipleAlignment *single)
1141 {
1142 BlockMultipleAlignment *currentMultiple =
1143 (sequenceViewer->GetCurrentAlignments().size() > 0) ?
1144 sequenceViewer->GetCurrentAlignments().front() : NULL;
1145 if (!currentMultiple) return;
1146
1147 // run the threader on the given alignment
1148 UpdateViewer::AlignmentList singleList, replacedList;
1149 BlockAligner::AlignmentList newAlignments;
1150
1151 singleList.push_back(single);
1152 if (CreateNewPairwiseAlignmentsByBlockExtension(*currentMultiple, singleList, &newAlignments)) {
1153 if (newAlignments.size() != 1) {
1154 ERRORMSG("AlignmentManager::ExtendUpdate() - returned alignment list size != 1!");
1155 return;
1156 }
1157
1158 // replace original alignment with extended one
1159 UpdateViewer::AlignmentList::const_iterator a, ae = updateViewer->GetCurrentAlignments().end();
1160 bool foundSingle = false; // sanity check
1161 for (a=updateViewer->GetCurrentAlignments().begin(); a!=ae; ++a) {
1162 if (*a == single) {
1163 replacedList.push_back(newAlignments.front());
1164 foundSingle = true;
1165 } else
1166 replacedList.push_back((*a)->Clone());
1167 }
1168 if (!foundSingle)
1169 ERRORMSG("AlignmentManager::ExtendUpdate() - changed alignment not found in update viewer!");
1170 updateViewer->ReplaceAlignments(replacedList);
1171 }
1172 }
1173
1174 ProgressMeter *g_refinerProgressMeter = NULL;
1175 #define PROGRESS_MAX 1000
1176
RefinementProgress(double progress)1177 void RefinementProgress(double progress)
1178 {
1179 if (g_refinerProgressMeter)
1180 g_refinerProgressMeter->SetValue((int) (progress * PROGRESS_MAX));
1181 }
1182
RefineAlignment(bool setUpOptionsOnly)1183 void AlignmentManager::RefineAlignment(bool setUpOptionsOnly)
1184 {
1185 // TODO: selection dialog for realigning specific rows
1186 // <can add this afterwards; refiner supports this sort of thing but need to improve api>
1187
1188 // get info on current multiple alignment
1189 const BlockMultipleAlignment *multiple = GetCurrentMultipleAlignment();
1190 if (!multiple || multiple->NRows() < 2 || multiple->HasNoAlignedBlocks()) {
1191 ERRORMSG("AlignmentManager::RefineAlignment() - can't refine alignments with < 2 rows or zero aligned blocks");
1192 return;
1193 }
1194
1195 // construct CDD from inputs
1196 CRef < CCdd > cdd(new CCdd);
1197 cdd->SetName("CDD input to refiner from Cn3D");
1198 CRef < CCdd_id > uid(new CCdd_id);
1199 uid->SetUid(0);
1200 cdd->SetId().Set().push_back(uid);
1201
1202 // construct Seq-entry from sequences in current alignment, starting with master
1203 CRef < CSeq_entry > seq(new CSeq_entry);
1204 seq->SetSeq().Assign(multiple->GetMaster()->bioseqASN.GetObject());
1205 cdd->SetSequences().SetSet().SetSeq_set().push_back(seq);
1206
1207 // construct Seq-annot from rows in the alignment
1208 CRef < CSeq_annot > seqAnnot(new CSeq_annot);
1209 cdd->SetSeqannot().push_back(seqAnnot);
1210 BlockMultipleAlignment::UngappedAlignedBlockList blocks;
1211 multiple->GetUngappedAlignedBlocks(&blocks);
1212 vector < unsigned int > rowOrder;
1213 sequenceViewer->GetCurrentDisplay()->GetRowOrder(multiple, &rowOrder);
1214
1215 // other info used by the refiner
1216 vector < string > rowTitles(multiple->NRows());
1217 vector < bool > rowIsStructured(multiple->NRows());
1218 vector < bool > blocksToRealign;
1219
1220 // fill out Seq-entry and Seq-annot based on current row ordering of the display (which may be different from BMA)
1221 for (unsigned int i=1; i<multiple->NRows(); ++i) {
1222 seq.Reset(new CSeq_entry);
1223 seq->SetSeq().Assign(multiple->GetSequenceOfRow(rowOrder[i])->bioseqASN.GetObject());
1224 cdd->SetSequences().SetSet().SetSeq_set().push_back(seq);
1225 CRef < CSeq_align > seqAlign(CreatePairwiseSeqAlignFromMultipleRow(multiple, blocks, rowOrder[i]));
1226 seqAnnot->SetData().SetAlign().push_back(seqAlign);
1227 rowTitles[i] = multiple->GetSequenceOfRow(rowOrder[i])->identifier->ToString();
1228 rowIsStructured[i] = (multiple->GetSequenceOfRow(rowOrder[i])->identifier->pdbID.size() > 0);
1229 }
1230
1231 // set blocks to realign, using marked blocks, if any; no marks -> realign all blocks
1232 vector < unsigned int > marks;
1233 multiple->GetMarkedBlockNumbers(&marks);
1234 blocksToRealign.resize(blocks.size(), (marks.size() == 0));
1235 for (unsigned int i=0; i<marks.size(); ++i) {
1236 TRACEMSG("refining block " << (marks[i]+1));
1237 blocksToRealign[marks[i]] = true;
1238 }
1239
1240 // set up refiner
1241 align_refine::BMARefinerInterface interface;
1242 static auto_ptr < align_refine::BMARefinerOptions > options;
1243 if ((options.get() && !interface.SetOptions(*options)) || // set these first since some are overridden by alignment inputs
1244 !interface.SetInitialAlignment(*cdd, blocks.size(), multiple->NRows()) ||
1245 !interface.SetRowTitles(rowTitles) ||
1246 !interface.SetRowsWithStructure(rowIsStructured) ||
1247 !interface.SetBlocksToRealign(blocksToRealign))
1248 {
1249 ERRORMSG("AlignmentManager::RefineAlignment() - failed to set up BMARefinerInterface");
1250 return;
1251 }
1252 if (!options.get() || setUpOptionsOnly) {
1253 if (!align_refine::BMARefinerOptionsDialog::SetRefinerOptionsViaDialog(NULL, interface))
1254 {
1255 WARNINGMSG("AlignmentManager::RefineAlignment() - failed to set refiner options via dialog, probably cancelled");
1256 return;
1257 }
1258 options.reset(interface.GetOptions());
1259 }
1260 if (setUpOptionsOnly)
1261 return;
1262
1263 // set up a progress meter
1264 g_refinerProgressMeter = new ProgressMeter(NULL, "Refinement in progress...", "Wait", PROGRESS_MAX);
1265 g_refinerProgressMeter->SetValue(0);
1266
1267 // actually run the refiner algorithm
1268 CCdd::TSeqannot results;
1269 SetDialogSevereErrors(false);
1270 wxSetCursor(*wxHOURGLASS_CURSOR);
1271 bool okay = interface.Run(results, RefinementProgress);
1272 wxSetCursor(wxNullCursor);
1273 SetDialogSevereErrors(true);
1274 SetDiagPostLevel(eDiag_Info); // may be changed by refiner
1275
1276 // delete progress meter
1277 g_refinerProgressMeter->Close(true);
1278 g_refinerProgressMeter->Destroy();
1279 g_refinerProgressMeter = NULL;
1280
1281 if (okay) {
1282
1283 // unpack results
1284 if (results.size() > 0 && results.front().NotEmpty()) {
1285
1286 // just use first returned alignment for now; convert asn data into BlockMultipleAlignment
1287 PairwiseAlignmentList pairs;
1288 if (!results.front()->IsSetData() || !results.front()->GetData().IsAlign()) {
1289 WARNINGMSG("AlignmentManager::RefineAlignment() - got result Seq-annot in unexpected format");
1290 } else {
1291 CSeq_annot::C_Data::TAlign::const_iterator l, le = results.front()->GetData().GetAlign().end();
1292 for (l=results.front()->GetData().GetAlign().begin(); l!=le; ++l)
1293 pairs.push_back(new MasterDependentAlignment(NULL, multiple->GetMaster(), **l));
1294 }
1295 auto_ptr < BlockMultipleAlignment > refined(CreateMultipleFromPairwiseWithIBM(pairs));
1296 DELETE_ALL_AND_CLEAR(pairs, PairwiseAlignmentList);
1297
1298 // feed result back into alignment window
1299 if (refined.get() && refined->NRows() == multiple->NRows())
1300 sequenceViewer->ReplaceAlignment(multiple, refined.release());
1301 else
1302 ERRORMSG("AlignmentManager::RefineAlignment() - problem converting refinement result");
1303
1304 } else {
1305 WARNINGMSG("AlignmentManager::RefineAlignment() - failed to make a refined alignment. Alignment unchanged.");
1306 }
1307
1308 } else {
1309 ERRORMSG("AlignmentManager::RefineAlignment() - refinement failed. Alignment unchanged.");
1310 }
1311 }
1312
1313 END_SCOPE(Cn3D)
1314