Context Navigation

source: tags/arb-6.0/SECEDIT/SEC_split.cxx

Visit:

Last change on this file was 8725, checked in by westram, 12 years ago
comments cleanup
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.7 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
17	using namespace std;
18
19	enum AngleBufferMode {
20	BUFFER_ABSOLUTE_ANGLES,
21	BUFFER_CENTER_RELATIVE,
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) {
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	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
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
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
138	SEC_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 0;
163	}
164
165	inline SEC_segment findSegmentContaining(SEC_root root, int pos, GB_ERROR& error) {
166	SEC_segment *result = 0;
167	error = 0;
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 = 0;
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 = 0;
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 = 0;
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 = 0;
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 = 0; // 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, 0);
248
249	SEC_helix *new_helix = 0;
250	SEC_loop *new_loop = 0;
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 = 0;
267	SEC_segment *seg3 = 0;
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 = 0;
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 = 0; // 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(NULL);
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.

Download in other formats: