ARB
SEC_split.cxx
Go to the documentation of this file.
1 // =============================================================== //
2 // //
3 // File : SEC_split.cxx //
4 // Purpose : split/unsplit loops (aka fold/unfold helices) //
5 // //
6 // Coded by Ralf Westram (coder@reallysoft.de) in August 2007 //
7 // Institute of Microbiology (Technical University Munich) //
8 // http://www.arb-home.de/ //
9 // //
10 // =============================================================== //
11 
12 #include "SEC_root.hxx"
13 #include "SEC_drawn_pos.hxx"
14 #include "SEC_iter.hxx"
15 #include <arb_msg.h>
16 
17 using namespace std;
18 
22 };
23 
24 class AngleBuffer { // stores the absolute values of some SEC_oriented
25  typedef std::map<SEC_helix*, Angle> AngleMap;
26 
27  AngleMap angles;
28  AngleBufferMode mode;
29 
30  Angle loop2helix(SEC_loop *loop, SEC_helix *helix) {
31  return Angle(loop->get_center(), helix->strandAwayFrom(loop)->get_fixpoint());
32  }
33 
34 public:
35  AngleBuffer(AngleBufferMode Mode) : mode(Mode) {}
36 
37  void store(SEC_helix *helix, SEC_loop *loop) {
38  switch (mode) {
40  angles[helix] = helix->get_abs_angle();
41  break;
43  angles[helix] = helix->get_abs_angle()-loop2helix(loop, helix);
44  break;
45  }
46  }
47 
48  void set_angle(SEC_helix *helix, const Angle& angle) { angles[helix] = angle; }
49 
50  void restore(SEC_helix *helix, SEC_loop *loop) {
51  switch (mode) {
53  helix->set_abs_angle(angles[helix]);
54  break;
56  if (helix->hasLoop(loop)) {
57  helix->set_abs_angle(angles[helix]+loop2helix(loop, helix));
58  }
59  break;
60  }
61  }
62 
63  void restoreAll(SEC_loop *loop) {
64  AngleMap::iterator e = angles.end();
65  for (AngleMap::iterator a = angles.begin(); a != e; ++a) {
66  restore(a->first, loop);
67  }
68  }
69 
70  void remove(SEC_helix *helix) { angles.erase(helix); }
71 
72  void storeAllHelices(SEC_loop *loop, SEC_helix *skip) {
73  for (SEC_strand_iterator strand(loop); strand; ++strand) {
74  SEC_helix *helix = strand->get_helix();
75  if (helix != skip) store(helix, loop);
76  }
77  }
78 };
79 
80 
81 // ---------------------
82 // moving root
83 
84 void SEC_loop::toggle_root(SEC_loop *old_root) {
85  // make this the new root loop
86  sec_assert(old_root != this);
87  SEC_helix *mid_helix = get_rootside_helix();
88 
89  // set root to loop behind mid_helix
90  {
91  SEC_loop *behind = mid_helix->rootsideLoop();
92  if (behind != old_root) {
93  behind->toggle_root(old_root);
94  old_root = behind;
95  }
96  }
97 
98  // store abs angles of all strands
99  Angle midAngle = mid_helix->get_abs_angle();
100 
102  AngleBuffer otherOldAbs(BUFFER_ABSOLUTE_ANGLES);
103 
104  thisOldAbs.storeAllHelices(this, mid_helix);
105  otherOldAbs.storeAllHelices(old_root, mid_helix);
106 
107  // modify structure
108  get_root()->set_root_loop(this);
109  mid_helix->flip();
110  set_fixpoint_strand(mid_helix->strandAwayFrom(this));
111  old_root->set_fixpoint_strand(mid_helix->strandAwayFrom(old_root));
112 
113  // calculate abs angles of loops and mid_helix
114  set_rel_angle(Angle(center, get_fixpoint()));
115  mark_angle_absolute(); // root-loop: rel == abs
116  mid_helix->set_abs_angle(midAngle.rotate180deg());
117  old_root->set_abs_angle(Angle(old_root->get_fixpoint(), old_root->get_center()));
118 
119  // restore angles of other helices
120  thisOldAbs.restoreAll(this);
121  otherOldAbs.restoreAll(old_root);
122 }
123 
124 void SEC_root::set_root(SEC_loop *loop) {
125  SEC_loop *old_root = get_root_loop();
126  if (loop != old_root) {
127  Vector new2old(loop->get_center(), old_root->get_center());
128  add_autoscroll(new2old);
129  loop->toggle_root(old_root);
130  recalc();
131  }
132 }
133 
134 
135 // ----------------------------------
136 // search segment by abspos
137 
139  SEC_helix_strand *start_strand = root_loop->get_fixpoint_strand();
140  SEC_helix_strand *strand = start_strand;
141  do {
142  SEC_region *reg = strand->get_region();
143  if (reg->contains_seq_position(pos)) return strand;
144 
145  SEC_helix_strand *other_strand = strand->get_other_strand();
146  SEC_region *oreg = other_strand->get_region();
147 
148  SEC_segment *seg;
149  if (SEC_region(reg->get_sequence_end(), oreg->get_sequence_start()).contains_seq_position(pos)) {
150  seg = other_strand->get_next_segment();
151  }
152  else {
153  if (oreg->contains_seq_position(pos)) return other_strand;
154  seg = strand->get_next_segment();
155  }
156 
157  if (seg->get_region()->contains_seq_position(pos)) return seg;
158  strand = seg->get_next_strand();
159  }
160  while (strand != start_strand);
161 
162  return NULp;
163 }
164 
165 inline SEC_segment *findSegmentContaining(SEC_root *root, int pos, GB_ERROR& error) {
166  SEC_segment *result = NULp;
167  error = NULp;
168 
169  SEC_base_part *found = root->find(pos);
170  if (found) {
171  if (found->parent()->getType() == SEC_LOOP) {
172  result = static_cast<SEC_segment*>(found);
173  }
174  else {
175  error = GBS_global_string("Position %i not in a segment", pos);
176  }
177  }
178  else {
179  error = GBS_global_string("Position %i is outside allowed range", pos);
180  }
181  return result;
182 }
183 
184 inline SEC_segment *findSegmentContaining(SEC_root *root, int start, int end, GB_ERROR& error) {
185  // end is position behind questionable position
186  error = NULp;
187 
188  SEC_segment *start_segment = findSegmentContaining(root, start, error);
189  if (start_segment) {
190  SEC_segment *end_segment;
191  if (end == start+1) {
192  end_segment = start_segment;
193  }
194  else {
195  end_segment = findSegmentContaining(root, end-1, error);
196  }
197 
198  if (end_segment) {
199  if (end_segment != start_segment) {
200  error = GBS_global_string("Positions %i and %i are in different segments", start, end);
201  start_segment = NULp;
202  }
203  }
204  }
205  sec_assert(contradicted(start_segment, error));
206  return start_segment;
207 }
208 
209 // --------------------
210 // split loop
211 
212 GB_ERROR SEC_root::split_loop(int start1, int end1, int start2, int end2) {
213  // end1/end2 are positions behind the helix-positions!
214  sec_assert(start1<end1);
215  sec_assert(start2<end2);
216 
217  if (start1>start2) {
218  return split_loop(start2, end2, start1, end1);
219  }
220 
221  GB_ERROR error = NULp;
222  if (start2<end1) {
223  error = GBS_global_string("Helices overlap (%i-%i and %i-%i)", start1, end1, start2, end2);
224  }
225 
226  if (!error) {
227  SEC_segment *seg1 = findSegmentContaining(this, start1, end1, error);
228  SEC_segment *seg2 = NULp;
229 
230  if (!error) seg2 = findSegmentContaining(this, start2, end2, error);
231 
232  if (!error) {
233  sec_assert(seg1 && seg2);
234  SEC_loop *old_loop = seg1->get_loop();
235  if (old_loop != seg2->get_loop()) {
236  error = "Positions are in different loops (no tertiary structures possible)";
237  }
238  else {
239  SEC_loop *setRootTo = NULp; // set root back afterwards?
240 
241  if (old_loop->is_root_loop()) {
242  set_root(seg1->get_next_strand()->get_destination_loop()); // another loop
243  setRootTo = old_loop;
244  }
245 
247  oldAngles.storeAllHelices(old_loop, NULp);
248 
249  SEC_helix *new_helix = NULp;
250  SEC_loop *new_loop = NULp;
251 
252  if (seg1 == seg2) { // split one segment
253  // \ .
254  // \ seg1 >>> .
255  // seg1 \ strand1 ....
256  // ______________ => \_________________. . seg2
257  // _________________. .
258  // / ....
259  // / strand2
260  // / seg3 <<<
261  // /
262  //
263  // seg1 is the old segment
264 
265  SEC_helix_strand *strand1 = seg1->split(start1, end1, &seg2);
266  SEC_helix_strand *strand2 = NULp;
267  SEC_segment *seg3 = NULp;
268 
269  if (seg1->get_region()->contains_seq_position(start2)) {
270  seg3 = seg2;
271  strand2 = strand1;
272  strand1 = seg1->split(start2, end2, &seg2);
273  }
274  else {
275  sec_assert(seg2->get_region()->contains_seq_position(start2));
276  strand2 = seg2->split(start2, end2, &seg3);
277  }
278 
279  sec_assert(are_adjacent_regions(seg1->get_region(), strand1->get_region()));
280  sec_assert(are_adjacent_regions(strand1->get_region(), seg2->get_region()));
281  sec_assert(are_adjacent_regions(seg2->get_region(), strand2->get_region()));
282  sec_assert(are_adjacent_regions(strand2->get_region(), seg3->get_region()));
283 
284  new_helix = new SEC_helix(this, strand2, strand1); // strands are responsible for memory
285 
286  strand1->set_next_segment(seg3);
287  strand2->set_next_segment(seg2);
288 
289  new_loop = new SEC_loop(this);
290 
291  seg2->set_loop(new_loop);
292 
293  strand2->set_origin_loop(new_loop);
294  new_loop->set_fixpoint_strand(strand2);
295  }
296  else { // split two segments
297  // \ /
298  // seg1 \ seg1 >>> /
299  // ______________ \ strand1 / seg3
300  // \_________________/
301  // => old_loop _________________ new_loop
302  // ______________ / \ .
303  // seg2 / strand2 \ seg2
304  // / seg4 <<< \ .
305  // / \ .
306  //
307  // seg1 and seg2 are the old segments
308 
309  // maybe swap seg1/seg2 (to ensure fixpoint-strand stays in old loop)
310  for (SEC_segment *s = seg1; s != seg2;) {
311  SEC_helix_strand *hs = s->get_next_strand();
312  if (!hs->isRootsideFixpoint()) { // fixpoint-strand is between seg1 -> seg2
313  // swap seg1<->seg2
314  swap(seg1, seg2);
315  swap(start1, start2);
316  swap(end1, end2);
317  break;
318  }
319  s = hs->get_next_segment();
320  }
321 
322  SEC_segment *seg3;
323  SEC_segment *seg4;
324  SEC_helix_strand *strand1 = seg1->split(start1, end1, &seg3);
325  SEC_helix_strand *strand2 = seg2->split(start2, end2, &seg4);
326 
327  new_helix = new SEC_helix(this, strand2, strand1);
328 
329  strand1->set_next_segment(seg4);
330  strand2->set_next_segment(seg3);
331 
332  new_loop = new SEC_loop(this);
333 
334  for (SEC_segment *s = seg3; ;) {
335  s->set_loop(new_loop);
336  SEC_helix_strand *h = s->get_next_strand();
337  h->set_origin_loop(new_loop);
338  if (s == seg2) break;
339  s = h->get_next_segment();
340  }
341 
342  new_loop->set_fixpoint_strand(strand2);
343  }
344 
345  // set angles of new helix and new loop
346  new_helix->set_rel_angle(0); // wrong, but relayout fails otherwise
347  new_loop->set_rel_angle(0);
348 
349  relayout();
350 
351  // correct angles of other helices
352  oldAngles.restoreAll(new_loop);
353  oldAngles.set_angle(new_helix, Angle(0));
354  oldAngles.restoreAll(old_loop);
355 
356  recalc();
357 
358  if (setRootTo) {
359  set_root(setRootTo); // restore root loop
360  }
361  }
362  }
363  }
364  return error;
365 }
366 
367 // -----------------------
368 // fold a strand
369 
370 GB_ERROR SEC_root::unsplit_loop(SEC_helix_strand *remove_strand) {
371  //
372  // \ before[0] / after[1]
373  // \ /
374  // \ >>>> /
375  // \ strand[0] /
376  // \_________________/
377  // loop[0] _________________ loop[1]
378  // / strand[1] \ .
379  // / <<<< \ .
380  // / \ .
381  // / \ .
382  // / after[0] \ before[1]
383  //
384  // The strands are removed and segments get connected.
385  // One loop is deleted.
386 
387  GB_ERROR error = NULp;
388  SEC_helix_strand *strand[2] = { remove_strand, remove_strand->get_other_strand() };
389  SEC_segment *before[2], *after[2];
390  SEC_loop *loop[2];
391 
392 #if defined(CHECK_INTEGRITY)
393  check_integrity(CHECK_STRUCTURE);
394 #endif // CHECK_INTEGRITY
395 
396  int s;
397  for (s = 0; s<2; s++) {
398  after[s] = strand[s]->get_next_segment();
399  before[s] = strand[s]->get_previous_segment();
400  loop[s] = strand[s]->get_origin_loop();
401 
402  sec_assert(before[s]->get_loop() == loop[s]);
403  sec_assert(after[s]->get_loop() == loop[s]);
404  }
405 
406  bool is_terminal_loop[2] = { before[0] == after[0], before[1] == after[1] };
407  int i0 = -1; // index of terminal loop (or -1)
408  bool unsplit = true;
409 
410  if (is_terminal_loop[0]) {
411  if (is_terminal_loop[1]) {
412  error = "You cannot delete the last helix";
413  unsplit = false;
414  }
415  else i0 = 0;
416  }
417  else {
418  if (is_terminal_loop[1]) i0 = 1;
419  }
420 
421  if (unsplit) {
422  int del = i0 >= 0 ? i0 : 1; // index of loop which will be deleted
423 
424  SEC_loop *setRootTo = NULp; // set root back to afterwards?
425 
426  {
427  // move away root-loop to make things easy
428  SEC_loop *rootLoop = get_root_loop();
429  if (loop[0] == rootLoop || loop[1] == rootLoop) {
430  SEC_loop *termLoop = is_terminal_loop[0] ? loop[0] : loop[1];
431  SEC_helix_strand *toTerm = strand[0]->get_helix()->strandTowards(termLoop);
432  SEC_helix_strand *toNextLoop = toTerm->get_next_segment()->get_next_strand();
433  SEC_loop *anotherLoop = toNextLoop->get_destination_loop();
434 
435  sec_assert(anotherLoop != loop[0] && anotherLoop != loop[1]);
436 
437  set_root(anotherLoop);
438  setRootTo = loop[1-del]; // afterwards set root back to non-deleted loop
439  }
440  }
441 
442  SEC_helix *removed_helix = strand[0]->get_helix();
444  oldAngles.storeAllHelices(loop[0], removed_helix);
445  oldAngles.storeAllHelices(loop[1], removed_helix);
446 
447  if (i0 >= 0) { // one loop is terminal
448  // i0 and i1 are indexes 0 and 1 in picture above.
449  // The left loop (loop[i0]) will be removed.
450  int i1 = 1-i0; // index of non-terminal loop
451 
452  before[i1]->mergeWith(after[i1], loop[i1]);
453 
454  sec_assert(after[i0] == before[i0]);
455  delete after[i0]; // delete the segment of the terminal-loop
456  }
457  else { // none of the loops is terminal
458  // keep loop[0], delete loop[1]
459 
460  SEC_helix_strand *rootsideStrand = strand[0]->get_helix()->rootsideLoop()->get_rootside_strand();
461 
462  before[1]->mergeWith(after[0], loop[0]);
463  before[0]->mergeWith(after[1], loop[0]);
464  // after[] segments are invalid now!
465 
466  // loop over all segments in loop[1] and relink them to loop[0]
467  SEC_segment *seg = before[0];
468  while (seg != before[1]) {
469  SEC_helix_strand *loop1_strand = seg->get_next_strand();
470  loop1_strand->set_origin_loop(loop[0]);
471 
472  seg = loop1_strand->get_next_segment();
473  seg->set_loop(loop[0]);
474  }
475 
476  loop[0]->set_fixpoint_strand(rootsideStrand);
477  }
478 
479  loop[del]->set_fixpoint_strand(NULp);
480  delete loop[del];
481 
482  for (s = 0; s<2; s++) strand[s]->unlink(false);
483  delete strand[0]; // delete both strands
484 
485  relayout();
486 
487  oldAngles.restoreAll(loop[1-del]);
488  recalc();
489 
490  if (setRootTo) set_root(setRootTo);
491  }
492 
493  return error;
494 }
const char * GB_ERROR
Definition: arb_core.h:25
string result
GB_ERROR unsplit_loop(SEC_helix_strand *delete_strand)
Definition: SEC_split.cxx:370
void restoreAll(SEC_loop *loop)
Definition: SEC_split.cxx:63
int get_sequence_start() const
Definition: SEC_root.hxx:107
GB_ERROR split_loop(int start1, int end1, int start2, int end2)
Definition: SEC_split.cxx:212
int get_sequence_end() const
Definition: SEC_root.hxx:108
bool are_adjacent_regions(const SEC_region *reg1, const SEC_region *reg2)
Definition: SEC_root.hxx:993
void restore(SEC_helix *helix, SEC_loop *loop)
Definition: SEC_split.cxx:50
SEC_segment * findSegmentContaining(SEC_root *root, int pos, GB_ERROR &error)
Definition: SEC_split.cxx:165
SEC_base * parent()
Definition: SEC_root.hxx:295
const char * GBS_global_string(const char *templat,...)
Definition: arb_msg.cxx:204
STL namespace.
bool contains_seq_position(int pos) const
Definition: SEC_root.hxx:112
virtual SEC_BASE_TYPE getType() const =0
Angle & rotate180deg()
static HelixNrInfo * start
AngleBufferMode
Definition: SEC_split.cxx:19
AliDataPtr after(AliDataPtr data, size_t pos)
Definition: insdel.cxx:593
AngleBuffer(AngleBufferMode Mode)
Definition: SEC_split.cxx:35
void store(SEC_helix *helix, SEC_loop *loop)
Definition: SEC_split.cxx:37
static void error(const char *msg)
Definition: mkptypes.cxx:96
CONSTEXPR_INLINE_Cxx14 void swap(unsigned char &c1, unsigned char &c2)
Definition: ad_io_inline.h:19
AliDataPtr before(AliDataPtr data, size_t pos)
Definition: insdel.cxx:592
SEC_base_part * find(int pos)
Definition: SEC_split.cxx:138
void set_angle(SEC_helix *helix, const Angle &angle)
Definition: SEC_split.cxx:48
#define NULp
Definition: cxxforward.h:97
void storeAllHelices(SEC_loop *loop, SEC_helix *skip)
Definition: SEC_split.cxx:72
void set_root(SEC_loop *loop)
Definition: SEC_split.cxx:124
#define sec_assert(cond)
Definition: SEC_defs.hxx:19
GB_write_int const char s
Definition: AW_awar.cxx:156