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
  40   CL=duplicate_constraint_list (iCL);
  41   for ( a=0; a<CL->ne; a++)
  42     {
  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]);
  44      }
  45    CL=sort_constraint_list_on_n_fields(CL, 0, CL->ne,WE,1);
  46    G=aln2graph (A);
  47
  48
  49   Start=G[0][0];
  50   start=0;
  51   for (a=start; a<CL->ne; a++)
  52     {
  53
  54       cycle++;
  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);
  58     }
  59    start=0;
  60   for (a=start; a<CL->ne; a++)
  61     {
  62
  63       cycle++;
  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);
  67     }
  68
  69   A=graph2aln (A,G[0][0], S);
  70
  71   return A;
  72
  73 }
  74
  75 void print_graph (CL_node *G, Sequence *S)
  76 {
  77   Alignment *A=NULL;
  78
  79   if (S==NULL) S=Seq;
  80   A=seq2aln (Seq, A, 1);
  81   A=graph2aln (A,G, Seq);
  82   print_aln (A);
  83 }
  84
  85 Alignment* graph2aln (Alignment *A, CL_node *G, Sequence *S)
  86 {
  87   int s, l, a;
  88   CL_node *Gi;
  89
  90   /*Rewind G*/
  91   while ( G->p)G=G->p;
  92   while ( G->l)G=G->l;
  93
  94   l=s=0;
  95   Gi=G;
  96   while (Gi){Gi=Gi->r;l++;}
  97   Gi=G;
  98   while (Gi){Gi=Gi->c;s++;}
  99
 100   A=realloc_alignment (A, l+1);
 101
 102
 103   l=0;
 104   while (G)
 105     {
 106       Gi=G;
 107       s=0;
 108       while ( G!=NULL)
 109         {
 110
 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];
 115
 116           G=G->c;
 117           s++;
 118         }
 119       G=Gi->r;
 120       l++;
 121     }
 122
 123
 124
 125   for ( a=0;a<s; a++)
 126     A->seq_al[a][l]='\0';
 127   A->len_aln=strlen (A->seq_al[0]);
 128   A->nseq=s;
 129
 130
 131   return A;
 132 }
 133 CL_node ***aln2graph (Alignment *A)
 134 {
 135   int a=0, b;
 136   static CL_node ***galn;
 137   CL_node *N, *iN, *pN;
 138   int res;
 139
 140   if ( !galn)
 141     {
 142       galn=calloc ( A->nseq, sizeof (CL_node**));
 143       for ( a=0; a< A->nseq; a++)
 144         {
 145           galn[a]=calloc (A->len_aln, sizeof (CL_node*));
 146         }
 147     }
 148   pN=iN=NULL;
 149   N=declare_cl_nodes(-1, a);
 150
 151
 152   for ( a=0; a<A->nseq; a++)
 153     {
 154       iN=N;
 155       for (res=A->order[a][1], b=0; b<A->len_aln; b++)
 156         {
 157         if (b<A->len_aln-1)
 158           {
 159             N->r=declare_cl_nodes(-1, a);
 160             (N->r)->l=N;
 161           }
 162         if ( pN)
 163           {
 164             N->p=pN;
 165             pN->c=N;
 166           }
 167
 168         N->seq=A->order[a][0];
 169         N->res=(is_gap(A->seq_al[a][b]))?-1:res++;
 170
 171         if (N->res!=-1)galn[a][res-1]=N;
 172         if ( pN)
 173           {
 174             N->p=pN;
 175             pN->c=N;
 176             pN=pN->r;
 177           }
 178         N=N->r;
 179       }
 180       if ( a<A->nseq-1)iN->c=declare_cl_nodes(-1, a);
 181       pN=iN;
 182       N=iN->c;
 183     }
 184
 185   return galn;
 186 }
 187 CL_node ***add_constraint2graph_aln (CL_node ***G, int s1, int r1, int s2, int r2)
 188 {
 189   CL_node *S, *E, *B;
 190   int d;
 191
 192
 193   S=G[s1][r1];
 194   E=G[s2][r2];
 195
 196
 197   d=get_node_distance (S,E);
 198   if (d<0){B=S;S=E;E=B;}
 199   d=(d<0)?-d:d;
 200
 201
 202   insert_gap_columns (E,d);
 203   shift_segment( S,d+1,d);
 204
 205   return G;
 206
 207 }
 208
 209
 210 CL_node *shift_segment ( CL_node *S, int segL, int shiftL)
 211 {
 212   int a;
 213   CL_node *E, *G;
 214
 215
 216   if ( !shiftL)return S;
 217
 218   /*find segment coordinates*/
 219   for (E=S, a=1; a< segL; a++)E=E->r;
 220
 221   /*Shift the gaps*/
 222
 223   G=swap_gap_in_graph (S, E);
 224   for (a=1; a< shiftL; a++)swap_gap_in_graph (S, E);
 225
 226   while (G!=E)G=remove_graph_gap_column (G);
 227   remove_graph_gap_column (E);
 228
 229   return G;
 230 }
 231
 232 int  is_graph_gap_column(CL_node *S)
 233 {
 234   while (S->p)S=S->p;
 235
 236   while (S)
 237     {
 238       if (S->res>=0)return 0;
 239       S=S->c;
 240     }
 241   return 1;
 242 }
 243 CL_node * remove_graph_gap_column (CL_node *S)
 244 {
 245   CL_node *R,*L, *P, *RV;
 246
 247   RV=S->r;
 248   while (S->p)
 249     {
 250
 251       S=S->p;
 252     }
 253
 254   if ( !is_graph_gap_column (S))return RV;
 255
 256
 257
 258   while (S)
 259     {
 260
 261       R=S->r;
 262       L=S->l;
 263       P=S->p;
 264
 265       if (L)L->r=S->r;
 266       if (R)R->l=S->l;
 267
 268       P=S;
 269       S=S->c;
 270       vfree_cl_node (P);
 271     }
 272   return RV;
 273 }
 274
 275 CL_node * swap_gap_in_graph ( CL_node*S, CL_node *E)
 276 {
 277   /*Moves gap AFTER End to BEFORE Start
 278     SxxxE-
 279     -xxxxx
 280    straightens the links in between
 281   */
 282   CL_node *G, *N, *iE, *iS, *SP, *SC, *SL;
 283
 284
 285   /*Preserve the E/S values*/
 286   iE=E;
 287   iS=S;
 288
 289
 290   /*prepare the parent/child links first*/
 291
 292   SP=S->p;
 293   SC=S->c;
 294   SL=S->l;
 295
 296   while ( S!=E->r)
 297     {
 298       N=S->r;
 299
 300       S->p=N->p;
 301       if (N->p)(S->p)->c=S;
 302
 303       S->c=N->c;
 304       if (N->c)(S->c)->p=S;
 305
 306       S=S->r;
 307     }
 308
 309   E=iE;
 310   S=iS;
 311
 312   /*Remove the gap*/
 313   G=E->r;
 314   if ( G->res>=0)fprintf ( stderr, "\nERROR: NOT a GAP");
 315
 316   E->r=G->r;
 317   if (E->r)(E->r)->l=E;
 318
 319   /*insert the gap*/
 320
 321   G->r=S;
 322   S->l=G;
 323
 324   G->l=SL;
 325   if (SL)SL->r=G;
 326
 327
 328   G->p=SP;
 329   if (SP)SP->c=G;
 330
 331   G->c=SC;
 332   if (SC)SC->p=G;
 333
 334   return G;
 335
 336 }
 337
 338 CL_node * declare_cl_nodes ( int len, int seq)
 339 {
 340   static CL_node **N;
 341   CL_node *IN;
 342   static int Nlen;
 343   int a;
 344
 345   if (len==-1)
 346     {
 347       IN=calloc ( 1, sizeof (CL_node));
 348       IN->res=-1;
 349       return IN;
 350     }
 351
 352
 353
 354   if ( len>Nlen)
 355     {
 356       free (N);
 357       N=calloc (len, sizeof (CL_node*));
 358     }
 359
 360   if ( len==0)return NULL;
 361
 362   for (a=0; a<len; a++)N[a]=calloc ( 1, sizeof (CL_node));
 363   for (a=0; a<len; a++)
 364     {
 365       (N[a])->res=-1;
 366       (N[a])->seq=seq;
 367       if (a!=0)(N[a])->l=N[a-1];
 368       if (a!=len-1)(N[a])->r=N[a+1];
 369     }
 370
 371   (N[0])->l=N[len-1];
 372   (N[len-1])->r=N[0];
 373
 374   return N[0];
 375 }
 376
 377 CL_node *insert_gap_columns (CL_node *S, int d)
 378 {
 379   CL_node *Gs,*Ge, *pGs, *Gi, *Si;
 380   int a;
 381
 382   if ( d==0)return S;
 383
 384   pGs=Gi=NULL;
 385   Si=S;
 386   while (S->p!=NULL)S=S->p;
 387
 388   while (S!=NULL)
 389     {
 390       Gs=declare_cl_nodes(d, S->seq);
 391       Ge=Gs->l;
 392
 393       Ge->r=S->r;
 394       if (Ge->r)(Ge->r)->l=Ge;
 395
 396       Gs->l=S;
 397       S->r=Gs;
 398
 399       if (pGs)
 400         {
 401           Gi=Gs;
 402           for (a=0; a< d; a++)
 403             {
 404               Gs->p=pGs;
 405               pGs->c=Gs;
 406               Gs=Gs->r;
 407               pGs=pGs->r;
 408             }
 409           pGs=Gi;
 410         }
 411       else
 412         {
 413           pGs=Gs;
 414         }
 415       S=S->c;
 416     }
 417   return Si;
 418 }
 419
 420 int get_node_distance ( CL_node *S, CL_node *E)
 421 {
 422   int distance=0;
 423   CL_node *iS,*B;
 424   int swap=1;
 425
 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;
 429
 430   /*Walk from E to S */
 431   iS=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;}
 434
 435   while ( S!=E)
 436     {
 437       S=S->r;
 438       distance+=swap;
 439     }
 440   return distance;
 441 }
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451 int check_graph ( CL_node *S, char *string)
 452 {
 453   CL_node *iS;
 454   static int n;
 455   int lr;
 456
 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;
 461   while ( S)
 462     {
 463       iS=S;
 464       lr=-1;
 465       while (iS)
 466         {
 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);}
 469           if (iS->res>0)
 470             {
 471               if (lr!=-1 && iS->res-lr!=1){fprintf ( stderr, "\n\t\tERROR: lost residues");myexit (EXIT_FAILURE);}
 472               lr=iS->res;
 473             }
 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);}
 477           iS=iS->r;
 478         }
 479       S=S->c;
 480     }
 481   return 1;
 482 }
 483
 484 CL_node * vfree_graph ( CL_node *S)
 485 {
 486   CL_node *Si;
 487
 488   while ( S->p!=NULL)S=S->p;
 489   while ( S->l!=NULL)S=S->l;
 490
 491   while ( S)
 492     {
 493       Si=S->c;
 494       while ( S)
 495         {
 496
 497           S=S->r;
 498           if (S)vfree_cl_node (S->l);
 499         }
 500       S=Si;
 501     }
 502   return S;
 503
 504 }
 505 CL_node *vfree_cl_node ( CL_node *N)
 506 {
 507   if ( N->free==1)crash("freeing free block");
 508   N->free=1;
 509   free (N);
 510   return N;
 511 }
 512
 513
 514 void light_nodes (CL_node *A, int va, CL_node*B, int vb, CL_node*C,int vc, char *string )
 515 {
 516   int ta=0, tb=0, tc=0;
 517
 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);}
 522   print_graph (A, 0);
 523   if ( A){A->res=ta;}
 524   if ( B){B->res=tb;}
 525   if ( C){C->res=tc;}
 526 }
 527 int check_link (CL_node ***G, int s1, int r1, int s2, int r2)
 528 {
 529   CL_node *S;
 530   CL_node *E;
 531
 532   S=G[s1][r1];
 533   E=G[s2][r2];
 534   while ( S->p)S=S->p;
 535   while ( S)
 536     {
 537       S=S->c;
 538       if ( S==E)return 1;
 539     }
 540   return 0;
 541 }
 542 /*********************************COPYRIGHT NOTICE**********************************/
 543 /*© Centro de Regulacio Genomica */
 544 /*and */
 545 /*Cedric Notredame */
 546 /*Tue Oct 27 10:12:26 WEST 2009. */
 547 /*All rights reserved.*/
 548 /*This file is part of T-COFFEE.*/
 549 /**/
 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.*/
 554 /**/
 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.*/
 559 /**/
 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 /*...............................................                                                                                                                                     |*/
 567 /**/
 568 /**/
 569 /*      */
 570 /*********************************COPYRIGHT NOTICE**********************************/