+++ /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"
-
-/******************************************************************/
-/* MOCA DRIVER */
-/* */
-/* */
-/******************************************************************/
-Constraint_list *prepare_cl_for_moca ( Constraint_list *CL)
- {
- int a, b, c;
- int tot_l, l;
- char **name, **seq;
- Sequence *NS=NULL;
-
- /*Prepare the constraint list*/
- CL->do_self=1;
- CL->get_dp_cost=moca_slow_get_dp_cost;
- CL->evaluate_residue_pair=moca_residue_pair_extended_list;
-
- /*Prepare the moca parameters*/
- (CL->moca)->evaluate_domain=evaluate_moca_domain;
- (CL->moca)->cache_cl_with_domain=cache_cl_with_moca_domain;
- (CL->moca)->make_nol_aln=make_moca_nol_aln;
-
- /*Prepare the packing of the sequences*/
- for ( a=0, b=1; a< (CL->S)->nseq; a++)b+=strlen ( (CL->S)->seq[a])+1;
-
- seq =declare_char ( 1,b+1);
- name=declare_char( 1,30);
- CL->packed_seq_lu =declare_int ( b, 2);
-
-
- for (tot_l=1,a=0; a< (CL->S)->nseq; a++)
- {
- strcat (seq[0], (CL->S)->seq[a]);
- strcat (seq[0], "X");
- l=strlen((CL->S)->seq[a]);
- for ( c=1; c<= l; c++, tot_l++)
- {
- CL->packed_seq_lu[tot_l][0]=a;
- CL->packed_seq_lu[tot_l][1]=c;
- }
- CL->packed_seq_lu[tot_l++][0]=UNDEFINED;
- }
- sprintf ( name[0], "catseq");
- NS=fill_sequence_struc(1, seq, name);
- CL->S=add_sequence (NS, CL->S, 0);
- free_char( seq, -1);
- free_char(name, -1);
- free_sequence (NS, NS->nseq);
-
-
- return CL;
- }
-
-Alignment ** moca_aln ( Constraint_list *CL)
- {
- /*
- function documentation: start
-
- Alignment ** moca_aln ( Constraint_list *CL)
-
- This function inputs CL and outputs a series of local multiple alignments
- contained in aln_list;
-
- The terminator of aln_list is set to NULL;
-
- function documentation: end
- */
-
-
- static int max_n_domains=1000;
- int n_domains=0;
-
- Alignment **aln_list;
-
-
-
- aln_list=vcalloc (max_n_domains, sizeof (Alignment *));
- if ((CL->moca)->moca_interactive)aln_list[n_domains++]=extract_domain ( CL);
- else
- {
- while ( (aln_list[n_domains++]=extract_domain ( CL))!=NULL)
- {
- if ((CL->moca)->moca_len)break;
- if ( n_domains==max_n_domains)
- {
- n_domains+=1000;
- aln_list=vrealloc (aln_list, max_n_domains*sizeof (Alignment*));
- }
- }
- }
- return aln_list;
- }
-
-Alignment * extract_domain ( Constraint_list *CL)
- {
- /*
- function documentation: start
- Alignment * extract_domain ( Constraint_list *CL)
-
- given a CL, this function extracts the next best scoring local multiple alignment
- It returns a CL where the aligned residues have been indicated in (CL->moca)->forbiden_residues;
-
- the local alignment is extracted with the dp function indicated by
- CL->dp_mode: (gotoh_sw_pair_wise)
- Evaluation:
- CL->get_dp_cost=slow_get_dp_cost;
- CL->evaluate_residue_pair=sw_residue_pair_extended_list;
- Continuation:
- (CL->moca)->evaluate_domain=evaluate_moca_domain;
- Cache of CL:
- (CL->moca)->cache_cl_with_domain=cache_cl_with_moca_domain;
- Domain post processing:
- (CL->moca)->make_nol_aln=make_moca_nol_aln;
- function documentation: end
- */
- int min_start, max_start, start,min_len, max_len, len, score;
- int step;
- Alignment *C=NULL;
- Alignment *RESULT=NULL;
- Alignment *EA=NULL;
-
-
-
-
- /*CASE 1: Non Automatic Domain Extraction*/
- if ((CL->moca)->moca_interactive)
- {
- return interactive_domain_extraction (CL);
- }
- else if ((CL->moca)->moca_len)
- {
- while ((C=extract_domain_with_coordinates (C,(CL->moca)->moca_start,(CL->moca)->moca_len,CL))->nseq==0)(CL->moca)->moca_scale=(CL->moca)->moca_scale*0.9;
- RESULT=copy_aln ( C, RESULT);
- unpack_seq_aln (RESULT, CL);
- output_format_aln ("mocca_aln",RESULT,EA=fast_coffee_evaluate_output(RESULT, CL),"stdout");
- free_aln(EA);
-
- return RESULT;
- }
- else if ( !(CL->moca)->moca_len)
- {
- analyse_sequence (CL);
- myexit (EXIT_FAILURE);
- }
-
- /*CASE 2: Automatic Domain Extraction: Find Coordinates*/
-
-
- start=500;
-
- step=10;
- min_start=0;
- max_start=strlen ((CL->S)->seq[0]);
- min_len=20;
- max_len=strlen ((CL->S)->seq[0]);
-
- C=extract_domain_with_coordinates (C,13,30,CL);
- C->output_res_num=1;
- print_aln (C);
-
- (CL->moca)->moca_scale=-180;
- C=add_seq2aln (CL,C, CL->S);
- print_aln (C);
-
- (CL->moca)->moca_scale=-160;
- C=add_seq2aln (CL,C, CL->S);
- print_aln (C);
-
- myexit (EXIT_FAILURE);
-
- while ( step>0)
- {
- C=approximate_domain (min_start,max_start,step,min_len,max_len, step,&start, &len, &score, CL);
- min_start=start-step;
- max_start=start+step;
- min_len=len-step;
- max_len=len+step;
- step=step/2;
- }
-
- C=extract_domain_with_coordinates (C,start-10, len+20,CL);
- C->output_res_num=1;
- print_aln (C);
-
- myexit (EXIT_FAILURE);
- return C;
-
-
- }
-Alignment * interactive_domain_extraction ( Constraint_list *CL)
- {
- int LEN=0;
- int START=1;
- int SCALE=2;
- int GOPP=3;
-
- int iteration=0;
- char *choice;
- int a,b, c;
- int index;
- char *s;
- char last_start[100];
- char out_format[100];
- Alignment *RESULT=NULL;
- Alignment *PREVIOUS=NULL;
- Alignment *C=NULL;
- Alignment *EA=NULL;
-
- int **parameters;
-
-
- choice=vcalloc ( 100, sizeof (char));
- parameters=declare_int (10000, 4);
-
- parameters[0][START]=(CL->moca)->moca_start;
- parameters[0][LEN]= (CL->moca)->moca_len;
- parameters[0][SCALE]=(CL->moca)->moca_scale;
- parameters[0][GOPP]=CL->gop;
- iteration=0;
- sprintf ( last_start, "%d", (CL->moca)->moca_start);
- sprintf ( out_format, "mocca_aln");
-
- print_moca_interactive_choices ();
- while ( !strm4 (choice, "Q","X", "q", "x" ))
- {
- c=choice[0];
-
- if (c=='b' || c=='B')
- {
- iteration-=atoi(choice+1)+1;
-
- if (iteration<0)iteration=1;
- }
- else
- {
- iteration++;
- parameters[iteration][START]=parameters[iteration-1][START];
- parameters[iteration][LEN]=parameters[iteration-1][LEN];
- parameters[iteration][SCALE]=parameters[iteration-1][SCALE];
- parameters[iteration][GOPP]=parameters[iteration-1][GOPP];
-
- if ( c=='>')parameters[iteration][LEN]=atoi(choice+1);
- else if ( c=='|')
- {
- sprintf ( last_start, "%s", choice);
- parameters[iteration][START]=0;
- s=strrchr(choice, ':');
-
- if (s==NULL)
- {
- parameters[iteration][START]=atoi(choice+1);
- }
- else
- {
-
- s[0]='\0';
-
- if((index=name_is_in_list (choice+1,(CL->S)->name,(CL->S)->nseq,100))==-1)
- {
- fprintf ( stderr, "\n\tERROR: %s NOT in Sequence Set",choice+1);
- continue;
- }
-
- for ( a=0; a< index; a++)
- {
- parameters[iteration][START]+=(CL->S)->len[a]+1;
- }
- parameters[iteration][START]+=atoi(s+1)-1;
- }
-
- }
- else if ( c=='C'||c=='c')parameters[iteration][SCALE]=atoi(choice+1);
- else if ( c=='G'||c=='g')
- {
- parameters[iteration][GOPP]=atoi(choice+1);
- CL->gop=parameters[iteration][GOPP];
- }
- else if ( c=='F'||c=='f')
- {
- sprintf ( out_format, "%s", choice+1);
- }
- else if ( c=='S'||c=='s')
- {
- if (choice[1]=='\0')sprintf ( choice, "default.domain_aln.%d", iteration);
- output_format_aln (out_format,RESULT,EA=fast_coffee_evaluate_output(RESULT, CL),choice+1);
- fprintf (stderr, "\tOutput file [%15s] in [%10s] format\n",choice+1,out_format);
- free_aln (EA);
- }
- else if (c=='\0')
- {
- if ( parameters[iteration][SCALE]>0)
- {
- fprintf ( stderr, "\nWARNING: THRESHOLD RESET to 0");
- parameters[iteration][SCALE]=0;
- }
-
- (CL->moca)->moca_scale=parameters[iteration][SCALE];
- CL->gop=parameters[iteration][GOPP];
-
- C=extract_domain_with_coordinates (C,parameters[iteration][START],parameters[iteration][LEN],CL);
-
- if ( C==NULL)
- {
- fprintf ( stderr, "\nERROR: ILLEGAL COORDINATES! SEQUENCE BOUNDARY CROSSED\n");
- for ( b=1,a=0; a< (CL->S)->nseq-1; a++)
- {
-
- fprintf ( stderr, "\n\t%15s=> Abs:[%d %d] Rel:[0 %d]", (CL->S)->name[a],b, b+(CL->S)->len[a]-1,(CL->S)->len[a]);
- b+=(CL->S)->len[a];
- }
- fprintf ( stderr, "\n");
- }
- else if (parameters[iteration][START]==0 && parameters[iteration][LEN]==0)
- {
- fprintf ( stderr, "\n\tEnter the following parameters:\n\n\t\tSTART value: |x [Return]\n\t\tLENgth value: >y [Return]\n\t\ttype [Return]\n\n");
- fprintf ( stderr, "\n\n\tSTART is measured on the total length of the concatenated sequences\n\tx and y are positive integers\n\n");
- }
-
- else if ( C->nseq==0)
- {
- fprintf ( stderr, "\nNO MATCH FOUND: LOWER THE SCALE (C)\n");
- }
- else
- {
- RESULT=copy_aln ( C, RESULT);
- unpack_seq_aln (RESULT, CL);
- RESULT->output_res_num=1;
-
- output_format_aln (out_format,RESULT,EA=fast_coffee_evaluate_output(RESULT, CL),"stdout");
- free_aln(EA);
- PREVIOUS=copy_aln ( RESULT, PREVIOUS);
- free_aln (C);
- print_moca_interactive_choices ();
-
- }
- }
-
- fprintf ( stderr, "\t[ITERATION %3d][START=%s][LEN=%3d][GOPP=%3d][SCALE=%4d]\t",iteration,last_start,parameters[iteration][LEN],parameters[iteration][GOPP],parameters[iteration][SCALE]);
- a=0;
- fprintf ( stderr, "Your Choice: ");
- while ( (c=fgetc(stdin))!='\n')choice[a++]=c;
- choice[a]=0;
- }
- }
-
- if (!RESULT)myexit(EXIT_SUCCESS);
- if ( RESULT)RESULT->output_res_num=0;
- return RESULT;
- }
-
-int print_moca_interactive_choices ()
-{
- fprintf ( stderr, "\n**************************************************************");
- fprintf ( stderr, "\n******************** MOCCA: %s ***********",VERSION);
- fprintf ( stderr, "\n**************************************************************");
-
-
-fprintf ( stderr, "\nMENU: Type Flag[number] and Return: ex |10");
-
-fprintf ( stderr, "\n\t|x -->Set the START to x");
-fprintf ( stderr, "\n\t 100 start=100 on concatenated sequences");
-fprintf ( stderr, "\n\t human:100 start=100 on human sequence");
-fprintf ( stderr, "\n\t>x -->Set the LEN to x");
-fprintf ( stderr, "\n\tGx -->Set the Gap Opening Penalty to x");
-fprintf ( stderr, "\n\tCx -->Set the sCale to x");
-fprintf ( stderr, "\n\tSname -->Save the Alignment ");
-fprintf ( stderr, "\n\tFformat -->Save the Alignment Format");
-
-fprintf ( stderr, "\n\treturn -->Compute the Alignment");
-
-fprintf ( stderr, "\n\tX -->eXit\n\n");
-
-return 0;
-}
-
-Alignment * approximate_domain ( int min_start, int max_start, int step_start,int min_len, int max_len, int step_len, int *best_start, int *best_len, int *best_score, Constraint_list *CL)
- {
- Alignment *C=NULL;
- int start;
- int len;
- int score;
-
- /*1 Extract the first*/
- best_score[0]=UNDEFINED;
- best_start[0]=min_start;
- best_len[0]=min_len;
-
- for (start=min_start; start< max_start; start+=step_start)
- {
- for ( len=min_len; len<max_len; len+=step_len)
- {
- C=extract_domain_with_coordinates (C,start,len,CL);
- if ( C==NULL)continue;
- score=((CL->moca)->evaluate_domain)(C, CL);
- fprintf ( stderr, "\nSTART=%d LEN=%3d SCORE=%5d [%d]",start,len,score, C->nseq);
-
-
- if ( best_score[0]==UNDEFINED)best_score[0]=score;
- if ( score>best_score[0])
- {
- best_score[0]=score;
- best_start[0]=start;
- best_len[0]=len;
- }
- }
- }
-
- C=extract_domain_with_coordinates (C,best_start[0], best_len[0],CL);
- C->output_res_num=1;
- return C;
- }
-int measure_domain_length ( Constraint_list *CL,Alignment *IN, int start, int min_len, int max_len, int step)
- {
- Alignment *C=NULL;
- int score, best_score,best_len,a, b, l;
- int *score_matrix, *len_matrix;
- int n_val, best_val;
-
- score_matrix=vcalloc ( max_len, sizeof (int));
- len_matrix=vcalloc ( max_len, sizeof (int));
-
-
- l=strlen ( (CL->S)->seq[0]);
-
- min_len=MAX(0, min_len);
- min_len=MIN(l-start, min_len);
-
- if ( !IN)C=extract_domain_with_coordinates (C,start,min_len, CL);
- else
- {
- C=copy_aln (IN, C);
- C->len_aln=min_len;
- for ( a=0; a< C->nseq; a++)C->seq_al[a][min_len]='\0';
- C=add_seq2aln (CL,C, CL->S);
- }
-
- best_score= score=((CL->moca)->evaluate_domain)(C, CL);
-
-
- min_len=MAX(0, min_len);
- for ( best_len=best_val=n_val=0,b=min_len; b<max_len && (start+b)<l; b+=step, n_val++)
- {
- if ( !IN)C=extract_domain_with_coordinates (C,start, b, CL);
- else
- {
- C=copy_aln (IN, C);
- C->len_aln=min_len;
- for ( a=0; a< C->nseq; a++)C->seq_al[a][b]='\0';
- C=add_seq2aln (CL,C, CL->S);
- }
- if ( C->len_aln>0 )score=((CL->moca)->evaluate_domain)(C, CL);
- else score=-1;
-
- if ( score< -3000)break;
-
- fprintf ( stderr, "\n\t%d %d=>%d (%d, %d)[%d]",start, b, score, C->nseq, C->len_aln, step);
- score_matrix[n_val]=score;
- len_matrix [n_val]=b;
- if ( score>best_score)
- {
- best_score=score;
- best_len=b;
- best_val=n_val;
- }
- }
- free_aln(C);
-
- for ( a=best_val; a<n_val;a++)
- {
- if (score_matrix[a]>best_score/2)best_len=len_matrix[a];
- else break;
- }
- vfree ( score_matrix);
- vfree ( len_matrix);
-
- return best_len;
- }
-
-Alignment *extract_domain_with_coordinates ( Alignment *RESULT,int start, int len, Constraint_list *CL)
-{
- int a;
- char *buf;
- Alignment *SEQ_DOMAIN=NULL;
-
-
-
-
-
- /*ADJUST THE DIRECTION OF THE DOMAIN: len<0:left and len>0:right*/
-
- if (len>0);
- else if (len<0)
- {
- len=len*-1;
- start=start-len+1;
- }
-
- /*CHECK THAT THE BOUNDARY CONDITIONS*/
-
-
- if (start<0 || (!CL->packed_seq_lu && (start+len)>strlen((CL->S)->seq[0])) ||(CL->packed_seq_lu && (start+len)>strlen((CL->S)->seq[(CL->S)->nseq-1])) )return NULL;
- else
- {
- for ( a=start; a< start+len; a++)
- {
- if ((CL->moca)->forbiden_residues && (CL->moca)->forbiden_residues[0][a+1]==UNDEFINED)
- {
- fprintf ( stderr, "*");
- return NULL;
- }
- }
- }
-
- /*EXTRACT THE DOMAIN*/
-
- SEQ_DOMAIN=add_seq2aln (CL,SEQ_DOMAIN, CL->S);
- buf=extract_char (SEQ_DOMAIN->seq_al[0], start, len);
-
- for (a=0; a<len; a++)if ( buf[a]=='X'){free_aln(SEQ_DOMAIN);return NULL;}
- sprintf ( SEQ_DOMAIN->seq_al[0], "%s", buf);
- SEQ_DOMAIN->order[0][1]=start;
- SEQ_DOMAIN=add_seq2aln (CL,SEQ_DOMAIN, CL->S);
-
-
-
- return SEQ_DOMAIN;
-}
-
-
-
-
-int get_starting_point ( Constraint_list *CL)
-{
- int a;
-
-
- int l;
-
-
-
-
-
- int **seq;
- int start;
- int *entry=NULL;
-
- l=strlen ( (CL->S)->seq[0]);
-
- seq=declare_int ( l, 2);
-
-
-
- while (entry=extract_entry (CL))
- {
- seq[entry[R1]][1]=entry[R1];
- seq[entry[R2]][1]=entry[R2];
- if ((CL->moca) && (CL->moca)->forbiden_residues && ((CL->moca)->forbiden_residues[0][entry[R1]]==UNDEFINED||(CL->moca)->forbiden_residues[0][entry[R2]]==UNDEFINED ))continue;
- else
- {
- seq[entry[R1]][0]+=entry[MISC];
- seq[entry[R2]][0]+=entry[MISC];
- }
- }
-
- sort_int_inv ( seq, 2, 0, 0, l-1);
- fprintf ( stderr, "\nStart=%d %d", seq[0][1], seq[0][0]);
- start=seq[0][1];
-
-
- free_int ( seq, -1);
- return start;
-
-
-}
-
-
-int * analyse_sequence ( Constraint_list *CL)
-{
- int a, p;
- int len, start, n_dots;
- int left, right, tw, r, w;
- int best_tw, best_start=0, best_len=0;
- int l;
- int max_len=200;
-
-
- l=strlen (( CL->S)->seq[0]);
-
- for ( best_tw=UNDEFINED,start=0; start<l; start++)
- {
- for ( len=10; len< max_len && start+len<l; len++)
- {
- left=start;
- right=start+len;
- for (tw=0, p=0; p<len; p++)
- {
- n_dots=CL->residue_index[0][p+start+1][0];
-
- for ( a=1; a<n_dots; a+=ICHUNK)
- {
-
- r=CL->residue_index[0][p+start+1][a+R2];
- w=CL->residue_index[0][p+start+1][a+WE];
-
- if (r<left || r>right)tw+=w;
- }
- }
-
- if ( tw> best_tw || best_tw==UNDEFINED)
- {
- best_tw=tw;
- best_start=start;
- best_len=len;
- }
- }
- }
- fprintf ( stderr, "\nStart=%d Len=%d", best_start, best_len);
- return NULL;
-}
-
-/*********************************COPYRIGHT NOTICE**********************************/
-/*© Centre National de la Recherche Scientifique (CNRS) */
-/*and */
-/*Cedric Notredame */
-/*Fri Aug 8 19:03:27 MDT 2003. */
-/*All rights reserved.*/
-/*NOTICE: |*/
-/* This file is an integral part of the */
-/* T-COFFEE Software. */
-/* Its content is protected and all */
-/* the conditions mentioned in the licensing */
-/* agreement of the software apply to this file.*/
-/*............................................... |*/
-/* If you need some more information, or if you */
-/* wish to obtain a full license, please contact: */
-/* cedric.notredame@europe.com*/
-/*............................................... |*/
-/**/
-/**/
-/* */
-/*********************************COPYRIGHT NOTICE**********************************/
-/******************************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