binaries/src/tcoffee/t_coffee_source/util_graph_maln.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6
   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);
  14 Sequence *Seq;
  15 CL_node *Start;
  16 static int cycle;
  17
  18
  19
  20 Alignment * add_constraint2aln ( Alignment *A, int s1, int r1, int s2, int r2)
  21 {
  22   /*Note sCL_node ***G;1 and r1 must be numbered from 0 to n-1*/
  23   CL_node ***G;
  24
  25
  26   G=aln2graph(A);
  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]);
  30   return A;
  31 }
  32 Alignment * graph_aln (Alignment *A, Constraint_list *iCL, Sequence *S)
  33 {
  34   CL_node ***G;
  35   int a,start;
  36   CLIST_TYPE *entry=NULL;
  37   Constraint_list *CL;
  38   Seq=S;
  39   HERE ("Function graph_aln is deprecated");
  40   return A;
  41
  42 }
  43
  44 void print_graph (CL_node *G, Sequence *S)
  45 {
  46   Alignment *A=NULL;
  47
  48   if (S==NULL) S=Seq;
  49   A=seq2aln (Seq, A, 1);
  50   A=graph2aln (A,G, Seq);
  51   print_aln (A);
  52 }
  53
  54 Alignment* graph2aln (Alignment *A, CL_node *G, Sequence *S)
  55 {
  56   int s, l, a;
  57   CL_node *Gi;
  58
  59   /*Rewind G*/
  60   while ( G->p)G=G->p;
  61   while ( G->l)G=G->l;
  62
  63   l=s=0;
  64   Gi=G;
  65   while (Gi){Gi=Gi->r;l++;}
  66   Gi=G;
  67   while (Gi){Gi=Gi->c;s++;}
  68
  69   A=realloc_alignment (A, l+1);
  70
  71
  72   l=0;
  73   while (G)
  74     {
  75       Gi=G;
  76       s=0;
  77       while ( G!=NULL)
  78         {
  79
  80           if ( G->res==-1)A->seq_al[s][l]='-';
  81           else if ( G->res==-2)A->seq_al[s][l]='*';
  82           else if ( G->res==-3)A->seq_al[s][l]='#';
  83           else if (G->res>=0)A->seq_al[s][l]=S->seq[G->seq][G->res];
  84
  85           G=G->c;
  86           s++;
  87         }
  88       G=Gi->r;
  89       l++;
  90     }
  91
  92
  93
  94   for ( a=0;a<s; a++)
  95     A->seq_al[a][l]='\0';
  96   A->len_aln=strlen (A->seq_al[0]);
  97   A->nseq=s;
  98
  99
 100   return A;
 101 }
 102 CL_node ***aln2graph (Alignment *A)
 103 {
 104   int a=0, b;
 105   static CL_node ***galn;
 106   CL_node *N, *iN, *pN;
 107   int res;
 108
 109   if ( !galn)
 110     {
 111       galn=calloc ( A->nseq, sizeof (CL_node**));
 112       for ( a=0; a< A->nseq; a++)
 113         {
 114           galn[a]=calloc (A->len_aln, sizeof (CL_node*));
 115         }
 116     }
 117   pN=iN=NULL;
 118   N=declare_cl_nodes(-1, a);
 119
 120
 121   for ( a=0; a<A->nseq; a++)
 122     {
 123       iN=N;
 124       for (res=A->order[a][1], b=0; b<A->len_aln; b++)
 125         {
 126         if (b<A->len_aln-1)
 127           {
 128             N->r=declare_cl_nodes(-1, a);
 129             (N->r)->l=N;
 130           }
 131         if ( pN)
 132           {
 133             N->p=pN;
 134             pN->c=N;
 135           }
 136
 137         N->seq=A->order[a][0];
 138         N->res=(is_gap(A->seq_al[a][b]))?-1:res++;
 139
 140         if (N->res!=-1)galn[a][res-1]=N;
 141         if ( pN)
 142           {
 143             N->p=pN;
 144             pN->c=N;
 145             pN=pN->r;
 146           }
 147         N=N->r;
 148       }
 149       if ( a<A->nseq-1)iN->c=declare_cl_nodes(-1, a);
 150       pN=iN;
 151       N=iN->c;
 152     }
 153
 154   return galn;
 155 }
 156 CL_node ***add_constraint2graph_aln (CL_node ***G, int s1, int r1, int s2, int r2)
 157 {
 158   CL_node *S, *E, *B;
 159   int d;
 160
 161
 162   S=G[s1][r1];
 163   E=G[s2][r2];
 164
 165
 166   d=get_node_distance (S,E);
 167   if (d<0){B=S;S=E;E=B;}
 168   d=(d<0)?-d:d;
 169
 170
 171   insert_gap_columns (E,d);
 172   shift_segment( S,d+1,d);
 173
 174   return G;
 175
 176 }
 177
 178
 179 CL_node *shift_segment ( CL_node *S, int segL, int shiftL)
 180 {
 181   int a;
 182   CL_node *E, *G;
 183
 184
 185   if ( !shiftL)return S;
 186
 187   /*find segment coordinates*/
 188   for (E=S, a=1; a< segL; a++)E=E->r;
 189
 190   /*Shift the gaps*/
 191
 192   G=swap_gap_in_graph (S, E);
 193   for (a=1; a< shiftL; a++)swap_gap_in_graph (S, E);
 194
 195   while (G!=E)G=remove_graph_gap_column (G);
 196   remove_graph_gap_column (E);
 197
 198   return G;
 199 }
 200
 201 int  is_graph_gap_column(CL_node *S)
 202 {
 203   while (S->p)S=S->p;
 204
 205   while (S)
 206     {
 207       if (S->res>=0)return 0;
 208       S=S->c;
 209     }
 210   return 1;
 211 }
 212 CL_node * remove_graph_gap_column (CL_node *S)
 213 {
 214   CL_node *R,*L, *P, *RV;
 215
 216   RV=S->r;
 217   while (S->p)
 218     {
 219
 220       S=S->p;
 221     }
 222
 223   if ( !is_graph_gap_column (S))return RV;
 224
 225
 226
 227   while (S)
 228     {
 229
 230       R=S->r;
 231       L=S->l;
 232       P=S->p;
 233
 234       if (L)L->r=S->r;
 235       if (R)R->l=S->l;
 236
 237       P=S;
 238       S=S->c;
 239       vfree_cl_node (P);
 240     }
 241   return RV;
 242 }
 243
 244 CL_node * swap_gap_in_graph ( CL_node*S, CL_node *E)
 245 {
 246   /*Moves gap AFTER End to BEFORE Start
 247     SxxxE-
 248     -xxxxx
 249    straightens the links in between
 250   */
 251   CL_node *G, *N, *iE, *iS, *SP, *SC, *SL;
 252
 253
 254   /*Preserve the E/S values*/
 255   iE=E;
 256   iS=S;
 257
 258
 259   /*prepare the parent/child links first*/
 260
 261   SP=S->p;
 262   SC=S->c;
 263   SL=S->l;
 264
 265   while ( S!=E->r)
 266     {
 267       N=S->r;
 268
 269       S->p=N->p;
 270       if (N->p)(S->p)->c=S;
 271
 272       S->c=N->c;
 273       if (N->c)(S->c)->p=S;
 274
 275       S=S->r;
 276     }
 277
 278   E=iE;
 279   S=iS;
 280
 281   /*Remove the gap*/
 282   G=E->r;
 283   if ( G->res>=0)fprintf ( stderr, "\nERROR: NOT a GAP");
 284
 285   E->r=G->r;
 286   if (E->r)(E->r)->l=E;
 287
 288   /*insert the gap*/
 289
 290   G->r=S;
 291   S->l=G;
 292
 293   G->l=SL;
 294   if (SL)SL->r=G;
 295
 296
 297   G->p=SP;
 298   if (SP)SP->c=G;
 299
 300   G->c=SC;
 301   if (SC)SC->p=G;
 302
 303   return G;
 304
 305 }
 306
 307 CL_node * declare_cl_nodes ( int len, int seq)
 308 {
 309   static CL_node **N;
 310   CL_node *IN;
 311   static int Nlen;
 312   int a;
 313
 314   if (len==-1)
 315     {
 316       IN=calloc ( 1, sizeof (CL_node));
 317       IN->res=-1;
 318       return IN;
 319     }
 320
 321
 322
 323   if ( len>Nlen)
 324     {
 325       free (N);
 326       N=calloc (len, sizeof (CL_node*));
 327     }
 328
 329   if ( len==0)return NULL;
 330
 331   for (a=0; a<len; a++)N[a]=calloc ( 1, sizeof (CL_node));
 332   for (a=0; a<len; a++)
 333     {
 334       (N[a])->res=-1;
 335       (N[a])->seq=seq;
 336       if (a!=0)(N[a])->l=N[a-1];
 337       if (a!=len-1)(N[a])->r=N[a+1];
 338     }
 339
 340   (N[0])->l=N[len-1];
 341   (N[len-1])->r=N[0];
 342
 343   return N[0];
 344 }
 345
 346 CL_node *insert_gap_columns (CL_node *S, int d)
 347 {
 348   CL_node *Gs,*Ge, *pGs, *Gi, *Si;
 349   int a;
 350
 351   if ( d==0)return S;
 352
 353   pGs=Gi=NULL;
 354   Si=S;
 355   while (S->p!=NULL)S=S->p;
 356
 357   while (S!=NULL)
 358     {
 359       Gs=declare_cl_nodes(d, S->seq);
 360       Ge=Gs->l;
 361
 362       Ge->r=S->r;
 363       if (Ge->r)(Ge->r)->l=Ge;
 364
 365       Gs->l=S;
 366       S->r=Gs;
 367
 368       if (pGs)
 369         {
 370           Gi=Gs;
 371           for (a=0; a< d; a++)
 372             {
 373               Gs->p=pGs;
 374               pGs->c=Gs;
 375               Gs=Gs->r;
 376               pGs=pGs->r;
 377             }
 378           pGs=Gi;
 379         }
 380       else
 381         {
 382           pGs=Gs;
 383         }
 384       S=S->c;
 385     }
 386   return Si;
 387 }
 388
 389 int get_node_distance ( CL_node *S, CL_node *E)
 390 {
 391   int distance=0;
 392   CL_node *iS,*B;
 393   int swap=1;
 394
 395   /*project the two points onto one sequence*/
 396   if (S->seq>E->seq){B=S;S=E;E=B;swap*=-1;}
 397   while (S->seq!=E->seq)S=S->c;
 398
 399   /*Walk from E to S */
 400   iS=S;
 401   while ( iS->res<0 && iS->r!=NULL){iS=iS->r;}
 402   if (iS->res<0 || iS->res>E->res){B=S; S=E; E=B;swap*=-1;}
 403
 404   while ( S!=E)
 405     {
 406       S=S->r;
 407       distance+=swap;
 408     }
 409   return distance;
 410 }
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420 int check_graph ( CL_node *S, char *string)
 421 {
 422   CL_node *iS;
 423   static int n;
 424   int lr;
 425
 426   if ( S==NULL)S=Start;
 427   fprintf ( stderr, "\n\tGRAPH Check %s #%d\n",string, ++n);
 428   while ( S->p!=NULL)S=S->p;
 429   while ( S->l!=NULL)S=S->l;
 430   while ( S)
 431     {
 432       iS=S;
 433       lr=-1;
 434       while (iS)
 435         {
 436           if (iS->l && (iS->l)->seq!=iS->seq){fprintf ( stderr, "\n\t\tSEq pb");myexit(EXIT_FAILURE);}
 437           if (iS->free==1){fprintf ( stderr, "\n\t\tFree Node read");myexit(EXIT_FAILURE);}
 438           if (iS->res>0)
 439             {
 440               if (lr!=-1 && iS->res-lr!=1){fprintf ( stderr, "\n\t\tERROR: lost residues");myexit (EXIT_FAILURE);}
 441               lr=iS->res;
 442             }
 443           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);}
 444           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);}
 445           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);}
 446           iS=iS->r;
 447         }
 448       S=S->c;
 449     }
 450   return 1;
 451 }
 452
 453 CL_node * vfree_graph ( CL_node *S)
 454 {
 455   CL_node *Si;
 456
 457   while ( S->p!=NULL)S=S->p;
 458   while ( S->l!=NULL)S=S->l;
 459
 460   while ( S)
 461     {
 462       Si=S->c;
 463       while ( S)
 464         {
 465
 466           S=S->r;
 467           if (S)vfree_cl_node (S->l);
 468         }
 469       S=Si;
 470     }
 471   return S;
 472
 473 }
 474 CL_node *vfree_cl_node ( CL_node *N)
 475 {
 476   if ( N->free==1)crash("freeing free block");
 477   N->free=1;
 478   free (N);
 479   return N;
 480 }
 481
 482
 483 void light_nodes (CL_node *A, int va, CL_node*B, int vb, CL_node*C,int vc, char *string )
 484 {
 485   int ta=0, tb=0, tc=0;
 486
 487   fprintf ( stderr, "\nCycle %d\n LIGHT NODE: %s", cycle,string);
 488   if ( A){ta=A->res; A->res=va;fprintf ( stderr, "\nA: seq %d res %d", A->seq, A->res);}
 489   if ( B){tb=B->res; B->res=vb;fprintf ( stderr, "\nB: seq %d res %d", B->seq, B->res);}
 490   if ( C){tc=C->res; C->res=vc;fprintf ( stderr, "\nC: seq %d res %d", C->seq, C->res);}
 491   print_graph (A, 0);
 492   if ( A){A->res=ta;}
 493   if ( B){B->res=tb;}
 494   if ( C){C->res=tc;}
 495 }
 496 int check_link (CL_node ***G, int s1, int r1, int s2, int r2)
 497 {
 498   CL_node *S;
 499   CL_node *E;
 500
 501   S=G[s1][r1];
 502   E=G[s2][r2];
 503   while ( S->p)S=S->p;
 504   while ( S)
 505     {
 506       S=S->c;
 507       if ( S==E)return 1;
 508     }
 509   return 0;
 510 }
 511 /******************************COPYRIGHT NOTICE*******************************/
 512 /*© Centro de Regulacio Genomica */
 513 /*and */
 514 /*Cedric Notredame */
 515 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
 516 /*All rights reserved.*/
 517 /*This file is part of T-COFFEE.*/
 518 /**/
 519 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 520 /*    it under the terms of the GNU General Public License as published by*/
 521 /*    the Free Software Foundation; either version 2 of the License, or*/
 522 /*    (at your option) any later version.*/
 523 /**/
 524 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 525 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 526 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 527 /*    GNU General Public License for more details.*/
 528 /**/
 529 /*    You should have received a copy of the GNU General Public License*/
 530 /*    along with Foobar; if not, write to the Free Software*/
 531 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 532 /*...............................................                                                                                      |*/
 533 /*  If you need some more information*/
 534 /*  cedric.notredame@europe.com*/
 535 /*...............................................                                                                                                                                     |*/
 536 /**/
 537 /**/
 538 /*      */
 539 /******************************COPYRIGHT NOTICE*******************************/