Context Navigation

source: tags/cvs_2_svn/SECEDIT/SEC_split.cxx

Visit:

Last change on this file was 4803, checked in by westram, 17 years ago
rewrote SECEDIT
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: