7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 #include "dp_lib_header.h"
11 int check_link (CL_node ***G, int S1, int r1, int s2, int r2);
12 void light_nodes (CL_node *A, int va, CL_node*B, int vb, CL_node*C,int vc, char *s);
13 void print_graph (CL_node *G, Sequence *S);
20 Alignment * add_constraint2aln ( Alignment *A, int s1, int r1, int s2, int r2)
22 /*Note sCL_node ***G;1 and r1 must be numbered from 0 to n-1*/
27 G=add_constraint2graph_aln (G,s1,r1, s2, r2);
28 A=graph2aln (A,G[0][0],aln2seq(A));
29 vfree_graph (G[0][0]);
32 Alignment * graph_aln (Alignment *A, Constraint_list *iCL, Sequence *S)
36 CLIST_TYPE *entry=NULL;
40 CL=duplicate_constraint_list (iCL);
41 for ( a=0; a<CL->ne; a++)
43 CL->L[a*CL->entry_len+WE]=CL->evaluate_residue_pair ( iCL, CL->L[a*CL->entry_len+SEQ1],CL->L[a+CL->entry_len+R1] ,CL->L[a*CL->entry_len+SEQ2],CL->L[a*CL->entry_len+R2]);
45 CL=sort_constraint_list_on_n_fields(CL, 0, CL->ne,WE,1);
51 for (a=start; a<CL->ne; a++)
55 entry=extract_entry(entry,a, CL);
56 fprintf ( stderr, "\r\tCompletion: [%5.2f%%]",(float)(a*100)/CL->ne);
57 G=add_constraint2graph_aln (G, entry[SEQ1], entry[R1]-1, entry[SEQ2],entry[R2]-1);
60 for (a=start; a<CL->ne; a++)
64 entry=extract_entry(entry,a, CL);
65 fprintf ( stderr, "\r\tCompletion: [%5.2f%%]",(float)(a*100)/CL->ne);
66 G=add_constraint2graph_aln (G, entry[SEQ1], entry[R1]-1, entry[SEQ2],entry[R2]-1);
69 A=graph2aln (A,G[0][0], S);
75 void print_graph (CL_node *G, Sequence *S)
80 A=seq2aln (Seq, A, 1);
81 A=graph2aln (A,G, Seq);
85 Alignment* graph2aln (Alignment *A, CL_node *G, Sequence *S)
96 while (Gi){Gi=Gi->r;l++;}
98 while (Gi){Gi=Gi->c;s++;}
100 A=realloc_alignment (A, l+1);
111 if ( G->res==-1)A->seq_al[s][l]='-';
112 else if ( G->res==-2)A->seq_al[s][l]='*';
113 else if ( G->res==-3)A->seq_al[s][l]='#';
114 else if (G->res>=0)A->seq_al[s][l]=S->seq[G->seq][G->res];
126 A->seq_al[a][l]='\0';
127 A->len_aln=strlen (A->seq_al[0]);
133 CL_node ***aln2graph (Alignment *A)
136 static CL_node ***galn;
137 CL_node *N, *iN, *pN;
142 galn=calloc ( A->nseq, sizeof (CL_node**));
143 for ( a=0; a< A->nseq; a++)
145 galn[a]=calloc (A->len_aln, sizeof (CL_node*));
149 N=declare_cl_nodes(-1, a);
152 for ( a=0; a<A->nseq; a++)
155 for (res=A->order[a][1], b=0; b<A->len_aln; b++)
159 N->r=declare_cl_nodes(-1, a);
168 N->seq=A->order[a][0];
169 N->res=(is_gap(A->seq_al[a][b]))?-1:res++;
171 if (N->res!=-1)galn[a][res-1]=N;
180 if ( a<A->nseq-1)iN->c=declare_cl_nodes(-1, a);
187 CL_node ***add_constraint2graph_aln (CL_node ***G, int s1, int r1, int s2, int r2)
197 d=get_node_distance (S,E);
198 if (d<0){B=S;S=E;E=B;}
202 insert_gap_columns (E,d);
203 shift_segment( S,d+1,d);
210 CL_node *shift_segment ( CL_node *S, int segL, int shiftL)
216 if ( !shiftL)return S;
218 /*find segment coordinates*/
219 for (E=S, a=1; a< segL; a++)E=E->r;
223 G=swap_gap_in_graph (S, E);
224 for (a=1; a< shiftL; a++)swap_gap_in_graph (S, E);
226 while (G!=E)G=remove_graph_gap_column (G);
227 remove_graph_gap_column (E);
232 int is_graph_gap_column(CL_node *S)
238 if (S->res>=0)return 0;
243 CL_node * remove_graph_gap_column (CL_node *S)
245 CL_node *R,*L, *P, *RV;
254 if ( !is_graph_gap_column (S))return RV;
275 CL_node * swap_gap_in_graph ( CL_node*S, CL_node *E)
277 /*Moves gap AFTER End to BEFORE Start
280 straightens the links in between
282 CL_node *G, *N, *iE, *iS, *SP, *SC, *SL;
285 /*Preserve the E/S values*/
290 /*prepare the parent/child links first*/
301 if (N->p)(S->p)->c=S;
304 if (N->c)(S->c)->p=S;
314 if ( G->res>=0)fprintf ( stderr, "\nERROR: NOT a GAP");
317 if (E->r)(E->r)->l=E;
338 CL_node * declare_cl_nodes ( int len, int seq)
347 IN=calloc ( 1, sizeof (CL_node));
357 N=calloc (len, sizeof (CL_node*));
360 if ( len==0)return NULL;
362 for (a=0; a<len; a++)N[a]=calloc ( 1, sizeof (CL_node));
363 for (a=0; a<len; a++)
367 if (a!=0)(N[a])->l=N[a-1];
368 if (a!=len-1)(N[a])->r=N[a+1];
377 CL_node *insert_gap_columns (CL_node *S, int d)
379 CL_node *Gs,*Ge, *pGs, *Gi, *Si;
386 while (S->p!=NULL)S=S->p;
390 Gs=declare_cl_nodes(d, S->seq);
394 if (Ge->r)(Ge->r)->l=Ge;
420 int get_node_distance ( CL_node *S, CL_node *E)
426 /*project the two points onto one sequence*/
427 if (S->seq>E->seq){B=S;S=E;E=B;swap*=-1;}
428 while (S->seq!=E->seq)S=S->c;
430 /*Walk from E to S */
432 while ( iS->res<0 && iS->r!=NULL){iS=iS->r;}
433 if (iS->res<0 || iS->res>E->res){B=S; S=E; E=B;swap*=-1;}
451 int check_graph ( CL_node *S, char *string)
457 if ( S==NULL)S=Start;
458 fprintf ( stderr, "\n\tGRAPH Check %s #%d\n",string, ++n);
459 while ( S->p!=NULL)S=S->p;
460 while ( S->l!=NULL)S=S->l;
467 if (iS->l && (iS->l)->seq!=iS->seq){fprintf ( stderr, "\n\t\tSEq pb");myexit(EXIT_FAILURE);}
468 if (iS->free==1){fprintf ( stderr, "\n\t\tFree Node read");myexit(EXIT_FAILURE);}
471 if (lr!=-1 && iS->res-lr!=1){fprintf ( stderr, "\n\t\tERROR: lost residues");myexit (EXIT_FAILURE);}
474 if ( iS->r && (iS->r)->l!=iS){fprintf ( stderr, "\n\t\tERROR: left != right: [%d %d][%d %d]", iS->seq, iS->res, (iS->l)->seq, (iS->r)->res);myexit (EXIT_FAILURE);}
475 if ( iS->p && (iS->p)->c!=iS){fprintf ( stderr, "\n\t\tERROR: parent != child: [%d %d][%d %d]", iS->seq, iS->res, (iS->p)->seq, (iS->p)->res);myexit (EXIT_FAILURE);}
476 if ( iS->c && (iS->c)->p!=iS){fprintf ( stderr, "\n\t\tERROR: parent != child: [%d %d][%d %d]", iS->seq, iS->res, (iS->c)->seq, (iS->c)->res);myexit (EXIT_FAILURE);}
484 CL_node * vfree_graph ( CL_node *S)
488 while ( S->p!=NULL)S=S->p;
489 while ( S->l!=NULL)S=S->l;
498 if (S)vfree_cl_node (S->l);
505 CL_node *vfree_cl_node ( CL_node *N)
507 if ( N->free==1)crash("freeing free block");
514 void light_nodes (CL_node *A, int va, CL_node*B, int vb, CL_node*C,int vc, char *string )
516 int ta=0, tb=0, tc=0;
518 fprintf ( stderr, "\nCycle %d\n LIGHT NODE: %s", cycle,string);
519 if ( A){ta=A->res; A->res=va;fprintf ( stderr, "\nA: seq %d res %d", A->seq, A->res);}
520 if ( B){tb=B->res; B->res=vb;fprintf ( stderr, "\nB: seq %d res %d", B->seq, B->res);}
521 if ( C){tc=C->res; C->res=vc;fprintf ( stderr, "\nC: seq %d res %d", C->seq, C->res);}
527 int check_link (CL_node ***G, int s1, int r1, int s2, int r2)
542 /*********************************COPYRIGHT NOTICE**********************************/
543 /*© Centro de Regulacio Genomica */
545 /*Cedric Notredame */
546 /*Tue Oct 27 10:12:26 WEST 2009. */
547 /*All rights reserved.*/
548 /*This file is part of T-COFFEE.*/
550 /* T-COFFEE is free software; you can redistribute it and/or modify*/
551 /* it under the terms of the GNU General Public License as published by*/
552 /* the Free Software Foundation; either version 2 of the License, or*/
553 /* (at your option) any later version.*/
555 /* T-COFFEE is distributed in the hope that it will be useful,*/
556 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
557 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
558 /* GNU General Public License for more details.*/
560 /* You should have received a copy of the GNU General Public License*/
561 /* along with Foobar; if not, write to the Free Software*/
562 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
563 /*............................................... |*/
564 /* If you need some more information*/
565 /* cedric.notredame@europe.com*/
566 /*............................................... |*/
570 /*********************************COPYRIGHT NOTICE**********************************/