source: tags/svn.1.5.4/SECEDIT/SEC_split.cxx

Last change on this file was 7272, checked in by westram, 14 years ago

merged from dev [7244] [7250] [7251] [7252] [7253] [7254] [7255] [7257] [7258] [7259] [7260] [7261] [7271]

  • extracted conversion logic of absolute ↔ relative position into BasePosition
    • refactored BI_ecoli_ref and ED4_base_position (using BasePosition)
    • corrected input range of abs_2_rel ([0..N-1] → [0..N])
  • clean up chaos caused by different semantics of positions (biological[1..N] vs informational[0..N-1])
    • commented (hopefully) all AWARs and variables which use biological positions
    • made semantic changes explicit by using bio2info() and info2bio()
  • positions passed to/reported by PTserver are now biological positions. Affects
    • probe design (absolute position, ecoli range restriction)
    • probe match (absolute position, ecoli position (bugfix see below)).
    • ProbeMatchParser does not change semantics
  • fix several bugs
    • ED4_jump_to_cursor_position limited position to screen before converting from relative to absolute position
    • SEC_root::paintEcoliPositions added 1 to relative pos (instead of including position at cursor)
    • ecoli position in probe match/PTserver (did not include ecoli-base at match position)
    • ecoli range restriction in probe design/PTserver (check was broken, see [7261])
    • base position for consensus terminals
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.9 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// --------------------
84
85void SEC_loop::toggle_root(SEC_loop *old_root) {
86    // make this the new root loop
87    sec_assert(old_root != this);
88    SEC_helix *mid_helix = get_rootside_helix();
89
90    // set root to loop behind mid_helix
91    {
92        SEC_loop *behind = mid_helix->rootsideLoop();
93        if (behind != old_root) {
94            behind->toggle_root(old_root);
95            old_root = behind;
96        }
97    }
98
99    // store abs angles of all strands
100    Angle midAngle = mid_helix->get_abs_angle();
101
102    AngleBuffer thisOldAbs(BUFFER_ABSOLUTE_ANGLES);
103    AngleBuffer otherOldAbs(BUFFER_ABSOLUTE_ANGLES);
104
105    thisOldAbs.storeAllHelices(this, mid_helix);
106    otherOldAbs.storeAllHelices(old_root, mid_helix);
107
108    // modify structure
109    get_root()->set_root_loop(this);
110    mid_helix->flip();
111    set_fixpoint_strand(mid_helix->strandAwayFrom(this));
112    old_root->set_fixpoint_strand(mid_helix->strandAwayFrom(old_root));
113
114    // calculate abs angles of loops and mid_helix
115    set_rel_angle(Angle(center, get_fixpoint()));
116    mark_angle_absolute();      // root-loop: rel == abs
117    mid_helix->set_abs_angle(midAngle.rotate180deg());
118    old_root->set_abs_angle(Angle(old_root->get_fixpoint(), old_root->get_center()));
119
120    // restore angles of other helices
121    thisOldAbs.restoreAll(this);
122    otherOldAbs.restoreAll(old_root);
123}
124
125void SEC_root::set_root(SEC_loop *loop) {
126    SEC_loop *old_root = get_root_loop();
127    if (loop != old_root) {
128        Vector new2old(loop->get_center(), old_root->get_center());
129        add_autoscroll(new2old);
130        loop->toggle_root(old_root);
131        recalc();
132    }
133}
134
135
136// ---------------------------------
137//      search segment by abspos
138// ---------------------------------
139
140SEC_base_part *SEC_root::find(int pos) {
141    SEC_helix_strand *start_strand = root_loop->get_fixpoint_strand();
142    SEC_helix_strand *strand       = start_strand;
143    do {
144        SEC_region *reg = strand->get_region();
145        if (reg->contains_seq_position(pos)) return strand;
146
147        SEC_helix_strand *other_strand = strand->get_other_strand();
148        SEC_region       *oreg         = other_strand->get_region();
149
150        SEC_segment *seg;
151        if (SEC_region(reg->get_sequence_end(), oreg->get_sequence_start()).contains_seq_position(pos)) {
152            seg = other_strand->get_next_segment();
153        }
154        else {
155            if (oreg->contains_seq_position(pos)) return other_strand;
156            seg = strand->get_next_segment();
157        }
158
159        if (seg->get_region()->contains_seq_position(pos)) return seg;
160        strand = seg->get_next_strand();
161    }
162    while (strand != start_strand);
163
164    return 0;
165}
166
167inline SEC_segment *findSegmentContaining(SEC_root *root, int pos, GB_ERROR& error) {
168    SEC_segment *result = 0;
169    error               = 0;
170
171    SEC_base_part *found  = root->find(pos);
172    if (found) {
173        if (found->parent()->getType() == SEC_LOOP) {
174            result = static_cast<SEC_segment*>(found);
175        }
176        else {
177            error = GBS_global_string("Position %i not in a segment", pos);
178        }
179    }
180    else {
181        error = GBS_global_string("Position %i is outside allowed range", pos);
182    }
183    return result;
184}
185
186inline SEC_segment *findSegmentContaining(SEC_root *root, int start, int end, GB_ERROR& error) {
187    // end is position behind questionable position
188    error = 0;
189
190    SEC_segment *start_segment = findSegmentContaining(root, start, error);
191    if (start_segment) {
192        SEC_segment *end_segment;
193        if (end == start+1) {
194            end_segment = start_segment;
195        }
196        else {
197            end_segment = findSegmentContaining(root, end-1, error);
198        }
199
200        if (end_segment) {
201            if (end_segment != start_segment) {
202                error         = GBS_global_string("Positions %i and %i are in different segments", start, end);
203                start_segment = 0;
204            }
205        }
206    }
207    sec_assert(contradicted(start_segment, error));
208    return start_segment;
209}
210
211// -------------------
212//      split loop
213// -------------------
214
215GB_ERROR SEC_root::split_loop(int start1, int end1, int start2, int end2) {
216    // end1/end2 are positions behind the helix-positions!
217    sec_assert(start1<end1);
218    sec_assert(start2<end2);
219
220    if (start1>start2) {
221        return split_loop(start2, end2, start1, end1);
222    }
223
224    GB_ERROR error = 0;
225    if (start2<end1) {
226        error = GBS_global_string("Helices overlap (%i-%i and %i-%i)", start1, end1, start2, end2);
227    }
228
229    if (!error) {
230        SEC_segment *seg1 = findSegmentContaining(this, start1, end1, error);
231        SEC_segment *seg2 = 0;
232
233        if (!error) seg2 = findSegmentContaining(this, start2, end2, error);
234
235        if (!error) {
236            sec_assert(seg1 && seg2);
237            SEC_loop *old_loop = seg1->get_loop();
238            if (old_loop != seg2->get_loop()) {
239                error = "Positions are in different loops (no tertiary structures possible)";
240            }
241            else {
242                SEC_loop *setRootTo = 0; // set root back afterwards?
243
244                if (old_loop->is_root_loop()) {
245                    set_root(seg1->get_next_strand()->get_destination_loop()); // another loop
246                    setRootTo = old_loop;
247                }
248
249                AngleBuffer oldAngles(BUFFER_CENTER_RELATIVE);
250                oldAngles.storeAllHelices(old_loop, 0);
251
252                SEC_helix *new_helix = 0;
253                SEC_loop  *new_loop  = 0;
254
255                if (seg1 == seg2) { // split one segment
256                    //                           \                                                   .
257                    //                            \ seg1     >>>                                     .
258                    //     seg1                    \       strand1     ....
259                    // ______________     =>        \_________________.    .   seg2
260                    //                               _________________.    .
261                    //                              /                  ....
262                    //                             /       strand2
263                    //                            / seg3     <<<
264                    //                           /
265                    //
266                    // seg1 is the old segment
267
268                    SEC_helix_strand *strand1 = seg1->split(start1, end1, &seg2);
269                    SEC_helix_strand *strand2 = 0;
270                    SEC_segment      *seg3    = 0;
271
272                    if (seg1->get_region()->contains_seq_position(start2)) {
273                        seg3    = seg2;
274                        strand2 = strand1;
275                        strand1 = seg1->split(start2, end2, &seg2);
276                    }
277                    else {
278                        sec_assert(seg2->get_region()->contains_seq_position(start2));
279                        strand2 = seg2->split(start2, end2, &seg3);
280                    }
281
282                    sec_assert(are_adjacent_regions(seg1->get_region(), strand1->get_region()));
283                    sec_assert(are_adjacent_regions(strand1->get_region(), seg2->get_region()));
284                    sec_assert(are_adjacent_regions(seg2->get_region(), strand2->get_region()));
285                    sec_assert(are_adjacent_regions(strand2->get_region(), seg3->get_region()));
286
287                    new_helix = new SEC_helix(this, strand2, strand1); // strands are responsible for memory
288
289                    strand1->set_next_segment(seg3);
290                    strand2->set_next_segment(seg2);
291
292                    new_loop = new SEC_loop(this);
293
294                    seg2->set_loop(new_loop);
295
296                    strand2->set_origin_loop(new_loop);
297                    new_loop->set_fixpoint_strand(strand2);
298                }
299                else { // split two segments
300                    //                                   \                       /
301                    //     seg1                           \ seg1     >>>        /
302                    // ______________                      \       strand1     / seg3
303                    //                                      \_________________/
304                    //                    =>     old_loop    _________________        new_loop
305                    // ______________                       /                 \                                     .
306                    //     seg2                            /       strand2     \ seg2
307                    //                                    / seg4     <<<        \                                   .
308                    //                                   /                       \                                  .
309                    //
310                    // seg1 and seg2 are the old segments
311
312                    // maybe swap seg1/seg2 (to ensure fixpoint-strand stays in old loop)
313                    for (SEC_segment *s = seg1; s != seg2;) {
314                        SEC_helix_strand *hs = s->get_next_strand();
315                        if (!hs->isRootsideFixpoint()) { // fixpoint-strand is between seg1 -> seg2
316                            // swap seg1<->seg2
317                            swap(seg1, seg2);
318                            swap(start1, start2);
319                            swap(end1, end2);
320                            break;
321                        }
322                        s = hs->get_next_segment();
323                    }
324
325                    SEC_segment *seg3;
326                    SEC_segment *seg4;
327                    SEC_helix_strand *strand1 = seg1->split(start1, end1, &seg3);
328                    SEC_helix_strand *strand2 = seg2->split(start2, end2, &seg4);
329
330                    new_helix = new SEC_helix(this, strand2, strand1);
331
332                    strand1->set_next_segment(seg4);
333                    strand2->set_next_segment(seg3);
334
335                    new_loop = new SEC_loop(this);
336
337                    for (SEC_segment *s = seg3; ;) {
338                        s->set_loop(new_loop);
339                        SEC_helix_strand *h = s->get_next_strand();
340                        h->set_origin_loop(new_loop);
341                        if (s == seg2) break;
342                        s = h->get_next_segment();
343                    }
344
345                    new_loop->set_fixpoint_strand(strand2);
346                }
347
348                // set angles of new helix and new loop
349                new_helix->set_rel_angle(0); // wrong, but relayout fails otherwise
350                new_loop->set_rel_angle(0);
351
352                relayout();
353
354                // correct angles of other helices
355                oldAngles.restoreAll(new_loop);
356                oldAngles.set_angle(new_helix, Angle(0));
357                oldAngles.restoreAll(old_loop);
358
359                recalc();
360
361                if (setRootTo) {
362                    set_root(setRootTo); // restore root loop
363                }
364            }
365        }
366    }
367    return error;
368}
369
370// ----------------------
371//      fold a strand
372// ----------------------
373
374GB_ERROR SEC_root::unsplit_loop(SEC_helix_strand *remove_strand) {
375    //
376    //          \ before[0]               / after[1]
377    //           \                       /
378    //            \        >>>>         /
379    //             \      strand[0]    /
380    //              \_________________/
381    //   loop[0]     _________________           loop[1]
382    //              /     strand[1]   \                             .
383    //             /       <<<<        \                            .
384    //            /                     \                           .
385    //           /                       \                          .
386    //          / after[0]                \ before[1]
387    //
388    // The strands are removed and segments get connected.
389    // One loop is deleted.
390
391    GB_ERROR error = 0;
392    SEC_helix_strand *strand[2] = { remove_strand, remove_strand->get_other_strand() };
393    SEC_segment *before[2], *after[2];
394    SEC_loop    *loop[2];
395
396#if defined(CHECK_INTEGRITY)
397    check_integrity(CHECK_STRUCTURE);
398#endif // CHECK_INTEGRITY
399
400    int s;
401    for (s = 0; s<2; s++) {
402        after[s]  = strand[s]->get_next_segment();
403        before[s] = strand[s]->get_previous_segment();
404        loop[s]   = strand[s]->get_origin_loop();
405
406        sec_assert(before[s]->get_loop() == loop[s]);
407        sec_assert(after[s]->get_loop() == loop[s]);
408    }
409
410    bool is_terminal_loop[2] = { before[0] == after[0], before[1] == after[1] };
411    int  i0      = -1;          // index of terminal loop (or -1)
412    bool unsplit = true;
413
414    if (is_terminal_loop[0]) {
415        if (is_terminal_loop[1]) {
416            error   = "You cannot delete the last helix";
417            unsplit = false;
418        }
419        else i0 = 0;
420    }
421    else {
422        if (is_terminal_loop[1]) i0 = 1;
423    }
424
425    if (unsplit) {
426        int del = i0 >= 0 ? i0 : 1; // index of loop which will be deleted
427
428        SEC_loop *setRootTo = 0; // set root back to afterwards?
429
430        {
431            // move away root-loop to make things easy
432            SEC_loop *rootLoop = get_root_loop();
433            if (loop[0] == rootLoop || loop[1] == rootLoop) {
434                SEC_loop         *termLoop    = is_terminal_loop[0] ? loop[0] : loop[1];
435                SEC_helix_strand *toTerm      = strand[0]->get_helix()->strandTowards(termLoop);
436                SEC_helix_strand *toNextLoop  = toTerm->get_next_segment()->get_next_strand();
437                SEC_loop         *anotherLoop = toNextLoop->get_destination_loop();
438
439                sec_assert(anotherLoop != loop[0] && anotherLoop != loop[1]);
440
441                set_root(anotherLoop);
442                setRootTo = loop[1-del]; // afterwards set root back to non-deleted loop
443            }
444        }
445
446        SEC_helix   *removed_helix = strand[0]->get_helix();
447        AngleBuffer  oldAngles(BUFFER_CENTER_RELATIVE);
448        oldAngles.storeAllHelices(loop[0], removed_helix);
449        oldAngles.storeAllHelices(loop[1], removed_helix);
450
451        if (i0 >= 0) { // one loop is terminal
452            // i0 and i1 are indexes 0 and 1 in picture above.
453            // The left loop (loop[i0]) will be removed.
454            int i1 = 1-i0; // index of non-terminal loop
455
456            before[i1]->mergeWith(after[i1], loop[i1]);
457
458            sec_assert(after[i0] == before[i0]);
459            delete after[i0];       // delete the segment of the terminal-loop
460        }
461        else { // none of the loops is terminal
462            // keep loop[0], delete loop[1]
463
464            SEC_helix_strand *rootsideStrand = strand[0]->get_helix()->rootsideLoop()->get_rootside_strand();
465
466            before[1]->mergeWith(after[0], loop[0]);
467            before[0]->mergeWith(after[1], loop[0]);
468            // after[] segments are invalid now!
469
470            // loop over all segments in loop[1] and relink them to loop[0]
471            SEC_segment *seg = before[0];
472            while (seg != before[1]) {
473                SEC_helix_strand *loop1_strand = seg->get_next_strand();
474                loop1_strand->set_origin_loop(loop[0]);
475
476                seg = loop1_strand->get_next_segment();
477                seg->set_loop(loop[0]);
478            }
479
480            loop[0]->set_fixpoint_strand(rootsideStrand);
481        }
482
483        loop[del]->set_fixpoint_strand(NULL);
484        delete loop[del];
485
486        for (s = 0; s<2; s++) strand[s]->unlink(false);
487        delete strand[0];       // delete both strands
488
489        relayout();
490
491        oldAngles.restoreAll(loop[1-del]);
492        recalc();
493
494        if (setRootTo) set_root(setRootTo);
495    }
496
497    return error;
498}
Note: See TracBrowser for help on using the repository browser.