binaries/src/tcoffee/t_coffee_source/util_dp_clean_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
  12 Alignment *clean_maln ( Alignment *A, Alignment *I, int T, int n_it)
  13   {
  14     Alignment *C=NULL;
  15     int a, b;
  16     int in_una, in_aln, in_gap, gap, una, aln;
  17     int Sstart,Rstate, Sstate;
  18     int n_segment=0;
  19     int **segment_list;
  20
  21
  22     add_warning ( stderr, "\nWARNING: -clean_aln is not supported anymore [PROGRAM:%s]\n", PROGRAM);
  23     return A;
  24
  25
  26
  27     /*Initialization*/
  28     a=0;
  29     in_una=a++;in_gap=a++;in_aln=a++; aln=a++;gap=a++;una=a++;
  30     segment_list=declare_int ( A->len_aln*A->nseq, 3);
  31
  32
  33     /*1: Identify the segments*/
  34     C=copy_aln(A, C);
  35     for ( a=0; a< A->nseq; a++)
  36       {
  37         Sstate=in_aln;
  38         for ( b=0; b<A->len_aln; b++)
  39           {
  40             if (is_gap(A->seq_al[a][b]))Rstate=gap;
  41             else if ( I->seq_al[a][b]<=T){Rstate=una;}
  42             else if ( I->seq_al[a][b]==NO_COLOR_RESIDUE)Rstate=una;
  43             else Rstate=aln;
  44
  45             if (Rstate==una)C->seq_al[a][b]='-';
  46
  47             if (Sstate==in_aln)
  48               {
  49                 if ( Rstate==gap)
  50                   {Sstate=in_gap;
  51                   Sstart=b;
  52                   }
  53                 else if ( Rstate==una)
  54                   {
  55                     Sstate=in_una;
  56                     Sstart=b;
  57                   }
  58                 else if ( Rstate==aln)
  59                   Sstate=in_aln;
  60               }
  61             else if ( Sstate==in_gap)
  62               {
  63                 if ( Rstate==gap);
  64                 else if ( Rstate==una)Sstate=in_una;
  65                 else if ( Rstate==aln)Sstate=in_aln;
  66               }
  67             else if ( Sstate==in_una)
  68               {
  69                 if ( Rstate==gap);
  70                 else if ( Rstate==una);
  71                 else if ( Rstate==aln)
  72                   {
  73                     segment_list[n_segment][0]=a;
  74                     segment_list[n_segment][1]=Sstart;
  75                     segment_list[n_segment][2]=b-Sstart;
  76                     Sstate=in_aln;
  77                     n_segment++;
  78                   }
  79               }
  80           }
  81         if (Sstate==in_una)
  82           {
  83            segment_list[n_segment][0]=a;
  84            segment_list[n_segment][1]=Sstart;
  85            segment_list[n_segment][2]=b-Sstart;
  86            Sstate=in_aln;
  87            n_segment++;
  88           }
  89       }
  90
  91     /*2 Realign the segments*/
  92
  93     for ( b=0; b< n_it; b++)
  94       {
  95         for ( a=0; a< n_segment; a++)
  96           {
  97             HERE ("1");
  98             A=realign_segment ( segment_list[a][0], segment_list[a][1], segment_list[a][2], A, C);
  99
 100           }
 101       }
 102     free_aln (C);
 103     free_int ( segment_list, -1);
 104     make_fast_generic_dp_pair_wise (NULL, NULL, NULL, NULL);
 105
 106     return A;
 107   }
 108 Alignment *realign_segment (int seq, int start, int len,Alignment *A, Alignment *C)
 109   {
 110     Alignment *S1=NULL, *S2=NULL, *S3=NULL;
 111     int  *ns, **ls;
 112     int a,b;
 113     static Constraint_list *CL;
 114
 115
 116     /*1 Prepare the Constraint list*/
 117     if ( !CL)
 118       {
 119         CL=vcalloc ( 1, sizeof (Constraint_list));
 120         CL->extend_jit=0;
 121         CL->pw_parameters_set=1;
 122         CL->M=read_matrice ("blosum62mt");
 123         CL->gop=-20;
 124         CL->gep=-1;
 125         CL->evaluate_residue_pair=evaluate_matrix_score;
 126         sprintf ( CL->dp_mode, "myers_miller_pair_wise");
 127       }
 128
 129     S1=copy_aln(A,S1);
 130     S1=extract_aln (S1,0,start);
 131     S2=copy_aln(A,S2);
 132     S2=extract_aln (S2, start, start+len);
 133     S3=copy_aln(A,S3);
 134     S3=extract_aln (S3, start+len,A->len_aln);
 135
 136
 137     /*for (a=0; a<S2->nseq; a++){S2->order[a][1]=0;S2->order[a][0]=a;}*/
 138
 139
 140     ungap ( S2->seq_al[seq]);
 141     CL->S=A->S;/*aln2seq(S2);*/
 142     /*3 Prepare Sequence Presentation*/
 143     ns=vcalloc (2, sizeof (int));
 144     ls=declare_int (2,S2->nseq);
 145
 146     ns[0]=A->nseq-1;
 147     for ( a=0,b=0; a< S2->nseq; a++)if (a!=seq)ls[0][b++]=a;
 148     ns[1]=1;
 149     ls[1][0]=seq;
 150
 151     pair_wise (S2, ns, ls, CL);
 152
 153     A=realloc_aln (A, strlen (S1->seq_al[0])+ strlen (S2->seq_al[0])+ strlen (S3->seq_al[0])+1);
 154     for ( a=0; a< A->nseq; a++)
 155       {
 156         sprintf ( A->seq_al[a], "%s%s%s", S1->seq_al[a], S2->seq_al[a], S3->seq_al[a]);
 157       }
 158
 159     free_aln (S1);
 160     free_aln (S2);
 161     free_aln (S3);
 162     vfree(ns);free_int(ls, -1);
 163
 164     return A;
 165   }
 166 Alignment *realign_segment_old (int seq, int start, int len,Alignment *A, Alignment *C)
 167   {
 168     Alignment *S=NULL;
 169     int  *ns, **ls;
 170     char *sub_seq;
 171     static Dp_Model *M=NULL;
 172     static Constraint_list *CL=NULL;
 173     Dp_Result *R=NULL;
 174     int a,b, c;
 175
 176
 177     /*1 Prepare the Constraint list*/
 178     if ( !CL)
 179       {
 180         CL=vcalloc ( 1, sizeof (Constraint_list));
 181         CL->extend_jit=0;
 182         CL->pw_parameters_set=1;
 183         CL->M=read_matrice ("blosum62mt");
 184         CL->gop=-20;
 185         CL->gep=-1;
 186         CL->evaluate_residue_pair=evaluate_matrix_score;
 187
 188       }
 189     S=copy_aln(C,S);
 190     S=extract_aln (S, start, start+len);
 191     S->len_aln=strlen(S->seq_al[0]);
 192     sub_seq=extract_char  (A->seq_al[seq], start, len);
 193
 194     ungap(sub_seq);
 195
 196     sprintf ( S->seq_al[seq],"%s", sub_seq);
 197     CL->S=aln2seq(S);
 198
 199
 200
 201     /*2 Prepare the Model*/
 202     M=initialize_seg2prf_model((start==0)?2:0,(start+len==A->len_aln)?2:0,CL);
 203     M->diag=vcalloc ( 2*len+1, sizeof (int));
 204     M->diag[0]=len+strlen (sub_seq)-1;
 205     for ( a=1; a<=M->diag[0]; a++)M->diag[a]=a;
 206
 207     /*3 Prepare Sequence Presentation*/
 208     ns=vcalloc (2, sizeof (int));
 209     ls=declare_int (2,A->nseq);
 210
 211     ns[0]=A->nseq-1;
 212     for ( a=0,b=0; a< A->nseq; a++)if (a!=seq)ls[0][b++]=a;
 213     ns[1]=1;
 214     ls[1][0]=seq;
 215
 216     if ( strlen (sub_seq)!=len)
 217       {
 218
 219
 220         R=make_fast_generic_dp_pair_wise(S, ns, ls, M);
 221
 222         for (c=0, b=1,a=start; a< start+len; b++,a++)
 223           {
 224             if (R->traceback[b]==0)
 225               {
 226                 A->seq_al[seq][a]=sub_seq[c];
 227                 C->seq_al[seq][a]=sub_seq[c];
 228                 c++;
 229               }
 230             else
 231               {
 232                 A->seq_al[seq][a]='-';
 233                 C->seq_al[seq][a]='-';
 234               }
 235           }
 236       }
 237
 238     free_dp_model  (M);
 239     free_aln (S);
 240     free_dp_result (R);
 241     vfree(sub_seq);
 242     vfree(ns);
 243     free_int (ls, -1);
 244     free_sequence (CL->S, (CL->S)->nseq);
 245
 246
 247     return A;
 248   }
 249
 250 Dp_Model * initialize_seg2prf_model(int left_tg_mode, int right_tg_mode, Constraint_list *CL)
 251   {
 252
 253     Dp_Model *M;
 254     int a, b, c,d;
 255
 256     M=vcalloc ( 1, sizeof (Dp_Model));
 257     M->nstate=2;
 258     M->START=M->nstate++;
 259     M->END  =M->nstate++;
 260
 261     M->TG_MODE=1;
 262     M->F_TG_MODE=0;
 263     M->gop=CL->gop*SCORE_K;
 264     M->gep=CL->gep*SCORE_K;
 265
 266     M->bounded_model=declare_int (M->nstate+1, M->nstate+1);
 267     M->model=declare_int (M->nstate+1, M->nstate+1);
 268     for ( a=0; a<=M->nstate; a++)
 269       for ( b=0; b<= M->nstate; b++)
 270         M->model[a][b]=UNDEFINED;
 271
 272     a=0;
 273     M->TYPE=a++;M->LEN_I=a++; M->LEN_J=a++; M->DELTA_I=a++;M->DELTA_J=a++; M->CODING0=a++;M->DELETION=a++;
 274     M->model_properties=declare_int ( M->nstate, 10);
 275
 276     a=0;
 277     M->EMISSION=a++;M->TERM_EMISSION=a++;M->START_EMISSION=a++;
 278     M->model_emission_function=vcalloc(M->nstate, sizeof (int (**)(Alignment*, int **, int, int*, int, int **, int, int*, int, struct Constraint_list *)));
 279     for ( a=0; a< M->nstate; a++)
 280        M->model_emission_function[a]=vcalloc(3, sizeof (int (*)(Alignment*, int **, int, int*, int, int **, int, int*, int, struct Constraint_list *)));
 281
 282
 283     /*Substitution*/
 284     M->model_properties[0][M->TYPE]=M->CODING0;
 285     M->model_properties[0][M->LEN_I]=1;
 286     M->model_properties[0][M->LEN_J]=1;
 287     M->model_properties[0][M->DELTA_I]=-1;
 288     M->model_properties[0][M->DELTA_J]= 0;
 289
 290     M->model_emission_function[0][M->EMISSION]      =cw_profile_get_dp_cost;
 291     M->model_emission_function[0][M->START_EMISSION]=get_start_gep_cost;
 292     M->model_emission_function[0][M->TERM_EMISSION] =get_start_gep_cost;
 293
 294     /*Deletions*/
 295     M->model_properties[1][M->TYPE]=M->DELETION;
 296     M->model_properties[1][M->LEN_I]=1;
 297     M->model_properties[1][M->LEN_J]=0;
 298     M->model_properties[1][M->DELTA_I]=-1;
 299     M->model_properties[1][M->DELTA_J]=+1;
 300     M->model_emission_function[1][M->EMISSION]=get_gep_cost;
 301
 302     if (left_tg_mode ==2)
 303       M->model_emission_function[1][M->START_EMISSION]=get_start_gep_cost;
 304     else M->model_emission_function[1][M->START_EMISSION]=get_gep_cost;
 305
 306     if (right_tg_mode ==2)
 307       M->model_emission_function[1][M->TERM_EMISSION]=get_term_gep_cost;
 308     else M->model_emission_function[1][M->TERM_EMISSION]=get_gep_cost;
 309
 310     /*Transitions*/
 311     M->model[0][M->END]=M->model[M->START][0]=ALLOWED;
 312     M->model[0][1]=M->gop;
 313     M->model[0][0]=ALLOWED;
 314
 315     M->model[1][M->END]=  (right_tg_mode==0)?0:-M->gop;
 316     M->model[M->START][1]=( left_tg_mode==0)?M->gop:0;
 317     M->model[1][1]=ALLOWED;
 318     M->model[1][0]=ALLOWED;
 319
 320
 321
 322
 323     /*Prune the model*/
 324
 325     for (c=0,a=0, d=0; a< M->START; a++)
 326       for ( b=0; b<M->START; b++, d++)
 327         {
 328           if (M->model[a][b]!=UNDEFINED)
 329             {
 330               M->bounded_model[b][1+M->bounded_model[b][0]++]=a;
 331               c++;
 332             }
 333         }
 334     M->CL=CL;
 335
 336     return M;
 337   }
 338
 339 int get_gep_cost (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
 340 {
 341   return CL->gep*SCORE_K;
 342 }
 343
 344 int get_start_gep_cost (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
 345 {
 346   return 0;
 347 }
 348 int get_term_gep_cost (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
 349 {
 350   return CL->gep*SCORE_K*-1;
 351 }
 352
 353
 354
 355
 356 /*********************************COPYRIGHT NOTICE**********************************/
 357 /*© Centro de Regulacio Genomica */
 358 /*and */
 359 /*Cedric Notredame */
 360 /*Tue Oct 27 10:12:26 WEST 2009. */
 361 /*All rights reserved.*/
 362 /*This file is part of T-COFFEE.*/
 363 /**/
 364 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 365 /*    it under the terms of the GNU General Public License as published by*/
 366 /*    the Free Software Foundation; either version 2 of the License, or*/
 367 /*    (at your option) any later version.*/
 368 /**/
 369 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 370 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 371 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 372 /*    GNU General Public License for more details.*/
 373 /**/
 374 /*    You should have received a copy of the GNU General Public License*/
 375 /*    along with Foobar; if not, write to the Free Software*/
 376 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 377 /*...............................................                                                                                      |*/
 378 /*  If you need some more information*/
 379 /*  cedric.notredame@europe.com*/
 380 /*...............................................                                                                                                                                     |*/
 381 /**/
 382 /**/
 383 /*      */
 384 /*********************************COPYRIGHT NOTICE**********************************/