source: trunk/SECEDIT/SEC_split.cxx

Last change on this file was 16763, checked in by westram, 7 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.8 KB
Line 
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
17using namespace std;
18
19enum AngleBufferMode {
20    BUFFER_ABSOLUTE_ANGLES,
21    BUFFER_CENTER_RELATIVE,
22};
23
24class 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
34public:
35    AngleBuffer(AngleBufferMode Mode) : mode(Mode) {}
36
37    void store(SEC_helix *helix, SEC_loop *loop) {
38        switch (mode) {
39            case BUFFER_ABSOLUTE_ANGLES:
40                angles[helix] = helix->get_abs_angle();
41                break;
42            case BUFFER_CENTER_RELATIVE:
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) {
52            case BUFFER_ABSOLUTE_ANGLES:
53                helix->set_abs_angle(angles[helix]);
54                break;
55            case BUFFER_CENTER_RELATIVE:
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
84void 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
101    AngleBuffer thisOldAbs(BUFFER_ABSOLUTE_ANGLES);
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
124void 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
138SEC_base_part *SEC_root::find(int pos) {
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
165inline 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
184inline 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
212GB_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
246                AngleBuffer oldAngles(BUFFER_CENTER_RELATIVE);
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
370GB_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();
443        AngleBuffer  oldAngles(BUFFER_CENTER_RELATIVE);
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}
Note: See TracBrowser for help on using the repository browser.