binaries/src/tcoffee/t_coffee_source/util_domain_constraints_list.c

   1 #include <stdlib.h>
   2 #include <stdio.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 /*********************************************************************/
  11 /*                                                                   */
  12 /*                         MASKING LIST FUNCTIONS                    */
  13 /*                                                                   */
  14 /*                                                                   */
  15 /*********************************************************************/
  16 Constraint_list * mask_list_with_aln (Alignment *A,int start, int len,Constraint_list *CL, int new_value)
  17      {
  18      int a, c,d;
  19      static int *entry;
  20      int s1, s2, r1, r2;
  21      int **pos;
  22      int **cache;
  23      int max_nseq;
  24      int max_len;
  25
  26
  27
  28      if ( A==NULL || A->len_aln==0 || A->nseq<=1)return CL;
  29      if ( entry==NULL) entry=vcalloc (CL->entry_len , CL->el_size);
  30
  31      cache=declare_int (return_max_int (A->order,read_size_int ( A->order,sizeof (int*)),0)+1,return_max_int (A->order,read_size_int ( A->order,sizeof (int)),1)+A->len_aln+1);
  32
  33      max_nseq=return_max_int (A->order,read_size_int ( A->order,sizeof (int*)),0)+1;
  34      max_len =return_max_int (A->order,read_size_int ( A->order,sizeof (int*)),1)+A->len_aln+1;
  35
  36      pos=aln2pos_simple(A, A->nseq);
  37
  38      for (a=0; a< A->nseq; a++)
  39          for (c=start; c< (start+len); c++)
  40              if ( pos[a][c]>0)cache[A->order[a][0]][pos[a][c]]=1;
  41
  42
  43      if (!CL->M)
  44           {
  45           for ( d=0; d<CL->ne; d++)
  46                      {
  47                      s1=vread_clist(CL, d, SEQ1);
  48                      s2=vread_clist(CL, d, SEQ2);
  49                      r1=vread_clist(CL, d, R1);
  50                      r2=vread_clist(CL, d, R2);
  51                      if ( s1>=max_nseq);
  52                      else if (r1>=max_len);
  53                      else if ( cache[s1][r1]==1)
  54                          {
  55                          mask_entry(CL,d,new_value);
  56                          }
  57
  58                      if ( s2>=max_nseq);
  59                      else if (r2>=max_len);
  60                      else if ( cache[s2][r2]==1)
  61                               {
  62                               mask_entry(CL,d,new_value);
  63                               }
  64                      }
  65
  66           sort_constraint_list_inv (CL,0, CL->ne);
  67           sort_constraint_list (CL,0, CL->ne);
  68           free_int ( cache, -1);
  69           }
  70      else if ( CL->M)
  71           {
  72
  73           for (a=0; a< A->nseq; a++)
  74               for (c=start; c< (start+len); c++)
  75                   if ( pos[a][c]>0)
  76                       {
  77                       vwrite_clist(CL,30+A->order[a][0],pos[a][c],UNDEFINED);
  78                       }
  79           }
  80      free_int (pos, -1);
  81      return CL;
  82      }
  83 Constraint_list* mask_list_with_aln_pair (Alignment *A,int start, int len ,Constraint_list *CL,int new_value)
  84      {
  85      int a, b, p;
  86      int *entry;
  87
  88
  89      int l1, l2, r1, r2, s1, s2;
  90      int x, y;
  91
  92      if ( A==NULL || A->len_aln==0 || A->nseq<=1)return CL;
  93
  94      if (CL->M)
  95         {
  96         fprintf ( stderr, "\nERROR: AA matrix cannot be masked with  mask_list_with_aln_pair");
  97         myexit (EXIT_SUCCESS);
  98         }
  99
 100
 101      entry=vcalloc (CL->entry_len, sizeof (int));
 102
 103      for ( a=0; a< A->nseq-1; a++)
 104          {
 105          for (l1=A->order[a][1]+1,p=0    ; p<start        ; p++)l1+=!is_gap(A->seq_al[a][p]);
 106          for (r1=l1-1            ,p=start; p<(start+len)  ; p++)r1+=!is_gap(A->seq_al[a][p]);
 107          s1=A->order[a][0];
 108
 109          for ( b=a+1; b< A->nseq; b++)
 110              {
 111              for (l2=A->order[b][1]+1,p=0    ; p<start; p++)l2+=!is_gap(A->seq_al[b][p]);
 112              for (r2=l2-1            ,p=start; p<(start+len)  ; p++)r2+=!is_gap(A->seq_al[b][p]);
 113              s2=A->order[b][0];
 114
 115
 116              for ( x=l1; x<=r1; x++)
 117                  {
 118
 119                  for ( y=l2; y<=r2; y++)
 120                      {
 121
 122                      set_int(entry,4,x,R1,y,R2,s1,SEQ1,s2,SEQ2);
 123                      if ( (main_search_in_list_constraint ( entry,&p,4,CL))!=NULL)
 124                           {
 125                            mask_entry(CL,p,new_value);
 126                           }
 127                      set_int(entry,4,y,R1,x,R2,s2,SEQ1,s1,SEQ2);
 128                      if ( (main_search_in_list_constraint ( entry,&p,4,CL))!=NULL)
 129                           {
 130
 131                           mask_entry(CL,p,new_value);
 132                           }
 133                      }
 134                  }
 135              }
 136          }
 137
 138
 139      vfree(entry);
 140
 141      return CL;
 142      }
 143
 144 Constraint_list *mask_entry( Constraint_list *CL, int p, int new_value)
 145      {
 146      vwrite_clist(CL, p, WE,   new_value);
 147      vwrite_clist(CL, p, CONS, new_value);
 148      vwrite_clist(CL, p, MISC, new_value);
 149      return CL;
 150      }
 151 /*********************************************************************/
 152 /*                                                                   */
 153 /*                         SEQUENCE CONCATENATION                    */
 154 /*                                                                   */
 155 /*                                                                   */
 156 /*********************************************************************/
 157 Constraint_list *prepare_list_and_seq4sw(Constraint_list *I, int n_seq, char **seq_name)
 158         {
 159         int a, b;
 160         int len,l;
 161         char **long_seq=NULL;
 162         int  **translation;
 163         char **name;
 164         int s1, s2;
 165         Constraint_list *Out=NULL;
 166         Sequence *S_Out;
 167
 168         translation=declare_int ( (I->S)->nseq,2);
 169         name=declare_char (1,STRING);
 170         long_seq=declare_char(1, STRING);
 171
 172         for (len=0,a=0; a< (I->S)->nseq; a++)
 173             {
 174             if((b=name_is_in_list ((I->S)->name[a],seq_name, n_seq, 100))!=-1)
 175                {
 176                l=strlen((I->S)->seq[a])+1;
 177                long_seq[0]=vrealloc(long_seq[0],(len+l+1)*sizeof(char));
 178                long_seq[0]=strcat(long_seq[0], (I->S)->seq[a]);
 179                long_seq[0]=strcat(long_seq[0], "O");
 180
 181                translation[a][0]=b;
 182                translation[a][1]=len;
 183                len+=l;
 184                }
 185             else translation[a][0]=-1;
 186             }
 187
 188         long_seq[0][len-1]='\0';
 189         len--;
 190         sprintf ( name[0], "concatenat");
 191         S_Out=fill_sequence_struc(1, long_seq, name);
 192         free_char(name, -1);
 193         free_char(long_seq, -1);
 194
 195
 196         if (!I->M)
 197            {
 198            if ( I->fp)    Out=declare_constraint_list(S_Out, NULL, NULL, 0, vtmpfile(),NULL);
 199            else if ( I->L)Out=declare_constraint_list(S_Out, NULL, NULL, 0, NULL     ,NULL);
 200
 201            for (a=0; a<I->ne; a++)
 202                {
 203                s1=vread_clist(I,a,SEQ1);
 204                s2=vread_clist(I,a,SEQ2);
 205
 206                if ( translation[s1][0]!=-1 &&  translation[s2][0]!=-1)
 207                    Out=add_list_entry2list(Out, Out->entry_len, SEQ1, 0, SEQ2, 0, R1,vread_clist(I,a,R1)+translation[s1][1], R2,vread_clist(I,a,R2)+translation[s2][1], WE,vread_clist(I,a,WE),CONS, vread_clist(I,a,CONS), MISC, vread_clist(I,a,MISC));
 208                }
 209
 210            for ( a=0; a<(I->S)->nseq; a++)
 211                {
 212                if (translation[a][0]!=-1 && translation[a][1]!=0)
 213                   {
 214                   for (b=1; b<=len; b++)
 215                       {
 216                       add_list_entry2list(Out,Out->entry_len, SEQ1, 0, SEQ2, 0, R1, translation[a][1], R2, b, WE, UNDEFINED, CONS, 0,MISC, vread_clist(I,a,MISC));
 217                       }
 218                   }
 219                }
 220            sort_constraint_list (Out, 0, Out->ne);
 221            }
 222         else if (I->M)
 223            {
 224            Out=declare_constraint_list(S_Out, NULL, NULL, 0, vtmpfile(),I->M);
 225            vfree((Out->M)[30]);
 226            (Out->M)[30]=vcalloc ( len+1, sizeof (int));
 227            for ( a=0; a<len; a++)
 228                {
 229                if ( long_seq[0][a]=='O')vwrite_clist(Out, 30, a+1,UNDEFINED);
 230                }
 231            Out->ne=SIZEOF_AA_MAT;
 232            }
 233         free_int (translation,-1);
 234         return Out;
 235         }
 236 /*********************************************************************/
 237 /*                                                                   */
 238 /*                         MISCEANELLOUS                             */
 239 /*                                                                   */
 240 /*                                                                   */
 241 /*********************************************************************/
 242 int ** get_undefined_list (Constraint_list *CL)
 243          {
 244          int **list;
 245          int a;
 246          CLIST_TYPE x;
 247
 248          list=declare_int ( (CL->S)->nseq+1, (CL->S)->max_len+1);
 249
 250          for ( a=0; a< CL->ne; a++)
 251              {
 252              x=vread_clist(CL, a, WE);
 253              list[vread_clist(CL, a, SEQ1)][vread_clist(CL, a, R1)]=(x==UNDEFINED);
 254              list[vread_clist(CL, a, SEQ2)][vread_clist(CL, a, R2)]=(x==UNDEFINED);
 255              }
 256          return list;
 257          }
 258 int      is_never_undefined (Constraint_list *CL,int r)
 259          {
 260          int a;
 261          for ( a=0; a< CL->ne; a++)
 262              {
 263              if ( (vread_clist(CL,a,R1)==r || vread_clist(CL,a,R2)==r) && vread_clist(CL,a,WE)==UNDEFINED)return 0;
 264              }
 265          return 1;
 266          }
 267
 268 int* do_analyse_list ( Constraint_list *CL)
 269         {
 270         int **seq_score;
 271         int **seq_score2;
 272         int  *pos_L;
 273         int a;
 274         int n_res;
 275
 276
 277         int n_it=4;
 278
 279         double sum, sum2, tot;
 280         int field=2;
 281         double z;
 282         int r1, r2;
 283         int max_we=0;
 284
 285
 286
 287         fprintf ( stderr, "\nDO ANALYSE");
 288
 289         n_res=(CL->S)->max_len;
 290
 291         pos_L =vcalloc     (n_res+2, sizeof (int));
 292         pos_L++;
 293         pos_L[-1]=n_res;
 294         seq_score=declare_int (n_res+1,n_it+1);
 295         seq_score2=declare_int(n_res+1,n_it+1);
 296
 297
 298
 299         for ( a=0; a< CL->ne; a++)
 300             {
 301             r1=vread_clist(CL, a, R1);
 302             r2=vread_clist(CL, a, R2);
 303
 304             seq_score[r1][0]+=vread_clist(CL, a, WE);
 305             seq_score[r2][0]+=vread_clist(CL, a, WE);
 306
 307
 308             seq_score[r1][1]+=vread_clist(CL, a, CONS);
 309             seq_score[r2][1]+=vread_clist(CL, a, CONS);
 310
 311
 312             seq_score[r1][2]+=vread_clist(CL, a, MISC);
 313             seq_score[r2][2]+=vread_clist(CL, a, MISC);
 314
 315             }
 316         for ( a=1; a<=n_res; a++)max_we=MAX(seq_score[a][0],max_we);
 317         for ( a=1; a<=n_res; a++)
 318             {
 319             if ( a!=n_res && seq_score[a][0]> seq_score[a-1][0]) fprintf ( stderr, "\n%4d %s", a, num2plot(seq_score[a][0],max_we,40));
 320             else fprintf ( stderr, "\n");
 321             }
 322
 323         for ( a=0; a< CL->ne; a++)if ( vread_clist(CL, a, MISC)>1000)fprintf ( stderr, "\n%4d %4d %4d", vread_clist(CL, a,R1), vread_clist(CL, a, R2), vread_clist(CL, a, MISC));
 324         myexit (EXIT_SUCCESS);
 325
 326         for ( a=0; a<CL->ne; a++)
 327             {
 328             if ( vread_clist(CL, a, WE)!=UNDEFINED &&vread_clist(CL, a, MISC)>=1 )
 329                 {
 330                 seq_score[vread_clist(CL, a, R1)][0]+=vread_clist(CL, a, CONS);
 331                 seq_score[vread_clist(CL, a, R1)][1]+=vread_clist(CL, a, WE);
 332                 seq_score[vread_clist(CL, a, R1)][2]+=vread_clist(CL, a, MISC);
 333
 334                 seq_score[vread_clist(CL, a, R2)][0]+=vread_clist(CL, a, CONS);
 335                 seq_score[vread_clist(CL, a, R2)][1]+=vread_clist(CL, a, WE);
 336                 seq_score[vread_clist(CL, a, R2)][2]+=vread_clist(CL, a, MISC);
 337                 }
 338             }
 339
 340         for (a=1, tot=0,sum=0, sum2=0; a<= n_res  ; a++)
 341             {
 342             if ( seq_score[a][2]>0)
 343                 {
 344                 sum +=seq_score[a][field];
 345                 sum2+=seq_score[a][field]*seq_score[a][field];
 346                 tot++;
 347                 }
 348             }
 349
 350         fprintf ( stderr, "\n");
 351
 352
 353
 354         for (a=1; a<= n_res  ; a++)
 355             {
 356             z=return_z_score (seq_score[a][field],sum, sum2, tot );
 357             if ( seq_score[a][2]>0)
 358                 {
 359
 360                 pos_L[a]=(int)(z*10);
 361                 pos_L[0]++;
 362
 363                 }
 364             }
 365
 366
 367         free_int (seq_score2,-1);
 368         free_int (seq_score,-1);
 369         return pos_L;
 370         }
 371 /*********************************COPYRIGHT NOTICE**********************************/
 372 /*© Centro de Regulacio Genomica */
 373 /*and */
 374 /*Cedric Notredame */
 375 /*Tue Oct 27 10:12:26 WEST 2009. */
 376 /*All rights reserved.*/
 377 /*This file is part of T-COFFEE.*/
 378 /**/
 379 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
 380 /*    it under the terms of the GNU General Public License as published by*/
 381 /*    the Free Software Foundation; either version 2 of the License, or*/
 382 /*    (at your option) any later version.*/
 383 /**/
 384 /*    T-COFFEE is distributed in the hope that it will be useful,*/
 385 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
 386 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
 387 /*    GNU General Public License for more details.*/
 388 /**/
 389 /*    You should have received a copy of the GNU General Public License*/
 390 /*    along with Foobar; if not, write to the Free Software*/
 391 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
 392 /*...............................................                                                                                      |*/
 393 /*  If you need some more information*/
 394 /*  cedric.notredame@europe.com*/
 395 /*...............................................                                                                                                                                     |*/
 396 /**/
 397 /**/
 398 /*      */
 399 /*********************************COPYRIGHT NOTICE**********************************/