+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdarg.h>
-#include <string.h>
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "define_header.h"
-
-#include "dp_lib_header.h"
-
-int domain_pair_wise (Alignment *A,int*in_ns, int **in_l_s,Constraint_list *CL )
- {
-/*******************************************************************************/
-/* SEQ_DOMAIN DP */
-/* */
-/* makes DP between the the ns[0] sequences and the ns[1] */
-/* */
-/* for MODE, see the function get_dp_cost */
-/*******************************************************************************/
-
- int scale, gop, gep, maximise;
- int a, b, i, j,l,x;
- int best_j;
-
- int lenal[2], len;
- int sub;
- int match;
-
- int *ns, **l_s;
-
-
- int *f;
- int *pf;
-
- int *dd;
- int *dd_len;
-
- int * e;
- int *pe;
- int * e_len;
- int *pe_len;
-
- int fop;
- int **pos0;
-
- int **al=NULL;
- int **pos_al=NULL;
-
-
-
- int ala,LEN;
- char *buffer;
- char *char_buf;
-
-
-/*trace back variables */
- TRACE_TYPE *buf_trace=NULL;
- TRACE_TYPE **trace;
- TRACE_TYPE k;
- TRACE_TYPE *tr;
- int **result_aln;
- int nseq;
-/*Test Varaibles*/
- int score;
-
-/*Prepare l_s and ns*/
-
-
- ns=vcalloc( 2, sizeof (int));
- l_s=vcalloc ( 2, sizeof (int*));
-
- ns[0]=in_ns[1];
- ns[1]=in_ns[0];
-
- l_s[0]=in_l_s[1];
- l_s[1]=in_l_s[0];
- lenal[0]=strlen (A->seq_al[l_s[0][0]]);
- lenal[1]=strlen (A->seq_al[l_s[1][0]]);
- len=lenal[0]+lenal[1]+2;
-/********************************/
-gop=(CL->gop)*SCORE_K;
-gep=(CL->gep)*SCORE_K;
-maximise=CL->maximise;
-scale=(lenal[1]*SCORE_K*(CL->moca)->moca_scale);
-/*******************************/
-
-
-
-/*DO MEMORY ALLOCATION FOR DP*/
-
-
-
-
-
- buf_trace=vcalloc ( len, sizeof (TRACE_TYPE));
- buffer=vcalloc ( 2*len, sizeof (char));
-
- al =declare_int (2, 2*len);
- pos_al=declare_int (2, 2*len);
- result_aln=declare_int (1,len);
- char_buf= vcalloc (2*len, sizeof (char));
-
- f =vcalloc (len, sizeof (int));
- pf =vcalloc (len, sizeof (int));
- e =vcalloc (len, sizeof (int));
- pe =vcalloc (len, sizeof (int));
- e_len =vcalloc (len, sizeof (int));
- pe_len =vcalloc (len, sizeof (int));
- dd =vcalloc (len, sizeof (int));
- dd_len=vcalloc (len, sizeof (int));
-
-
- trace=declare_int (lenal[0]+2, lenal[1]+2);
-
-/*END OF MEMORY ALLOCATION*/
-
-
- /*
- 0(s) +(dd)
- \ |
- \ |
- \ |
- \ |
- \ |
- \ |
- \|
- -(e)----O
- */
-
- pos0=aln2pos_simple ( A,-1, ns, l_s);
-
- for ( i=0; i<=lenal[0]+1; i++)
- {
- tr=trace[i];
-
-
- for ( sub=0,j=0; j<=lenal[1]; j++)
- {
- if (i==0 && j==0){tr[j]=1;pe[j]=dd[j]=gop;}
- else if (i==0) {e[j]=pe[j]=dd[j]=gop;dd_len[j]=e_len[j]=pe_len[j]=f[j]=pf[j]=0;tr[j]=-1;}
- else if (j==0)
- {
- for (f[j]=pf[0],best_j=0,a=1; a<=lenal[1]; a++)
- {
- if (f[j]!=MAX(pf[a]+scale,f[j]))
- {
- f[j]=pf[a]+scale;
- best_j=a;
- }
- }
-
-
- dd [j]=e[j]=pe[j]=gop;
- dd_len[j]=e_len[j]=0;
- tr [j]=best_j;
-
- }
- else if (i>lenal[0]);
- else
- {
- sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
-
-
-
- match=pf[j-1]+sub;
- if (a_better_than_b(pf[j]+gop, pe[j]+gep,maximise))
- {
- e [j]=pf[j]+gop;
- e_len[j]=1;
- }
- else
- {
- e [j]=pe [j]+gep;
- e_len[j]=pe_len[j]+1;
- }
-
-
- if (a_better_than_b(f[j-1]+gop, dd[j-1]+gep,maximise))
- {
- dd [j]=f[j-1]+gop;
- dd_len[j]=1;
- }
- else
- {
- dd [j]=dd [j-1]+gep;
- dd_len[j]=dd_len[j-1]+1 ;
- }
-
-
-
- if ( sub!=UNDEFINED)
- {
-
- f[j] =best_int(4,maximise,&fop,e[j],match,dd[j],f[0]);
- fop-=1;
-
- if (fop==-1)fop= e_len[j]*fop;
- else if (fop==1 )fop=dd_len[j]*fop;
- else if (fop==2 )fop=UNDEFINED;
-
-
-
- }
- else
- {
- dd[j]=e[j]=match=-10000;
- f[j]=f[0];
- fop=UNDEFINED;
- }
- pe [j] =e [j];
- pf [j-1]=f [j-1];
- pe_len[j] =e_len[j];
- tr[j] =fop;
- }
-
- }
-
- pf [j-1]=f [j-1];
- pe [j] =e [j];
- pe_len[j] =e_len[j];
- }
- score=f[0];
- i=lenal[0]+1;
- j=0;
- ala=0;
-
-
- while (i!=0)
- {
-
-
- k=trace[i][j];
- if (j==0 && i<=lenal[0])
- {
-
- pos_al[0][ala]=i;
- pos_al[1][ala]=j;
- al[0][ala]=MATCH;
- al[1][ala]=UNALIGNED;
- i--;
- j=k;
- ala++;
- }
- else if ( j==0 && i>lenal[0])
- {
- j=k;
- i--;
-
- }
- else if (k==0)
- {
- pos_al[0][ala]=i;
- pos_al[1][ala]=j;
- al[0][ala]=MATCH;
- al[1][ala]=MATCH;
- i--;
- j--;
-
- ala++;
- }
- else if (k!=UNDEFINED && k<0 && i>0)
- {
- for (x=0; x< -k && i>0; x++)
- {
- pos_al[0][ala]=i;
- pos_al[1][ala]=j;
- al[0][ala]=MATCH;
- al[1][ala]=GAP;
- i--;
- ala++;
- }
- }
- else if (k!=UNDEFINED && k>0 && j>0)
- {
- for ( x=0; x< k && j>0; x++)
- {
- pos_al[0][ala]=i;
- pos_al[1][ala]=j;
- al[0][ala]=GAP;
- al[1][ala]=MATCH;
- j--;
- ala++;
- }
- }
- else if ( k==UNDEFINED){j=0;}
-
- }
-
-
- LEN=ala;
-
- invert_list_int ( pos_al[0], LEN);
- invert_list_int ( pos_al[1], LEN);
- invert_list_int ( al[0], LEN);
- invert_list_int ( al[1], LEN);
-
-
- /*O: TARGET SEQUENCE (long)*/
- /*1: PATTERN SEQUENCE (short)*/
-
-
-
- for ( b=0; b<lenal[1]; b++)result_aln[0][b]=-1;
- for (l=0, nseq=0, a=0; a< LEN; a++)
- {
- i=pos_al[0][a];
- j=pos_al[1][a];
-
- if (al[1][a]==UNALIGNED && l>0)
- {
- result_aln=realloc_int ( result_aln, read_size_int ( result_aln,sizeof (int*)),len, 1, 0);
- nseq++;
- l=0;
- for ( b=0; b<lenal[1]; b++)result_aln[nseq][b]=-1;
- }
-
- else if (al[1][a]==MATCH && al[0][a]==MATCH ){l++;result_aln[nseq][j-1]=i-1;}
- else if (al[1][a]==MATCH && al[0][a]==GAP ){l++;result_aln[nseq][j-1]=-1;}
-
- }
- if ( l>0)nseq++;
-
-
-
-
- A=domain_match_list2aln ( A,ns,l_s,result_aln,nseq,lenal[1]);
-
-
- vfree (f);
- vfree (pf);
- vfree (e);
- vfree (pe);
- vfree (e_len);
- vfree (pe_len);
- vfree (dd_len);
- vfree (dd);
- free_int (pos0, -1);
- vfree (buffer);
- vfree (char_buf);
- vfree (buf_trace);
- free_int ( pos_al, -1);
- pos_al=NULL;
-
- free_int ( al, -1);
-
- free_int (trace,-1);
- free_int ( result_aln, -1);
- return score;
- }
-
-
-Alignment *domain_match_list2aln ( Alignment *A,int *ns,int **l_s,int **ml, int nseq, int len)
- {
-
- /*
- function documentation: start
-
- This function edits the alignment given the results obtained by DP
- ns: ns[0]->number of sequences serarched (TARGET)
- ns[1]->number of sequences in the pattern (PATTERN SEQ)
- l_s:
- l_s[0]->list of sequences in the TARGET...
-
- nseq: number of occurences of PATTERN in TARGET
- len: length of the PATTERN
-
- ml: detail of the nseq matches
- ml[x][y]=k-> residue k of TARGET matches residue y of pattern
- -> k=-1 means a gap;
-
- NOTE: This implementation can only match ONE target sequence with the PATTERN
- The Pattern can either be one sequence or a profile.
-
- function documentation: end
- */
-
-
-
-
- int a, b, c, d, e;
- Alignment *B=NULL;
- int **new_ml;
- int *max_ml;
- int *start_ml;
- int tot_nseq;
- int max_len,seq;
- char *buf;
-
- if ( len==0 || nseq==0)
- {
- A->nseq=0;
- A->len_aln=0;
- }
- else
- {
- B=copy_aln(A, B);
- /*1 Extract the sequence used as a pattern, put it on the top*/
-
-
- A=shrink_aln (A, ns[1], l_s[1]);
- A=realloc_aln2(A, ns[1]+ns[0]*nseq,len+A->len_aln+1);
-
-
-
- new_ml =declare_int ( nseq, 3*len);
- max_ml =vcalloc ( nseq, sizeof (int));
- for ( a=0; a<nseq; a++)max_ml[a]=-1;
-
- start_ml=vcalloc ( nseq, sizeof (int));
-
- for ( b=0,a=0; a< len; a++, b+=3)
- {
- for ( max_len=0,c=0; c< nseq; c++)
- {
- if ( max_ml[c]<0 && ml[c][a]<0)
- {
- new_ml[c][b]=ml[c][a];
- new_ml[c][b+1]=0;
- }
- else if ( ml[c][a]>=0)
- {
- new_ml[c][b]=ml[c][a];
- if ( max_ml[c]<0) start_ml[c]=max_ml[c]=ml[c][a];
- for ( d=a+1; d<len; d++){if ( ml[c][d]>=0){ max_ml[c]= ml[c][d];break;}}
- if (max_ml[c]!=new_ml[c][b])
- {
- new_ml[c][b+1]=max_ml[c]-new_ml[c][b]-1;
- }
- max_len=MAX( max_len, new_ml[c][b+1]);
- }
- else
- {
- new_ml[c][b]=ml[c][a];
- new_ml[c][b+1]=0;
- }
- }
-
- for ( c=0; c< nseq; c++){new_ml[c][b+2]=max_len;}
- }
-
- tot_nseq=ns[1]+ns[0]*nseq;
-
- for ( a=0, b=0; a< len ;a++)
- {
-
- /*1: Place the Match Column*/
- for ( c=0; c< ns[1]; c++)
- {
- A->seq_al[c][b]=B->seq_al[l_s[1][c]][a];
- A->seq_al[c][b+1]='\0';
- }
- for ( e=0,c=ns[1]; c<tot_nseq; c+=ns[0],e++)
- for ( d=0; d< ns[0]; d++)
- {
- if ( new_ml[e][3*a]!=-1)
- A->seq_al[c+d][b]=B->seq_al[l_s[0][d]][new_ml[e][3*a]];
- else
- A->seq_al[c+d][b]='-';
- A->seq_al[c+d][b+1]='\0';
-
- }
- b++;
-
- /*2: Add the Gaps before the next_column*/
- if ( new_ml[0][3*a+2]>0)
- {
- for ( c=0; c< ns[1]; c++)
- {
- buf=generate_null(new_ml[0][3*a+2]);
- strcat ( A->seq_al[c],buf);
- vfree (buf);
- }
- for (e=0,c=ns[1];c< tot_nseq; c+=ns[0], e++)
- {
- buf=extract_char (B->seq_al[l_s[0][0]], new_ml[e][3*a]+1, new_ml[e][3*a+1]);
- strcat ( A->seq_al[c],buf);
- vfree (buf);
- buf=generate_null(new_ml[e][3*a+2]-new_ml[e][3*a+1]);
- strcat ( A->seq_al[c],buf);
-
- vfree (buf);
-
- }
- }
- b+=new_ml[0][3*a+2];
- }
-
- for (e=0,a=ns[1]; a< tot_nseq; a+=ns[0],e++)
- {
- for ( b=0; b<ns[0]; b++)
- {
- seq=l_s[0][b];
- A->order[a+b][0]=B->order[seq][0];
- A->order[a+b][1]=B->order[seq][1];
- for ( c=0; c<start_ml[e];c++)A->order[a+b][1]+=!is_gap(B->seq_al[seq][c]);
- sprintf ( A->name[a+b], "Repeat_%d", a+b);
- }
- }
-
- free_aln(B);
- A->nseq=tot_nseq;
- A->len_aln=strlen ( A->seq_al[0]);
-
- }
- return A;
- }
-Alignment * domain_seq2domain (Constraint_list *CL,int scale,int gop,int gep,Alignment *SEQ_DOMAIN, Alignment *TARGET)
- {
- static Alignment *A;
- int *n_groups;
- int **group_list;
- int a,b,c;
-
- A=copy_aln (TARGET, A);
- A=stack_aln( A, SEQ_DOMAIN);
-
-
- n_groups=vcalloc ( 2, sizeof (int));
- group_list=declare_int (2, A->nseq);
-
- n_groups[0]=TARGET->nseq;
- n_groups[1]=SEQ_DOMAIN->nseq;
- for (c=0, a=0; a< 2; a++)
- {
- for (b=0; b< n_groups[a]; b++, c++)
- {
- group_list[a][b]=c;
- }
- }
- A->score_aln=domain_pair_wise (A, n_groups, group_list,CL);
-
- SEQ_DOMAIN=copy_aln (A, SEQ_DOMAIN);
- vfree (n_groups);
- free_int (group_list,-1);
- return SEQ_DOMAIN;
- }
-
-
-/******************************COPYRIGHT NOTICE*******************************/
-/*© Centro de Regulacio Genomica */
-/*and */
-/*Cedric Notredame */
-/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
-/*All rights reserved.*/
-/*This file is part of T-COFFEE.*/
-/**/
-/* T-COFFEE is free software; you can redistribute it and/or modify*/
-/* it under the terms of the GNU General Public License as published by*/
-/* the Free Software Foundation; either version 2 of the License, or*/
-/* (at your option) any later version.*/
-/**/
-/* T-COFFEE is distributed in the hope that it will be useful,*/
-/* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
-/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
-/* GNU General Public License for more details.*/
-/**/
-/* You should have received a copy of the GNU General Public License*/
-/* along with Foobar; if not, write to the Free Software*/
-/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
-/*............................................... |*/
-/* If you need some more information*/
-/* cedric.notredame@europe.com*/
-/*............................................... |*/
-/**/
-/**/
-/* */
-/******************************COPYRIGHT NOTICE*******************************/
git://source.jalview.org / jabaws.git / blobdiff