7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 /*********************************************************************/
12 /* MASKING LIST FUNCTIONS */
15 /*********************************************************************/
16 Constraint_list * mask_list_with_aln (Alignment *A,int start, int len,Constraint_list *CL, int new_value)
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);
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);
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;
36 pos=aln2pos_simple(A, A->nseq);
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;
45 for ( d=0; d<CL->ne; d++)
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);
52 else if (r1>=max_len);
53 else if ( cache[s1][r1]==1)
55 mask_entry(CL,d,new_value);
59 else if (r2>=max_len);
60 else if ( cache[s2][r2]==1)
62 mask_entry(CL,d,new_value);
66 sort_constraint_list_inv (CL,0, CL->ne);
67 sort_constraint_list (CL,0, CL->ne);
68 free_int ( cache, -1);
73 for (a=0; a< A->nseq; a++)
74 for (c=start; c< (start+len); c++)
77 vwrite_clist(CL,30+A->order[a][0],pos[a][c],UNDEFINED);
83 Constraint_list* mask_list_with_aln_pair (Alignment *A,int start, int len ,Constraint_list *CL,int new_value)
89 int l1, l2, r1, r2, s1, s2;
92 if ( A==NULL || A->len_aln==0 || A->nseq<=1)return CL;
96 fprintf ( stderr, "\nERROR: AA matrix cannot be masked with mask_list_with_aln_pair");
97 myexit (EXIT_SUCCESS);
101 entry=vcalloc (CL->entry_len, sizeof (int));
103 for ( a=0; a< A->nseq-1; a++)
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]);
109 for ( b=a+1; b< A->nseq; b++)
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]);
116 for ( x=l1; x<=r1; x++)
119 for ( y=l2; y<=r2; y++)
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)
125 mask_entry(CL,p,new_value);
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)
131 mask_entry(CL,p,new_value);
144 Constraint_list *mask_entry( Constraint_list *CL, int p, int new_value)
146 vwrite_clist(CL, p, WE, new_value);
147 vwrite_clist(CL, p, CONS, new_value);
148 vwrite_clist(CL, p, MISC, new_value);
151 /*********************************************************************/
153 /* SEQUENCE CONCATENATION */
156 /*********************************************************************/
157 Constraint_list *prepare_list_and_seq4sw(Constraint_list *I, int n_seq, char **seq_name)
161 char **long_seq=NULL;
165 Constraint_list *Out=NULL;
168 translation=declare_int ( (I->S)->nseq,2);
169 name=declare_char (1,STRING);
170 long_seq=declare_char(1, STRING);
172 for (len=0,a=0; a< (I->S)->nseq; a++)
174 if((b=name_is_in_list ((I->S)->name[a],seq_name, n_seq, 100))!=-1)
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");
182 translation[a][1]=len;
185 else translation[a][0]=-1;
188 long_seq[0][len-1]='\0';
190 sprintf ( name[0], "concatenat");
191 S_Out=fill_sequence_struc(1, long_seq, name);
193 free_char(long_seq, -1);
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);
201 for (a=0; a<I->ne; a++)
203 s1=vread_clist(I,a,SEQ1);
204 s2=vread_clist(I,a,SEQ2);
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));
210 for ( a=0; a<(I->S)->nseq; a++)
212 if (translation[a][0]!=-1 && translation[a][1]!=0)
214 for (b=1; b<=len; b++)
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));
220 sort_constraint_list (Out, 0, Out->ne);
224 Out=declare_constraint_list(S_Out, NULL, NULL, 0, vtmpfile(),I->M);
226 (Out->M)[30]=vcalloc ( len+1, sizeof (int));
227 for ( a=0; a<len; a++)
229 if ( long_seq[0][a]=='O')vwrite_clist(Out, 30, a+1,UNDEFINED);
231 Out->ne=SIZEOF_AA_MAT;
233 free_int (translation,-1);
236 /*********************************************************************/
241 /*********************************************************************/
242 int ** get_undefined_list (Constraint_list *CL)
248 list=declare_int ( (CL->S)->nseq+1, (CL->S)->max_len+1);
250 for ( a=0; a< CL->ne; a++)
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);
258 int is_never_undefined (Constraint_list *CL,int r)
261 for ( a=0; a< CL->ne; a++)
263 if ( (vread_clist(CL,a,R1)==r || vread_clist(CL,a,R2)==r) && vread_clist(CL,a,WE)==UNDEFINED)return 0;
268 int* do_analyse_list ( Constraint_list *CL)
279 double sum, sum2, tot;
287 fprintf ( stderr, "\nDO ANALYSE");
289 n_res=(CL->S)->max_len;
291 pos_L =vcalloc (n_res+2, sizeof (int));
294 seq_score=declare_int (n_res+1,n_it+1);
295 seq_score2=declare_int(n_res+1,n_it+1);
299 for ( a=0; a< CL->ne; a++)
301 r1=vread_clist(CL, a, R1);
302 r2=vread_clist(CL, a, R2);
304 seq_score[r1][0]+=vread_clist(CL, a, WE);
305 seq_score[r2][0]+=vread_clist(CL, a, WE);
308 seq_score[r1][1]+=vread_clist(CL, a, CONS);
309 seq_score[r2][1]+=vread_clist(CL, a, CONS);
312 seq_score[r1][2]+=vread_clist(CL, a, MISC);
313 seq_score[r2][2]+=vread_clist(CL, a, MISC);
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++)
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");
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);
326 for ( a=0; a<CL->ne; a++)
328 if ( vread_clist(CL, a, WE)!=UNDEFINED &&vread_clist(CL, a, MISC)>=1 )
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);
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);
340 for (a=1, tot=0,sum=0, sum2=0; a<= n_res ; a++)
342 if ( seq_score[a][2]>0)
344 sum +=seq_score[a][field];
345 sum2+=seq_score[a][field]*seq_score[a][field];
350 fprintf ( stderr, "\n");
354 for (a=1; a<= n_res ; a++)
356 z=return_z_score (seq_score[a][field],sum, sum2, tot );
357 if ( seq_score[a][2]>0)
360 pos_L[a]=(int)(z*10);
367 free_int (seq_score2,-1);
368 free_int (seq_score,-1);
371 /*********************************COPYRIGHT NOTICE**********************************/
372 /*© Centro de Regulacio Genomica */
374 /*Cedric Notredame */
375 /*Tue Oct 27 10:12:26 WEST 2009. */
376 /*All rights reserved.*/
377 /*This file is part of T-COFFEE.*/
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.*/
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.*/
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 /*............................................... |*/
399 /*********************************COPYRIGHT NOTICE**********************************/