+++ /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 count_threshold_nodes (Alignment *A, NT_node P, int t);
-int set_node_score (Alignment *A, NT_node P, char *mode);
-char *split_nodes_nseq (Alignment *A, NT_node P, int max_nseq, char *list);
-char *split_nodes_idmax (Alignment *A, NT_node P, int max_id, char *list);
-/******************************************************************/
-/* MAIN DRIVER */
-/* */
-/* */
-/******************************************************************/
-int method_uses_structure(TC_method *M)
-{
- if ( strchr (M->seq_type, 'P'))return 1;
- else return 0;
-}
-Constraint_list *profile2list (Job_TC *job, int nprf)
-{
- int a,b, s1, s2, r1, r2,w2;
- Constraint_list *CL, *RCL;
- Sequence *S;
- Alignment *A1, *A2;
- int *iA1, *iA2;
- TC_method *M;
- int ***RI, ***NRI;
- Sequence *RS;
- int Rne;
-
- int *seqlist;
- int **cache;
- static int *entry;
-
-
- if (!entry)entry=vcalloc ( ICHUNK+3, sizeof (int));
-
- CL=(job->io)->CL;
- //1- buffer CL
-
- RS=CL->S;
- RI=CL->residue_index;
- Rne=CL->ne;
- M=(job->param)->TCM;
-
-
-
-
-
-
- //Index
- seqlist=string2num_list ((job->param)->seq_c)+1;
- A1=seq2profile(CL->S, seqlist[1]);
- A2=seq2profile(CL->S, seqlist[2]);
-
- S=merge_seq (A1->S, NULL );
- S=merge_seq (A2->S, S);
-
- iA1=get_name_index (A1->name,A1->nseq, S->name, S->nseq);
- iA2=get_name_index (A2->name,A2->nseq, S->name, S->nseq);
- cache=vcalloc ( S->nseq, sizeof (int*));
- for (a=0; a<A1->nseq; a++)cache[iA1[a]]=seq2inv_pos (A1->seq_al[a]);
- for (a=0; a<A2->nseq; a++)cache[iA2[a]]=seq2inv_pos (A2->seq_al[a]);
-
- //Declare local CL
- CL->S=S;
- CL->residue_index=declare_residue_index(S);
- CL->ne=0;
-
- //Compute lib
- for (a=0; a<A1->nseq; a++)
- for ( b=0; b<A2->nseq; b++)
- {
- char buf[1000];
- sprintf ( buf, "2 %d %d",iA1[a], iA2[b]);
- job=print_lib_job (job, "param->seq_c=%s", buf);
- CL=seq2list (job);
- }
-
- //restaure CL;
- CL->S=RS;
- NRI=CL->residue_index;
- CL->residue_index=RI;
- CL->ne=Rne;
-
-
-
- //incorporate new lib
- entry[SEQ1]=seqlist[1];
- entry[SEQ2]=seqlist[2];
- for (a=0; a<A1->nseq; a++)
- {
- s1=iA1[a];
-
- for (r1=1; r1<=S->len[s1];r1++)
- {
- for (b=1; b<NRI[s1][r1][0]; b+=ICHUNK)
- {
- int s2=NRI[s1][r1][b+SEQ2];
- int r2=NRI[s1][r1][b+R2];
- int w2=NRI[s1][r1][b+WE];
-
- entry[R1]=cache[s1][r1];
- entry[R2]=cache[s2][r2];
- entry[WE]=w2;
- entry[CONS]=1;
- entry[MISC]=1;
- add_entry2list (entry,CL);
- }
- }
- }
-
- free_int (cache, -1);
- free_arrayN (NRI, 3);
- vfree (entry);
- free_sequence (S, S->nseq);
- return (job->io)->CL=CL;
-}
-
-Constraint_list *seq2list ( Job_TC *job)
- {
- char *mode;
- Alignment *A=NULL;
- Constraint_list *PW_CL;
- Constraint_list *RCL=NULL;
-
- int full_prf, nprf;
- int *seqlist;
-
-
- Sequence *S, *STL;
- Constraint_list *CL;
-
-
- char *seq;
- char *weight;
- TC_method *M;
-
-
- M=(job->param)->TCM;
- mode=M->executable;
- weight=M->weight;
-
- PW_CL=M->PW_CL;
-
- CL=(job->io)->CL;
- seq=(job->param)->seq_c;
-
-
-
- S=(CL)?CL->S:NULL;
- STL=(CL)?CL->STRUC_LIST:NULL;
-
- seqlist=string2num_list (seq)+1;
-
-
-
-/*Proteins*/
-
-
- if ( strncmp (CL->profile_comparison, "full", 4)==0)
- {
- full_prf=1;
- if ( CL->profile_comparison[4])nprf=atoi ( CL->profile_comparison+4);
- else
- nprf=0;
- }
- else
- {
- full_prf=0;
- }
-
- if ((method_uses_structure (M)) && profile2P_template_file (CL->S, seqlist[1]) && profile2P_template_file (CL->S, seqlist[2]))
- {
- RCL=profile2list (job, nprf);
- }
- else if ( strm (mode, "ktup_msa"))
- {
- RCL=hasch2constraint_list (CL->S, CL);
- }
- else if ( strm (mode, "test_pair") || strm ( mode,"fast_pair") || strm (mode, "ifast_pair") \
- || strm ( mode, "diag_fast_pair")|| strm (mode, "idiag_fast_pair")\
- || strm ( mode, "blast_pair") || strm (mode, "lalign_blast_pair") \
- || strm ( mode, "viterbi_pair") || strm (mode, "slow_pair") || strm(mode, "glocal_pair") || strm (mode, "biphasic_pair") \
- || strm ( mode, "islow_pair") || strm (mode, "tm_slow_pair") || strm (mode, "r_slow_pair") \
- || strm ( mode, "lalign_id_pair")|| strm (mode, "tm_lalign_id_pair") || strm (mode , "lalign_len_pair") \
- || strm (mode, "prrp_aln") || strm ( mode, "test_pair") \
- || strm (mode, "cdna_fast_pair") || strm (mode, "diaa_slow_pair") || strm (mode, "monoaa_slow_pair")\
- || strncmp (mode,"cdna_fast_pair",14)==0 \
- )
- {
-
- A=fast_pair (job);
- RCL=aln2constraint_list ((A->A)?A->A:A, CL,weight);
- }
-
- else if ( strm ( mode, "subop1_pair") || strm ( mode, "subop2_pair") )
- {
- A=fast_pair (job);
- RCL=A->CL;
- }
- else if ( strm ( mode, "proba_pair") )
- {
- A=fast_pair (job);
- RCL=A->CL;
- }
- else if ( strm ( mode, "best_pair4prot"))
- {
- RCL=best_pair4prot (job);
- }
- else if ( strm ( mode, "best_pair4rna"))
- {
- RCL=best_pair4rna (job);
- }
- else if ( strm ( mode, "exon2_pair"))
- {
- char weight2[1000];
-
- A=fast_pair (job);
- sprintf ( weight2, "%s_subset_objOBJ-",weight);
- RCL=aln2constraint_list (A, CL,weight2);
- }
- else if ( strm ( mode, "exon_pair"))
- {
- A=fast_pair (job);
- RCL=aln2constraint_list (A, CL,weight);
-
- }
- else if ( strm ( mode, "exon3_pair"))
- {
- char weight2[1000];
-
- A=fast_pair (job);
- sprintf ( weight2, "%s_subset_objOBJ-",weight);
- RCL=aln2constraint_list (A, CL,weight2);
- }
-
-/*STRUCTURAL METHODS*/
-
- else if ( strm (mode, "seq_msa"))
- {
- RCL=seq_msa(M, seq, CL);
- }
-/*STRUCTURAL METHODS*/
- else if (strm (mode, "profile_pair") || strm (mode, "hh_pair"))
- {
- RCL=profile_pair (M, seq, CL);
- }
-
- else if ( strm (mode, "sap_pair"))
- {
- RCL=sap_pair (seq, weight, CL);
- }
- else if ( strm (mode, "thread_pair"))
- {
- RCL=thread_pair (M,seq, CL);
- }
- else if ( strm (mode, "pdb_pair"))
- {
- RCL=pdb_pair (M,seq, CL);
- }
- else if (strm (mode, "rna_pair"))
- {
- RCL=rna_pair(M, seq, CL);
- }
- else if ( strm (mode, "pdbid_pair"))
- {
- RCL=pdbid_pair (M,seq, CL);
- }
- else if ( strm (mode, "fugue_pair"))
- {
- RCL=thread_pair (M,seq, CL);
- }
- else if ( strm (mode, "lsqman_pair"))
- {
- RCL=lsqman_pair(seq, CL);
- }
- else if ( strm ( mode, "align_pdb_pair"))
- {
- RCL=align_pdb_pair ( seq,"gotoh_pair_wise", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
- }
- else if ( strm ( mode, "lalign_pdb_pair"))
- {
- RCL=align_pdb_pair ( seq,"sim_pair_wise_lalign", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
- }
- else if ( strm ( mode, "align_pdb_pair_2"))
- {
- RCL=align_pdb_pair_2 ( seq, CL);
- }
- else
- {
- fprintf ( CL->local_stderr, "\nERROR: THE FUNCTION %s DOES NOT EXIST [FATAL:%s]\n", mode, PROGRAM);crash("");
- }
- add_method_output2method_log (NULL,NULL, (A&&A->len_aln)?A:NULL,RCL, NULL);
- RCL=(RCL==NULL)?CL:RCL;
-
-
- vfree ( seqlist-1);
- free_aln (A);
- return RCL;
- }
-
-Constraint_list *method2pw_cl (TC_method *M, Constraint_list *CL)
- {
- char *mode;
- Constraint_list *PW_CL=NULL;
- Sequence *S;
- char mat[100], *m;
- char group_mat[100];
-
-
-
- mode=M->executable;
- PW_CL=copy_constraint_list ( CL, SOFT_COPY);
- PW_CL->pw_parameters_set=1;
-
-
-
- S=(PW_CL)?PW_CL->S:NULL;
-
- /*DNA or Protein*/
- m=PW_CL->method_matrix;
- if ( strm ((PW_CL->S)->type, "PROTEIN"))
- {
-
- sprintf ( mat, "%s", (strm(m, "default"))?"blosum62mt":m);
- sprintf (group_mat, "vasiliky");
- PW_CL->ktup=2;
-
- }
- else if ( strm ((PW_CL->S)->type, "DNA") || strm ((PW_CL->S)->type, "RNA") )
- {
-
- sprintf(group_mat, "idmat");
- sprintf ( mat, "%s", (strm(m, "default"))?"dna_idmat":m);
- PW_CL->ktup=5;
-
- }
- if ( M->matrix[0])sprintf ( mat, "%s", M->matrix);
-
- PW_CL->M=read_matrice (mat);
-
-
- if ( M->gop!=UNDEFINED) {PW_CL->gop=M->gop;}
- else
- {
- PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
- }
-
- if ( M->gep!=UNDEFINED)PW_CL->gep=M->gep;
- else PW_CL->gep=-1;
-
- if (M->extend_seq ==1) PW_CL->extend_seq=1;
- if (M->reverse_seq==1)PW_CL->reverse_seq=1;
-
-
- if ( strm2 ( mode,"fast_pair", "ifast_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- if ( !PW_CL->use_fragments)PW_CL->diagonal_threshold=0;
- else PW_CL->diagonal_threshold=6;
-
- sprintf (PW_CL->dp_mode, "fasta_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
-
- if ( strm ( mode, "fast_pair"))
- {
- PW_CL->residue_index=NULL;
-
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- }
- else if ( strm2 ( mode,"diag_fast_pair","idiag_fast_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->S=CL->S;
-
- PW_CL->use_fragments=1;
- PW_CL->diagonal_threshold=3;
-
- sprintf (PW_CL->dp_mode, "fasta_pair_wise");
- PW_CL->ktup=1;
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
-
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
-
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode,"blast_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
-
- PW_CL->use_fragments=0;
-
- PW_CL->pair_wise=gotoh_pair_wise;
- PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode,"lalign_blast_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
-
- PW_CL->use_fragments=0;
- PW_CL->pair_wise=sim_pair_wise_lalign;
- PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
- PW_CL->lalign_n_top=10;
-
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode,"viterbi_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "viterbi_pair_wise");
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode,"glocal_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "glocal_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
-
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode,"test_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "test_pair_wise");
- }
- else if ( strm ( mode,"sticky_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=cw_profile_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp_sticky");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- }
-
- else if ( strm ( mode,"slow_pair")|| strm (mode, "islow_pair" ) )
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
-
- if ( strm ( "islow_pair", mode))
- {
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=residue_pair_extended_list;
- PW_CL->extend_jit=1;
- }
- else if ( strm ("slow_pair", mode) )
- {
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=cw_profile_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- }
- else if ( strm (mode, "subop1_pair"))
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "subop1_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm (mode, "biphasic_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "biphasic_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=cw_profile_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm (mode, "proba_pair"))
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "proba_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
-
- else if ( strm (mode, "diaa_slow_pair"))
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_diaa_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm (mode, "r_slow_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- PW_CL->reverse_seq=1;
- }
- else if ( strm (mode, "tm_slow_pair"))
- {
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_tm_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm (mode, "monoaa_slow_pair"))
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_monoaa_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm (mode, "subop2_pair"))
- {
-
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "subop2_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
-
- else if (strm ( mode, "exon2_pair"))
- {
- int a;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
-
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
-
- for ( a=0; a<60; a++)
- {
- PW_CL->M['x'-'A'][a]=0;
- PW_CL->M[a]['x'-'A']=0;
- PW_CL->M['X'-'A'][a]=0;
- PW_CL->M[a]['X'-'A']=0;
- }
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if (strm ( mode, "exon3_pair"))
- {
- int a;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
-
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->residue_index=NULL;
-
- for ( a=0; a<60; a++)
- {
- PW_CL->M['x'-'A'][a]=0;
- PW_CL->M[a]['x'-'A']=0;
- PW_CL->M['X'-'A'][a]=0;
- PW_CL->M[a]['X'-'A']=0;
- }
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if (strm ( mode, "exon_pair"))
- {
- int a;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
-
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
- sprintf (PW_CL->matrix_for_aa_group, "%s",group_mat);
- PW_CL->residue_index=NULL;
-
- for ( a=0; a<60; a++)
- {
- PW_CL->M['x'-'A'][a]=0;
- PW_CL->M[a]['x'-'A']=0;
- PW_CL->M['X'-'A'][a]=0;
- PW_CL->M[a]['X'-'A']=0;
- }
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode , "lalign_len_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- PW_CL->pair_wise=sim_pair_wise_lalign;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->lalign_n_top=CL->lalign_n_top;
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode , "lalign_id_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- PW_CL->pair_wise=sim_pair_wise_lalign;
- PW_CL->evaluate_residue_pair=evaluate_matrix_score;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->lalign_n_top=CL->lalign_n_top;
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->extend_jit=0;
- }
- else if ( strm ( mode , "tm_lalign_id_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- PW_CL->pair_wise=sim_pair_wise_lalign;
- PW_CL->evaluate_residue_pair=evaluate_tm_matrix_score;
- PW_CL->get_dp_cost=slow_get_dp_cost;
- PW_CL->lalign_n_top=CL->lalign_n_top;
- sprintf (PW_CL->matrix_for_aa_group,"%s", group_mat);
- PW_CL->extend_jit=0;
- }
-/*CDNA*/
- else if ( strm ( mode, "cdna_cfast_pair"))
- {
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->S=CL->S;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "cfasta_cdna_pair_wise");
-
- PW_CL->M=read_matrice (strcpy ( mat, "blosum62mt"));
- PW_CL->extend_jit=0;
-
- PW_CL->f_gop=CL->f_gop;
- PW_CL->f_gep=CL->f_gep;
- PW_CL->get_dp_cost=get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
- PW_CL->ktup=1;
- }
- else if ( strm ( mode, "cdna_fast_pair") || strncmp (mode,"cdna_fast_pair",14)==0)
- {
-
- PW_CL->residue_index=NULL;
- PW_CL->maximise=1;
- PW_CL->TG_MODE=1;
- PW_CL->use_fragments=0;
- sprintf (PW_CL->dp_mode, "fasta_cdna_pair_wise");
-
- PW_CL->extend_jit=0;
- PW_CL->gop=-5;
- PW_CL->gep=-1;
- PW_CL->f_gop=-15;
- PW_CL->f_gep=0;
-
- PW_CL->get_dp_cost=get_dp_cost;
- PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
- PW_CL->ktup=1;
- }
- else
- {
- free_constraint_list (PW_CL);
- PW_CL=NULL;
- }
-
-
- if (!strm (CL->method_evaluate_mode, "default"))
- {
- choose_extension_mode (CL->method_evaluate_mode, PW_CL);
- }
- return PW_CL;
- }
-/******************************************************************/
-/* MULTIPLE ALIGNMENTS */
-/* */
-/* */
-/******************************************************************/
-Alignment * compute_prrp_aln (Alignment *A, Constraint_list *CL)
- {
- char *tmpseq=NULL;
- char *tmpaln=NULL;
- char command[10000];
- Sequence *S;
-
- tmpseq=vtmpnam(NULL);
- tmpaln=vtmpnam(NULL);
-
-
- A=seq2aln ( CL->S, A, 1);
- output_gotoh_seq (tmpseq, A);
- sprintf ( command, "prrp -E/dev/null -o%s -F9 %s >/dev/null", tmpaln, tmpseq);
- my_system (command);
- if (!check_file_exists(tmpaln)){return NULL;}
- S=get_fasta_sequence (tmpaln, NULL);
-
- S->contains_gap=0;
- A=seq2aln(S, A,0);
- free_sequence (S, S->nseq);
-
- return A;
- }
-
-Alignment *seq2clustalw_aln (Sequence *S)
-{
- return aln2clustalw_aln (seq2aln (S, NULL,RM_GAP), NULL);
-}
-
-Alignment * aln2clustalw_aln (Alignment *B, Constraint_list *CL)
-{
- char *seq=NULL,*aln=NULL, command[1000];
-
- output_fasta_seq (seq=vtmpnam (NULL), B);
- sprintf ( command, "clustalw -infile=%s -outorder=input -outfile=%s %s", seq, aln=vtmpnam (NULL), TO_NULL_DEVICE);
- my_system (command);
-
- if (!check_file_exists(aln))
- return NULL;
- else{B->nseq=0;return main_read_aln(aln,B);}
-}
-Alignment * compute_tcoffee_aln_quick (Alignment *A, Constraint_list *CL)
- {
- char *tmpseq=NULL;
- char *tmpaln=NULL;
- char command[10000];
-
-
- tmpseq=vtmpnam(NULL);
- tmpaln=vtmpnam(NULL);
-
- if ( CL)A=seq2aln ( CL->S, A, 1);
- output_fasta_seq (tmpseq, A);
-
- sprintf ( command, "t_coffee -seq %s -very_fast -outfile %s -quiet ",tmpseq,tmpaln);
-
- my_system (command);
- if (!check_file_exists(tmpaln))return NULL;
- A->nseq=0;
- A=main_read_aln(tmpaln,A);
-
- vremove( tmpseq);
- vremove (tmpaln);
- return A;
- }
-
-Alignment * compute_clustalw_aln (Alignment *A, Constraint_list *CL)
- {
- char *tmpseq=NULL;
- char *tmpaln=NULL;
- char command[10000];
-
-
- tmpseq=vtmpnam(NULL);
- tmpaln=vtmpnam(NULL);
-
- A=seq2aln ( CL->S, A, 1);
- output_fasta_seq (tmpseq, A);
-
- sprintf ( command, "clustalw %sinfile=%s %soutfile=%s %s",CWF, tmpseq,CWF, tmpaln,TO_NULL_DEVICE);
-
- my_system (command);
- if (!check_file_exists(tmpaln))return NULL;
- A->nseq=0;
- A=main_read_aln(tmpaln,A);
-
- vremove( tmpseq);
- vremove (tmpaln);
- return A;
- }
-
-Alignment * realign_block ( Alignment *A, int col1, int col2, char *pg)
-{
- /*Uses pg: (pg -infile=<infile> -outfile=<outfile> to realign the block [col1 col2[
- Only guaranteed if pg can handle empty sequences
- set pg to NULL to use the default program
- */
-
-
- Alignment *L, *M, *R;
- char *seq_name;
- char *aln_name;
- char command[1000], script[1000];
-
-
-
- seq_name=vtmpnam(NULL);
- aln_name=vtmpnam (NULL);
-
- L=copy_aln (A, NULL);
- M=copy_aln (A, NULL);
- R=copy_aln (A, NULL);
-
- L=extract_aln ( L, 0, col1);
- M=extract_aln ( M, col1, col2);
- R=extract_aln ( R, col2, A->len_aln);
- output_fasta_seq (seq_name, M);
-
- sprintf ( script, "%s", (pg==NULL)?"t_coffee":pg);
-
- sprintf ( command, "%s -infile=%s -outfile=%s %s", script,seq_name, aln_name, TO_NULL_DEVICE);
- my_system ( command);
- free_aln (M);
- M=main_read_aln (aln_name, NULL);
-
-
- M=reorder_aln (M, L->name,L->nseq);
- L=aln_cat (L, M);
- L=aln_cat (L, R);
- A=copy_aln (L, A);
- free_aln (L);free_aln (M); free_aln (R);
-
- return A;
-}
-
-
-
-
-
-/******************************************************************/
-/* DNA */
-/* */
-/* */
-/******************************************************************/
-
-
-/******************************************************************/
-/* STRUCTURES */
-/* */
-/* */
-/******************************************************************/
-
-
-
-
-
-Constraint_list * align_pdb_pair_2 (char *seq, Constraint_list *CL)
- {
- char *tmp_name=NULL;
- int s1, s2;
-
-
- static char *command;
- static char *program;
-
- tmp_name=vtmpnam ( NULL);
-
- if ( !program)program=vcalloc ( LONG_STRING, sizeof (char));
- if ( !command)command=vcalloc ( LONG_STRING, sizeof (char));
-
-#ifndef ALIGN_PDB_4_TCOFFEE
- if ( getenv ( "ALIGN_PDB_4_TCOFFEE")==NULL)crash ("ALIGN_PDB_4_TCOFFEE IS NOT DEFINED");
- else sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
-#else
- if ( getenv ( "ALIGN_4_TCOFFEE")==NULL)sprintf (program, "%s", ALIGN_PDB_4_TCOFFEE);
- else sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
-#endif
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- sprintf ( command , "%s -in P%s P%s -gapopen=-40 -max_delta=2.5 -gapext=0 -scale=0 -hasch_mode=hasch_ca_trace_bubble -maximum_distance=10 -output pdb_constraint_list -outfile stdout> %s%s",program, (CL->S)->file[s1], (CL->S)->file[s2], get_cache_dir(),tmp_name);
-
- my_system ( command);
- CL=read_constraint_list_file(CL, tmp_name);
-
-
- vremove ( tmp_name);
-
-
- return CL;
- }
-
-Constraint_list *align_pdb_pair (char *seq_in, char *dp_mode,char *evaluate_mode, char *file, Constraint_list *CL, Job_TC *job)
- {
- int s1, s2;
- char seq[1000];
- char name1[1000];
- char name2[1000];
-
-
- Constraint_list *PWCL;
- Alignment *F;
-
- sprintf ( seq, "%s",seq_in);
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
-
-
- sprintf (name1, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s1], (CL->S)->name[s2], dp_mode);
- sprintf (name2, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s2], (CL->S)->name[s1], dp_mode);
-
-
- if ( check_file_exists (name1) && is_lib(name1))CL=read_constraint_list_file(CL,name1);
- else if ( check_file_exists (name2) && is_lib(name2))CL=read_constraint_list_file(CL,name2);
- else
- {
- PWCL=set_constraint_list4align_pdb ( CL,s1,dp_mode, evaluate_mode, NULL);
- PWCL=set_constraint_list4align_pdb ( CL,s2,dp_mode, evaluate_mode, NULL);
- ((job->param)->TCM)->PW_CL=PWCL;
- F=fast_pair (job);
- output_constraints (name1, "100", F);
- CL=aln2constraint_list (F, CL, "100");
- free_aln (F);
- }
- return CL;
- }
-
-Constraint_list * profile_pair (TC_method *M , char *in_seq, Constraint_list *CL)
- {
-
- char seq[1000];
- int a, s1, s2;
- char *result,*prf1_file, *prf2_file;
- Alignment *F=NULL, *A1, *A2;
- FILE *fp;
- char command[10000];
- char *param;
-
- if ( M->executable2[0]=='\0')
- fprintf ( stderr, "\nERROR: profile_pair requires a method: thread_pair@EP@executable2@<method> [FATAL:%s]\n", PROGRAM);
-
-
- sprintf ( seq, "%s", in_seq);
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- A1=seq2R_template_profile(CL->S,s1);
- A2=seq2R_template_profile(CL->S,s2);
-
-
- prf1_file=vtmpnam (NULL);
- fp=vfopen (prf1_file, "w");
- if ( A1)
- {
- fprintf (fp, ">%s\n%s\n",(CL->S)->name[s1], aln2cons_seq_mat(A1, "blosum62mt"));
- for ( a=0; a< A1->nseq; a++)fprintf (fp, ">prf_seq1_%d\n%s\n", a, A1->seq_al[a]);
- }
- else
- {
- fprintf ( fp, ">%s\n%s\n", (CL->S)->name[s1], (CL->S)->seq[s1]);
- }
- vfclose (fp);
-
- prf2_file=vtmpnam (NULL);
- fp=vfopen (prf2_file, "w");
- if (A2)
- {
- fprintf (fp, ">%s\n%s\n",(CL->S)->name[s2], aln2cons_seq_mat(A2, "blosum62mt"));
- for ( a=0; a< A2->nseq; a++)fprintf (fp, ">prf_seq2_%d\n%s\n", a, A2->seq_al[a]);
- }
- else
- {
- fprintf ( fp, ">%s\n%s\n", (CL->S)->name[s2], (CL->S)->seq[s2]);
- }
- vfclose (fp);
-
- result=vtmpnam (NULL);
- if ( M->param)
- {
- param=vcalloc(strlen (M->param)+1, sizeof (char));
- sprintf ( param, "%s", M->param);
- param=substitute ( param, " ", "");
- param=substitute ( param, "\n", "");
- }
-
- sprintf ( command, "tc_generic_method.pl -mode=profile_pair -method=%s %s%s %s%s %s%s -param=%s -tmpdir=%s", M->executable2,M->in_flag,prf1_file, M->in_flag2,prf2_file,M->out_flag, result, param, get_tmp_4_tcoffee());
-
- my_system ( command);
-
-
-
- if ( !check_file_exists (result))
- {
- fprintf ( stderr, "\n\tprofile_pair/%s failed:\n\t%s\n",M->executable2, command);
- myexit (EXIT_FAILURE);
- }
- else if ( is_lib (result))
- {
- CL=read_constraint_list_file(CL,result);
- }
- else if ( is_aln (result))
- {
- F=main_read_aln (result, NULL);
- char *name1, *name2;
- name1=(CL->S)->name[s1];
- name2=(CL->S)->name[s2];
-
- fp=vfopen (result, "w");
- for ( a=0; a< F->nseq; a++)
- if (strm ( F->name[a], name1) || strm (F->name[a], name2))
- fprintf ( fp, ">%s\n%s\n", F->name[a], F->seq_al[a]);
- vfclose (fp);
- free_aln (F);
- F=main_read_aln (result, NULL);
- CL=aln2constraint_list (F, CL, "sim");
- free_aln (F);
- }
- return CL;
- }
-Constraint_list * pdbid_pair (TC_method *M , char *in_seq, Constraint_list *CL)
- {
-
- char seq[1000];
-
- int s1, s2;
- char *result, *pdb1, *pdb2;
- Alignment *F=NULL;
- char *command;
-
-
- if ( M->executable2[0]=='\0')
- {
- fprintf ( stderr, "\nERROR: pdbid_pair requires a structural alignment method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
- myexit (EXIT_FAILURE);
- }
- sprintf ( seq, "%s", in_seq);
-
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- pdb1=seq2P_pdb_id(CL->S,s1);
- pdb2=seq2P_pdb_id(CL->S,s2);
-
- if (!is_pdb_name (pdb1) || !is_pdb_name(pdb2))
- {
- return CL;
- }
-
-
- result=vtmpnam (NULL);
- command = vcalloc ( 1000, sizeof (char));
-
- sprintf ( command, "tc_generic_method.pl -mode=pdbid_pair -method=%s %s%s %s%s %s%s -email=%s -cache=%s -tmpdir=%s", M->executable2,M->in_flag,pdb1, M->in_flag2,pdb2,M->out_flag, result,getenv ("EMAIL"),get_cache_dir(), get_tmp_4_tcoffee());
- my_system ( command);
- vfree (command);
- if (file_is_empty (result))return CL;
- else
- {
- F=main_read_aln (result, NULL);
-
- if ( !F)
- {
- fprintf ( stderr, "\n\tpdb_pair/%s failed:\n\t%s\n",M->executable2, command);
- }
- else
- {
-
- sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
- sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
- CL=aln2constraint_list (F, CL, "sim");
- }
- free_aln (F);
- }
- return CL;
- }
-
-Constraint_list * pdb_pair (TC_method *M , char *in_seq, Constraint_list *CL)
- {
-
- char seq[1000];
- int s1, s2;
- char *result, *pdb1,*pdb1_file, *pdb2, *pdb2_file;
- Alignment *F=NULL;
-
-
- char command[10000];
-
- if ( M->executable2[0]=='\0')
- {
- fprintf ( stderr, "\nERROR: pdb_pair requires a structural alignment method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
- myexit (EXIT_FAILURE);
- }
-
- sprintf ( seq, "%s", in_seq);
-
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- pdb1=seq2P_template_file(CL->S,s1);
- pdb2=seq2P_template_file(CL->S,s2);
- if ( !pdb1 || !pdb2) return CL;
-
-
- pdb1_file=vtmpnam (NULL);
- pdb2_file=vtmpnam (NULL);
-
-
-
- sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -nodiagnostic > %s", pdb1, pdb1_file);
- my_system (command);
-
-
-
- sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -nodiagnostic > %s", pdb2, pdb2_file);
- my_system (command);
-
-
- result=vtmpnam (NULL);
-
- sprintf ( command, "tc_generic_method.pl -mode=pdb_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pdb1_file, M->in_flag2,pdb2_file,M->out_flag, result, get_tmp_4_tcoffee());
- my_system ( command);
-
- F=main_read_aln (result, NULL);
-
- if ( !F)
- {
- fprintf ( stderr, "\n\tpdb_pair/%s failed:\n\t%s\n",M->executable2, command);
- }
- else
- {
-
- sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
- sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
- CL=aln2constraint_list (F, CL, "sim");
- }
-
-
- free_aln (F);
- return CL;
- }
-
-Constraint_list * seq_msa (TC_method *M , char *in_seq, Constraint_list *CL)
-{
- char seq[1000];
- char db[1000];
- char *infile, *outfile;
- int a, n, s;
- Alignment *F=NULL;
- FILE *fp;
- char command[1000];
- static char *db_file;
-
-
- if (!db_file)
- {
- db_file=vtmpnam (NULL);
- fp=vfopen (db_file, "w");
- for (a=0; a<(CL->S)->nseq; a++)fprintf ( fp, ">%s\n%s\n", (CL->S)->name[a], (CL->S)->seq[a]);
- vfclose (fp);
- }
-
- infile=vtmpnam (NULL);
- outfile=vtmpnam (NULL);
-
- sprintf ( seq, "%s", in_seq);
-
- n=atoi(strtok (seq,SEPARATORS));
-
- fp=vfopen (infile, "w");
- for ( a=0; a<n; a++)
- {
- s=atoi(strtok (NULL,SEPARATORS));
- fprintf (fp, ">%s\n%s\n", (CL->S)->name[s], (CL->S)->seq[s]);
- }
- vfclose (fp);
-
- if (strstr (M->executable2, "blast"))sprintf ( db, " -database=%s", db_file);
- else db[0]='\0';
-
- sprintf ( command, "t_coffee -other_pg tc_generic_method.pl -mode=%s -method=%s %s %s%s %s %s%s -tmpdir=%s %s %s", M->executable, M->executable2, M->param1,M->in_flag,infile,M->param2,M->out_flag,outfile,get_tmp_4_tcoffee(),db,M->param);
-
-
-
-
- //HERE ("%s", command);exit (0);
- //sprintf ( command, "t_coffee -other_pg tc_generic_method.pl -mode=seq_msa -method=%s %s%s %s%s -tmpdir=%s %s", M->executable2, M->in_flag, infile, M->out_flag, outfile, get_tmp_4_tcoffee(), M->param);
- my_system (command);
-
-
- if ( strm (M->out_mode, "aln") || strm (M->out_mode, "A"))
- {
- F=main_read_aln (outfile, NULL);
- if ( !F)
- {
- fprintf ( stderr, "\n\tseq_msa/%s failed:\n\t%s\n", M->executable2,command);
- }
- else
- {
- CL=aln2constraint_list (F, CL, "sim");
- }
- free_aln (F);
- }
- else if ( strm (M->out_mode, "fL")|| strm (M->out_mode, "lib"))
- {
- Constraint_list *NCL;
- NCL=read_constraint_list_file(CL,outfile);
- if ( !NCL)
- {
- fprintf ( stderr, "\n\tseq_msa/%s failed:\n\t%s\n", M->executable2,command);
- }
- else
- {
- CL=NCL;
- }
- }
- return CL;
-}
-Constraint_list * thread_pair (TC_method *M , char *in_seq, Constraint_list *CL)
- {
-
- char seq[1000];
- int s1, s2;
-
- if ( M->executable2[0]=='\0')
- {
- fprintf ( stderr, "\nERROR: thread_pair requires a threading method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
- myexit (EXIT_FAILURE);
- }
-
- sprintf ( seq, "%s", in_seq);
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- CL=thread_pair2(M,s1, s2, CL);
- CL=thread_pair2(M,s2, s1, CL);
-
- return CL;
- }
-
-
-
-Constraint_list* thread_pair2 ( TC_method *M, int s1, int s2, Constraint_list *CL)
- {
- char *result, *pep, *pdb, *pdb1;
- Alignment *F=NULL;
- Sequence *STL;
- FILE *fp;
- char command[10000];
-
- STL=(CL)?CL->STRUC_LIST:NULL;
-
-
- if ( !(CL->S) || !((CL->S)->T[s1]) || !((CL->S)->T[s1])->P || !seq2P_template_file(CL->S,s1))return CL;
- else pdb1=seq2P_template_file(CL->S,s1);
-
-
- pdb=vtmpnam (NULL);
- result=vtmpnam (NULL);
- pep=vtmpnam (NULL);
-
- sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -nodiagnostic > %s", pdb1, pdb);
- my_system (command);
-
-
-
- fp=vfopen (pep, "w");
- fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s2],(CL->S)->seq[s2] );
- vfclose (fp);
- sprintf ( command, "tc_generic_method.pl -mode=thread_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pep, M->in_flag2,pdb,M->out_flag, result, get_tmp_4_tcoffee());
- my_system ( command);
- F=main_read_aln (result, NULL);
-
- if ( !F)
- {
- fprintf ( stderr, "\n\tthread_pair/%s failed:\n\t%s\n", M->executable2,command);
- }
- else
- {
- sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
- sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
- CL=aln2constraint_list (F, CL, "sim");
- }
-
-
- free_aln (F);
- return CL;
- }
-
-Constraint_list * lsqman_pair ( char *in_seq, Constraint_list *CL)
- {
- FILE *fp;
- static CLIST_TYPE *entry;
- char command[STRING];
- char seq[1000];
- int s1, s2;
- char *seq_file, *lsqman_result, *tmp_name1;
- Alignment *F=NULL;
- int n_failure=0;
-
- sprintf ( seq, "%s", in_seq);
-
-
- if ( !entry)entry=vcalloc ( LIST_N_FIELDS, sizeof ( CLIST_TYPE ));
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
-
- tmp_name1=vcalloc (100, sizeof (char));
- sprintf ( tmp_name1, "%s_%s.lsqman_aln", (CL->S)->name[s1], (CL->S)->name[s2]);
- if ( check_file_exists ( tmp_name1) && (F=main_read_aln(tmp_name1, NULL))!=NULL)
- {
- free_aln(F);
- lsqman_result=tmp_name1;
- }
-
- else
- {
- seq_file=vtmpnam (NULL);
- lsqman_result=tmp_name1;
- fp=vfopen (seq_file, "w");
- fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s1],(CL->S)->seq[s2] );
- vfclose (fp);
- sprintf ( command, "%s -pdb %s -pep %s > %s%s", LSQMAN_4_TCOFFEE, (CL->S)->name[s1], seq_file,get_cache_dir(), lsqman_result);
-
- while (!F)
- {
- my_system ( command);
- F=main_read_aln (lsqman_result, NULL);
- if ( !F)
- {
- fprintf ( stderr, "\n\tlsqman failed: will be retried");
- if ( n_failure==0)fprintf ( stderr, "\n\t%s", command);
- n_failure++;
- if ( n_failure==10)
- {
- fprintf ( stderr, "\nCould not run Fugue: will replace it with slow_pair\n");
- vremove (lsqman_result);
- return NULL;
- }
- }
- free_aln (F);
- }
- vremove( seq_file);
-
- }
-
-
- F=main_read_aln(lsqman_result, NULL);
- /*sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
- sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
- */
- CL=aln2constraint_list (F, CL, "100");
- free_aln (F);
- return CL;
- }
-
-Constraint_list * sap_pair (char *seq, char *weight, Constraint_list *CL)
- {
- register int a;
- FILE *fp;
- char full_name[FILENAMELEN];
- char *tmp_pdb1, *tmp_pdb2;
- char *sap_seq1, *sap_seq2;
- char *sap_lib, *tmp_name, *tmp_name1, *tmp_name2;
- char *buf=NULL;
- int s1, s2, r1=0, r2=0;
- int sim=0, tot=0, score=0;
-
- char program[STRING];
- char *string1, *string2, *string3, *string4, *string5;
- int max_struc_len=10000;
- char *template1, *template2;
- int c;
-
- /*check_program_is_installed ( "sap" ,SAP_4_TCOFFEE, "SAP_4_TCOFFEE",MAIL, IS_FATAL);*/
-
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- template1=seq2T_value(CL->S,s1, "template_name", "_P_");
- template2=seq2T_value(CL->S,s2, "template_name", "_P_");
-
-
- if (!template1 || !template2) return CL;
-
-
- declare_name (string1);
- declare_name (string2);
- declare_name (string3);
- declare_name (string4);
- declare_name (string5);
-
-#ifndef SAP_4_TCOFFEE
- if ( getenv ( "SAP_4_TCOFFEE")==NULL)crash ("SAP_4_TCOFFEE IS NOT DEFINED");
- else sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
-#else
- if ( getenv ( "SAP_4_TCOFFEE")==NULL)sprintf (program, "%s", SAP_4_TCOFFEE);
- else sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
-#endif
-
-
-
- tmp_name1=vcalloc (100, sizeof (char));
- sprintf ( tmp_name1, "%s_%s.sap_results",template1,template2);
- tmp_name2=vcalloc (100, sizeof (char));
- sprintf ( tmp_name2, "%s_%s.sap_results",template2,template1);
-
-
-
- if (is_sap_file (tmp_name1))
- {
- tmp_name=tmp_name1;
- }
- else if ( is_sap_file (tmp_name2))
- {
- tmp_name=tmp_name2;
- SWAP (s1, s2);
- }
- else
- {
- char local1[100];
- char local2[100];
- char cdir[1000];
-
- tmp_name=tmp_name1;
-
- tmp_pdb1=normalize_pdb_file(seq2P_template_file(CL->S,s1),(CL->S)->seq[s1], vtmpnam (NULL));
- tmp_pdb2=normalize_pdb_file(seq2P_template_file(CL->S,s2),(CL->S)->seq[s2], vtmpnam (NULL));
- sprintf ( full_name, "%s%s", get_cache_dir (), tmp_name);
-
- //pb: sap crashes when file names are too long
- //solution: create shorter names
- //To ease server, create files in tmp dir
-
- getcwd (cdir, sizeof(char)*1000);
- chdir (get_cache_dir());
-
- sprintf (local1, "%d_1.%d.sap_tmp", getpid(), rand()%10000);
- sprintf (local2, "%d_2.%d.sap_tmp", getpid(), rand()%10000);
- printf_system ("cp %s %s", tmp_pdb1, local1);
- printf_system ("cp %s %s", tmp_pdb2, local2);
- printf_system ("%s %s %s >%s 2>/dev/null::IGNORE_FAILURE::",program,local1,local2, full_name);
- printf_system ("rm %s", local1);
- printf_system ("rm %s", local2);
- chdir (cdir);
-
- if ( !check_file_exists (full_name) || !is_sap_file(full_name))
- {
- add_warning ( stderr, "SAP failed to align: %s against %s [%s:WARNING]\n", seq2P_template_file(CL->S,s1),seq2P_template_file(CL->S,s2), PROGRAM);
- if ( check_file_exists (full_name))add2file2remove_list (full_name);
- return CL;
- }
- if ( flag_file2remove_is_on())add2file2remove_list (full_name);
- remove ("super.pdb");
- }
-
-
-
- sap_seq1=vcalloc (max_struc_len, sizeof (char));
- sap_seq2=vcalloc (max_struc_len, sizeof (char));
-
- fp=find_token_in_file ( tmp_name, NULL, "Percent");
- fp=find_token_in_file ( tmp_name, fp , "Percent");
- while ( (c=fgetc (fp))!='\n' && c!=EOF)fprintf ( fp, "%c", c);
- while ((buf=vfgets (buf, fp)))
- {
-
- if ( !strstr (buf, "eighted") && !strstr (buf, "RMSd"))
- {
- remove_charset (buf, "!alnum");
- r1=buf[0];
- r2=buf[strlen(buf)-1];
- }
- else
- continue;
-
- sim+=(r1==r2)?1:0;
- if ( tot>max_struc_len)
- {max_struc_len+=max_struc_len;
- sap_seq1=vrealloc ( sap_seq1, sizeof(char)*max_struc_len);
- sap_seq2=vrealloc ( sap_seq2, sizeof(char)*max_struc_len);
- }
- sap_seq1[tot]=r1;
- sap_seq2[tot]=r2;
- tot++;
- }
- sim=(sim*100)/tot;
-
- if ( is_number (weight))score=atoi(weight);
- else if ( strstr ( weight, "OW"))
- {
- int ow;
- sscanf ( weight, "OW%d", &ow);
- score=sim*ow;
- }
- else
- {
- score=sim;
- }
-
-
- vfclose (fp);
- sap_seq1[tot]=sap_seq2[tot]='\0';
-
-
- fp=vfopen ( sap_lib=vtmpnam(NULL), "w");
- fprintf (fp, "! TC_LIB_FORMAT_01\n");
- fprintf (fp, "2\n");
- fprintf (fp, "%s %d %s\n", (CL->S)->name[s2],(int)strlen (sap_seq1), sap_seq1);
- fprintf (fp, "%s %d %s\n", (CL->S)->name[s1],(int)strlen (sap_seq2), sap_seq2);
- fprintf (fp, "#1 2\n");
-
-
-
- //HERE ("\n%s\n%s\n%s\n\n%s\n%s\n%s", (CL->S)->name[s1],(CL->S)->seq[s1], sap_seq1, (CL->S)->name[s2],(CL->S)->seq[s2], sap_seq2);exit (0);
- for ( a=0; a< tot; a++)
- {
- fprintf (fp, "%d %d %d 1 0\n", a+1, a+1, score);
- }
-
- fprintf (fp, "! CPU 0\n");
- fprintf (fp, "! SEQ_1_TO_N\n");
- vfclose (fp);
-
- CL=read_constraint_list_file(CL,sap_lib);
- vremove (sap_lib);
-
- vfree ( string1);vfree ( string2);vfree ( string3);vfree ( string4);vfree ( string5);
- vfree (sap_seq1); vfree(sap_seq2);vfree (tmp_name1); vfree(tmp_name2);
- vfree (buf);
- return CL;
- }
-
-
-
-Constraint_list *rna_pair (TC_method *M ,
- char *in_seq,
- Constraint_list *CL)
-{
- char seq[1000];
- int s1, s2;
- char *result, *pdb1, *pdb2, *pdb1_file, *pdb2_file;
- Alignment *F=NULL;
-
-
-
- char command[10000];
-
- if ( M->executable2[0]=='\0')
- {
- fprintf ( stderr, "\nERROR: rna_pair requires a structural alignment method: rna_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
- myexit (EXIT_FAILURE);
- }
-
- sprintf ( seq, "%s", in_seq);
-
-
- atoi(strtok (seq,SEPARATORS));
- s1=atoi(strtok (NULL,SEPARATORS));
- s2=atoi(strtok (NULL,SEPARATORS));
-
- pdb1=seq2P_template_file(CL->S,s1);
- pdb1_file=vtmpnam (NULL);
- sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -nodiagnostic > %s", pdb1, pdb1_file);
- my_system (command);
- //
- pdb2=seq2P_template_file(CL->S,s2);
- pdb2_file=vtmpnam (NULL);
- sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST -nodiagnostic > %s", pdb2, pdb2_file);
- my_system (command);
-
- result=vtmpnam (NULL);
-
- sprintf ( command, "tc_generic_method.pl -mode=rna_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pdb1_file, M->in_flag2,pdb2_file,M->out_flag, result, get_tmp_4_tcoffee());
-
- my_system ( command);
-
- F=main_read_aln (result, NULL);
-
-
- if ( !F)
- {
- fprintf ( stderr, "\n\trna_pair/%s failed:\n\t%s\n",M->executable2, command);
- }
- else
- {
- sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
- sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
- CL=aln2constraint_list (F, CL, "sim");
- }
-
-
- free_aln (F);
- return CL;
-}
-
-
-/******************************************************************/
-/* GENERIC PAIRWISE METHODS */
-/* */
-/* */
-/******************************************************************/
-
-
-
-Constraint_list * best_pair4rna(Job_TC *job)
-{
- int n,a;
-
-
- static char *seq;
- Alignment *A;
- Constraint_list *PW_CL;
- Constraint_list *CL, *RCL;
- char *seq_in;
- Sequence *S;
- TC_method *M, *sara_pairM, *proba_pairM;
- int*seqlist;
-
- int s1, s2;
- Template *T1, *T2;
- int ml=0;
- struct X_template *r1, *r2, *p1, *p2;
- static int **blosum;
-
- if (!seq)seq=vcalloc (100, sizeof (char));
-
- A=(job->io)->A;
- M=(job->param)->TCM;
- CL=(job->io)->CL;
- S=CL->S;
- for (a=0; a<S->nseq; a++)ml=MAX(ml, strlen (S->name[a]));
-
-
- if ( !strm ( retrieve_seq_type(), "RNA") && !strm ( retrieve_seq_type(), "DNA"))
- printf_exit (EXIT_FAILURE, stderr, "ERROR: RNA Sequences Only with best4rna_pair [FATAL:%s]\n",PROGRAM);
-
-
- seq_in=(job->param)->seq_c;
- sprintf (seq, "%s", seq_in);
- seqlist=string2num_list (seq);
- n=seqlist[1];
- if ( n!=2){fprintf ( stderr, "\nERROR: best_pair can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
- s1=seqlist[2];
- s2=seqlist[3];
-
- T1=S->T[s1];
- T2=S->T[s2];
- r1=T1->R;
- r2=T2->R;
- p1=T1->P;
- p2=T2->P;
-
- PW_CL=((job->param)->TCM)->PW_CL;
- CL=(job->io)->CL;
-
- if (!blosum)blosum=read_matrice ("blosum62mt");
-
-// id=idscore_pairseq (S->seq[s1], S->seq[s2],-10,-1,blosum, "sim");
-
-
- proba_pairM=method_file2TC_method(method_name2method_file ("proba_pair"));
- proba_pairM->PW_CL=method2pw_cl(proba_pairM, CL);
-
- sara_pairM=method_file2TC_method(method_name2method_file ("sara_pair"));
- sara_pairM->PW_CL=method2pw_cl(sara_pairM, CL);
-
- if ( p1 && p2)
- {
- //Avoid Structural Tem
- T1->R=NULL;
- T2->R=NULL;
- fprintf ( stderr, "\n\t%-*s %-*s: Structure Based Alignment\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=sara_pairM;
- }
- else
- {
- fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=proba_pairM;
- }
-
- RCL=seq2list (job);
- T1->R=r1;
- T2->R=r2;
-
- return RCL;
-}
-
-Constraint_list * best_pair4prot (Job_TC *job)
-{
- int n,a;
-
-
- static char *seq;
- Alignment *A;
- Constraint_list *PW_CL;
- Constraint_list *CL, *RCL;
- char *seq_in;
- Sequence *S;
- TC_method *M, *sap_pairM, *proba_pairM;
- int*seqlist;
-
- int id, s1, s2;
- Template *T1, *T2;
- int ml=0;
- struct X_template *r1, *r2, *p1, *p2;
- static int **blosum;
-
- if (!seq)seq=vcalloc (100, sizeof (char));
-
- A=(job->io)->A;
- M=(job->param)->TCM;
- CL=(job->io)->CL;
- S=CL->S;
- for (a=0; a<S->nseq; a++)ml=MAX(ml, strlen (S->name[a]));
-
-
- if ( strm ( retrieve_seq_type(), "DNA") ||strm ( retrieve_seq_type(), "RNA") )printf_exit (EXIT_FAILURE, stderr, "ERROR: Protein Sequences Only with bestprot_pair [FATAL:%s]\n",PROGRAM);
-
-
- seq_in=(job->param)->seq_c;
- sprintf (seq, "%s", seq_in);
- seqlist=string2num_list (seq);
- n=seqlist[1];
- if ( n!=2){fprintf ( stderr, "\nERROR: best_pair can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
- s1=seqlist[2];
- s2=seqlist[3];
-
- T1=S->T[s1];
- T2=S->T[s2];
- r1=T1->R;
- r2=T2->R;
- p1=T1->P;
- p2=T2->P;
-
- PW_CL=((job->param)->TCM)->PW_CL;
- CL=(job->io)->CL;
-
- if (!blosum)blosum=read_matrice ("blosum62mt");
-
- id=idscore_pairseq (S->seq[s1], S->seq[s2],-10,-1,blosum, "sim");
-
-
- proba_pairM=method_file2TC_method(method_name2method_file ("proba_pair"));
- proba_pairM->PW_CL=method2pw_cl(proba_pairM, CL);
-
- sap_pairM=method_file2TC_method(method_name2method_file ("sap_pair"));
- sap_pairM->PW_CL=method2pw_cl(sap_pairM, CL);
-
- if ( id>80)
- {
- //Hide The Template
- T1->R=NULL;
- T2->R=NULL;
- fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=proba_pairM;
- }
- else if ( p1 && p2)
- {
- //Avoid Structural Tem
- T1->R=NULL;
- T2->R=NULL;
- fprintf ( stderr, "\n\t%-*s %-*s: Structure Based Alignment\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=sap_pairM;
- }
- else if ( r1 || r2)
- {
- fprintf ( stderr, "\n\tt%-*s %-*s: PSIBLAST Profile Alignment\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=proba_pairM;
- }
- else
- {
- fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment (No Profile)\n", ml,S->name[s1], ml,S->name[s2]);
- (job->param)->TCM=proba_pairM;
- }
-
- RCL=seq2list (job);
- T1->R=r1;
- T2->R=r2;
-
- return RCL;
-}
-
-
-Alignment * fast_pair (Job_TC *job)
- {
- int s, n,a;
- int score;
- static int **l_s;
- static int *ns;
- char seq[1000];
- Alignment *A;
- Constraint_list *PW_CL;
- Constraint_list *CL;
- char *seq_in;
- Sequence *S;
- TC_method *M;
- int*seqlist;
- char **buf;
-
- A=(job->io)->A;
-
- M=(job->param)->TCM;
- PW_CL=((job->param)->TCM)->PW_CL;
- CL=(job->io)->CL;
- seq_in=(job->param)->seq_c;
-
-
- sprintf (seq, "%s", seq_in);
- seqlist=string2num_list (seq);
- n=seqlist[1];
- if ( n!=2){fprintf ( stderr, "\nERROR: fast_pw_aln can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
-
- S=(CL)->S;
-
- if (!A) {A=declare_aln (CL->S);}
- if ( !ns)
- {
- ns=vcalloc ( 2, sizeof (int));
- l_s=declare_int (2,(CL->S)->nseq);
- }
- buf=vcalloc ( S->nseq, sizeof (char*));
-
- for ( a=0; a< n; a++)
- {
- s=seqlist[a+2];
- if ( strm (M->seq_type, "G"))
- {
- buf[s]=S->seq[s];
- S->seq[s]=((((S->T[s])->G)->VG)->S)->seq[0];
- }
- else
- buf[s]=S->seq[s];
-
-
- sprintf ( A->seq_al[a], "%s",S->seq[s]);
- sprintf ( A->name[a], "%s", (CL->S)->name[s]);
- A->order[a][0]=s;
- }
-
- A->S=CL->S;
- PW_CL->S=CL->S;
- A->CL=CL;
- A->nseq=n;
- ns[0]=ns[1]=1;
- l_s[0][0]=0;
- l_s[1][0]=1;
-
- //Preprocessing of the sequences
- if (PW_CL->reverse_seq)
- {
- invert_string2(A->seq_al[0]);
- invert_string2(A->seq_al[1]);
- invert_string2 ((CL->S)->seq[A->order[0][0]]);
- invert_string2 ((CL->S)->seq[A->order[1][0]]);
- }
- if (PW_CL->extend_seq)//use te alphabet extension for nucleic acids
- {
-
- extend_seqaln (A->S,NULL);
- extend_seqaln (NULL,A);
- }
-
- score=pair_wise ( A, ns, l_s, PW_CL);
- //PostProcessing of the sequences
- if (PW_CL->reverse_seq)
- {
-
- invert_string2(A->seq_al[0]);
- invert_string2(A->seq_al[1]);
- invert_string2 ((CL->S)->seq[A->order[0][0]]);
- invert_string2 ((CL->S)->seq[A->order[1][0]]);
- }
- if (PW_CL->extend_seq)
- {
- unextend_seqaln (A->S,NULL);
- unextend_seqaln (NULL,A);
- }
- A->nseq=n;
-
- for ( a=0; a<S->nseq; a++)
- {
- if ( !buf[a] || buf[a]==S->seq[a]);
- else S->seq[a]=buf[a];
- }
- vfree (buf);vfree (seqlist);
- return A;
-
- }
-Alignment * align_two_aln ( Alignment *A1, Alignment *A2, char *in_matrix, int gop, int gep, char *in_align_mode)
- {
- Alignment *A=NULL;
- Constraint_list *CL;
- Sequence *S;
- int a;
- int *ns;
- int **ls;
- static char *matrix;
- static char *align_mode;
-
- if (!matrix)matrix=vcalloc ( 100, sizeof (char));
- if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
-
-
-
- sprintf ( matrix, "%s", in_matrix);
- sprintf ( align_mode, "%s", in_align_mode);
-
- CL=vcalloc ( 1, sizeof (Constraint_list));
- CL->pw_parameters_set=1;
- CL->M=read_matrice (matrix);
- CL->matrices_list=declare_char (10, 10);
-
- CL->evaluate_residue_pair=evaluate_matrix_score;
- CL->get_dp_cost=consensus_get_dp_cost;
- CL->normalise=1;
-
- CL->extend_jit=0;
- CL->maximise=1;
- CL->gop=gop;
- CL->gep=gep;
- CL->TG_MODE=2;
- sprintf (CL->matrix_for_aa_group, "vasiliky");
- CL->use_fragments=0;
- CL->ktup=5;
- if ( !CL->use_fragments)CL->diagonal_threshold=0;
- else CL->diagonal_threshold=6;
-
- sprintf (CL->dp_mode, "%s", align_mode);
-
- A=copy_aln (A1, A);
- A=stack_aln (A, A2);
- CL->S=fill_sequence_struc(A->nseq, A->seq_al,A->name);
-
- ns=vcalloc ( 2, sizeof(int));
- ls=declare_int ( 2,A->nseq);
- ns[0]=A1->nseq;
- ns[1]=A2->nseq;
- for ( a=0; a<ns[0]; a++)
- ls[0][a]=a;
- for ( a=0; a< ns[1]; a++)
- ls[1][a]=a+A1->nseq;
-
- A->score_aln=pair_wise (A, ns, ls,CL);
-
- vfree (ns);
- free_int (ls, -1);
- S=free_constraint_list (CL);
- free_sequence (S,-1);
- A->S=NULL;
- return A;
- }
-
-static int align_two_seq_keep_case;
-void toggle_case_in_align_two_sequences(int value)
-{
- align_two_seq_keep_case=value;
-}
-Alignment * align_two_sequences ( char *seq1, char *seq2, char *in_matrix, int gop, int gep, char *in_align_mode)
- {
- static Alignment *A;
- Constraint_list *CL;
- Sequence *S;
-
- int *ns;
- int **l_s;
-
- char **seq_array;
- char **name_array;
- static char *matrix;
- static int **M;
-
- static char *align_mode;
-
- if (!matrix)matrix=vcalloc ( 100, sizeof (char));
- if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
- sprintf ( align_mode, "%s", in_align_mode);
-
- CL=vcalloc ( 1, sizeof (Constraint_list));
-
- CL->pw_parameters_set=1;
-
- CL->matrices_list=declare_char (10, 10);
-
-
- if ( !strm (matrix, in_matrix))
- {
- sprintf ( matrix,"%s", in_matrix);
- M=CL->M=read_matrice (matrix);
-
- }
- else
- {
- CL->M=M;
- }
-
- if (strstr (in_align_mode, "cdna"))
- CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
- else
- CL->evaluate_residue_pair=evaluate_matrix_score;
-
- CL->get_dp_cost=get_dp_cost;
- CL->extend_jit=0;
- CL->maximise=1;
- CL->gop=gop;
- CL->gep=gep;
- CL->TG_MODE=2;
- sprintf (CL->matrix_for_aa_group, "vasiliky");
- CL->use_fragments=0;
- CL->ktup=3;
- if ( !CL->use_fragments)CL->diagonal_threshold=0;
- else CL->diagonal_threshold=6;
-
- sprintf (CL->dp_mode, "%s", align_mode);
-
- seq_array=declare_char ( 2, MAX(strlen(seq1), strlen (seq2))+1);
- sprintf (seq_array[0], "%s",seq1);
- sprintf (seq_array[1],"%s", seq2);
- ungap_array(seq_array,2);
- if (align_two_seq_keep_case !=KEEP_CASE)string_array_lower(seq_array,2);
-
- name_array=declare_char (2, STRING);
- sprintf ( name_array[0], "A");
- sprintf ( name_array[1], "B");
-
-
- ns=vcalloc ( 2, sizeof(int));
- l_s=declare_int ( 2, 1);
- ns[0]=ns[1]=1;
- l_s[0][0]=0;
- l_s[1][0]=1;
-
-
-
- CL->S=fill_sequence_struc(2, seq_array, name_array);
-
- A=seq2aln(CL->S, NULL, 1);
-
- ungap (A->seq_al[0]);
- ungap (A->seq_al[1]);
-
-
-
- A->score_aln=pair_wise (A, ns, l_s,CL);
-
- vfree (ns);
- free_int (l_s, -1);
- free_char (name_array, -1);free_char ( seq_array,-1);
-
- CL->M=NULL;
- S=free_constraint_list (CL);
- free_sequence (S,-1);
- A->S=NULL;
- return A;
- }
-
-
-NT_node make_root_tree ( Alignment *A,Constraint_list *CL,int gop, int gep,Sequence *S, char *tree_file,int maximise)
-{
- NT_node **T=NULL;
- T=make_tree (A, CL, gop, gep,S,tree_file,maximise);
- (T[3][0])->nseq=S->nseq;
- return T[3][0];
-}
-NT_node ** make_tree ( Alignment *A,Constraint_list *CL,int gop, int gep,Sequence *S, char *tree_file,int maximise)
- {
- int a, b, ra, rb;
- NT_node **T=NULL;
- int **distances;
- int out_nseq;
- char **out_seq_name;
- char **out_seq;
-
-
- if ( !CL || !CL->tree_mode || !CL->tree_mode[0])
- {
- fprintf ( stderr, "\nERROR: No CL->tree_mode specified (make_tree::util_dp_drivers.c [FATAL:%s]", PROGRAM);
- myexit (EXIT_FAILURE);
- }
- else
- fprintf (CL->local_stderr , "\nMAKE GUIDE TREE \n\t[MODE=%s][",CL->tree_mode);
-
-
- if ( A->nseq==2)
- {
- int tot_node;
- char *tmp;
- FILE *fp;
- fprintf (CL->local_stderr, "---Two Sequences Only: Make Dummy Pair-Tree ---]");
- tmp=vtmpnam (NULL);
- fp=vfopen (tmp,"w");
- fprintf ( fp, "(%s:0.1, %s:0.1):0.1;\n",S->name[0], S->name[1]);
- vfclose (fp);
- T=read_tree (tmp, &tot_node, (CL->S)->nseq,(CL->S)->name);
-
- return T;
- }
- else if ( strm (CL->tree_mode, "cwph"))
- {
- return seq2cw_tree ( S, tree_file);
- }
-
- else if (strm ( CL->tree_mode, "upgma") || strm ( CL->tree_mode, "nj"))
- {
- out_nseq=S->nseq;
- out_seq_name=S->name;
- out_seq=S->seq;
-
- CL->DM=cl2distance_matrix (CL, NOALN,NULL,NULL,0);
-
- if ( CL->S!=S)
- {
- /*Shrink the distance matrix so that it only contains the required sequences*/
- distances=declare_int (S->nseq, S->nseq);
- for (a=0; a< S->nseq; a++)
- {
- ra=name_is_in_list ((S)->name[a],(CL->S)->name, (CL->S)->nseq, 100);
- for ( b=0; b< S->nseq; b++)
- {
- rb=name_is_in_list ((S)->name[b],(CL->S)->name, (CL->S)->nseq, 100);
- distances[a][b]=(CL->DM)->score_similarity_matrix[ra][rb];
- }
- }
- }
- else
- {
- distances=duplicate_int ( (CL->DM)->score_similarity_matrix, -1, -1);
- }
-
-
-
- distances=sim_array2dist_array (distances, MAXID*SCORE_K);
- distances=normalize_array (distances, MAXID*SCORE_K, 100);
- if ( strm (CL->tree_mode, "order"))
- {
- for ( a=0; a< S->nseq; a++)
- for ( b=0; b< S->nseq; b++)
- distances[b][a]=100;
- T=make_nj_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
- }
- else if ( strm (CL->tree_mode, "nj"))
- {
- T=make_nj_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
- }
- else if ( strm (CL->tree_mode, "upgma"))
- T=make_upgma_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
- else
- {
- printf_exit (EXIT_FAILURE, stderr, "ERROR: %s is an unknown tree computation mode [FATAL:%s]", CL->tree_mode, PROGRAM);
- }
- free_int (distances, out_nseq);
-
- }
-
- fprintf (CL->local_stderr , "DONE]\n");
- return T;
- }
-
-
-Alignment *recompute_local_aln (Alignment *A, Sequence *S,Constraint_list *CL, int scale, int gep)
- {
- int **coor;
- int a;
- Alignment *B;
-
- sort_constraint_list (CL, 0, CL->ne);
- coor=declare_int (A->nseq, 3);
- for ( a=0; a< A->nseq; a++)
- {
- coor[a][0]=A->order[a][0];
- coor[a][1]=A->order[a][1]+1;
- coor[a][2]=strlen(S->seq[A->order[a][0]])-coor[a][1];
- }
- B=stack_progressive_nol_aln_with_seq_coor(CL,0,0,S,coor,A->nseq);
- A=copy_aln ( B, A);
-
- free_Alignment(B);
- return A;
- }
-
-
-Alignment *stack_progressive_nol_aln_with_seq_coor(Constraint_list *CL,int gop, int gep,Sequence *S, int **seq_coor, int nseq)
- {
-
- static int ** local_coor1;
- static int ** local_coor2;
- if ( local_coor1!=NULL)free_int (local_coor1, -1);
- if ( local_coor2!=NULL)free_int (local_coor2, -1);
-
- local_coor1=get_nol_seq ( CL,seq_coor, nseq, S);
- local_coor2=minimise_repeat_coor ( local_coor1, nseq, S);
-
- return stack_progressive_aln_with_seq_coor(CL,gop, gep,S, local_coor2,nseq);
- }
-
-
-Alignment *stack_progressive_aln_with_seq_coor (Constraint_list*CL,int gop, int gep, Sequence *S, int **coor, int nseq)
- {
- Alignment *A=NULL;
-
- A=seq_coor2aln (S,NULL, coor, nseq);
-
- return stack_progressive_aln ( A,CL, gop, gep);
- }
-
-Alignment *est_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
- {
- int a,n;
- int**group_list;
- int *n_groups;
- char *seq;
- n_groups=vcalloc ( 2, sizeof (int));
- group_list=declare_int ( 2, A->nseq);
-
- n=A->nseq;
-
- n_groups[0]=1;
- n_groups[1]=1;
- group_list[0][0]=0;
- group_list[0][1]=1;
-
- group_list[1][0]=1;
- fprintf ( stderr, "\n");
- for ( a=1; a<n; a++)
- {
- sprintf ( A->seq_al[1], "%s", A->seq_al[a]);
- fprintf ( stderr, "\t[%30s]->[len=%5d]", A->name[a],(int)strlen ( A->seq_al[0]));
- pair_wise ( A,n_groups, group_list, CL);
-
- seq=dna_aln2cons_seq(A);
-
- sprintf ( A->seq_al[0], "%s", seq);
- vfree (seq);
- fprintf ( stderr, "\n");
- }
-
- A->nseq=1;
- return A;
- }
-
-void analyse_seq ( Alignment *A, int s)
- {
- int a, b, c;
- int r;
-
- int len=0;
-
-
- int state=0;
- int pstate=-1;
- float score=0;
-
- for ( a=0; a< A->len_aln; a++)
- {
- for ( b=0, c=0; b< s; b++)
- if ( !is_gap(A->seq_al[b][a])){c=1; break;}
-
- r=!is_gap(A->seq_al[s][a]);
-
- if ( r && c) state=1;
- else if ( !r && !c) state=2;
- else if ( !r && c) state=3;
- else if ( r && !c) state=4;
-
- if ( state !=pstate)
- {
- score+=len*len;
- len=0;
- }
- len+=r;
- pstate=state;
- }
- score=score/(float)(((A->S)->len[s]*(A->S)->len[s]));
- fprintf ( stderr, "[%.2f]", score);
-
- return;
- }
-
-Alignment *realign_aln ( Alignment*A, Constraint_list *CL)
-{
- int a, b, c;
- int *ns, **ls;
- A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
-
- ns=vcalloc (2, sizeof(int));
- ls=declare_int ( 2, A->nseq);
-
- for (a=0; a< A->nseq; a++)
- {
- ns[0]=A->nseq-1;
- for (c=0,b=0; b<A->nseq; b++)if (b!=a)ls[0][c++]=b;
- ungap_sub_aln ( A, ns[0], ls[0]);
-
- ns[1]=1;
- ls[1][0]=a;
- ungap_sub_aln ( A, ns[1], ls[1]);
- A->score_aln=pair_wise (A, ns, ls,CL);
- }
-
- vfree (ns); free_int (ls, -1);
- return A;
-}
-
-Alignment *realign_aln_random_bipart ( Alignment*A, Constraint_list *CL)
-{
- int *ns;
- int **l_s;
-
- int a,g;
-
- ns=vcalloc (2, sizeof (int));
- l_s=declare_int (2,A->nseq);
-
- for ( a=0; a< A->nseq; a++)
- {
- g=rand()%2;
- l_s[g][ns[g]++]=a;
- }
-
- fprintf ( stderr, "\n");
- ungap_sub_aln ( A, ns[0], l_s[0]);
- ungap_sub_aln ( A, ns[1], l_s[1]);
-
- /* //Display Groups
- for (a=0;a<2; a++)
- for (b=0; b<ns[a]; b++)
- fprintf ( stderr, "[%d-%d]", a, l_s[a][b]);
- */
- A->score_aln=pair_wise (A, ns, l_s,CL);
-
- vfree(ns);free_int(l_s, -1);
- return A;
-}
-Alignment *realign_aln_random_bipart_n ( Alignment*A, Constraint_list *CL, int n)
-{
- int *ns;
- int **ls;
- int *used;
-
- int a,b,c, p;
-
- if (n>=A->nseq)n=A->nseq/2;
- used=vcalloc (A->nseq, sizeof (int));
- c=0;
- while (c<n)
- {
- p=rand()%A->nseq;
- if (!used[p]){used[p]=1;c++;}
- }
- ns=vcalloc (2, sizeof (int));
- ls=declare_int (2,A->nseq);
- for (a=0; a<2; a++)
- {
- for (b=0; b<A->nseq; b++)
- if (used[b]==a)ls[a][ns[a]++]=b;
- }
- ungap_sub_aln ( A, ns[0], ls[0]);
- ungap_sub_aln ( A, ns[1], ls[1]);
-
-
- A->score_aln=pair_wise (A, ns, ls,CL);
- vfree(ns);free_int(ls, -1);vfree (used);
- return A;
-}
-int ** seq2ecl_mat (Constraint_list *CL);
-int ** seq2ecl_mat (Constraint_list *CL)
-
-{
- int a, b, n;
-
- Alignment *A;
- int *ns, **ls;
- int **dm;
-
- ns=vcalloc (2, sizeof (int));
- ls=declare_int ((CL->S)->nseq, 2);
-
- A=seq2aln (CL->S,NULL, RM_GAP);
- n=(CL->S)->nseq;
- dm=declare_int (n, n);
- for (a=0; a<(CL->S)->nseq-1; a++)
- for (b=a+1; b<(CL->S)->nseq; b++)
- {
- ns[0]=ns[1]=1;
- ls[0][0]=a;
- ls[1][0]=b;
- ungap (A->seq_al[a]);
- ungap (A->seq_al[b]);
- dm[a][b]=dm[b][a]=linked_pair_wise (A, ns, ls, CL);
- }
-
- return dm;
-}
-Alignment *realign_aln_clust ( Alignment*A, Constraint_list *CL)
-{
- int *ns;
- int **ls;
-
- int a, b, c,n;
- static int **rm, **dm, **target;
- int score;
-
-
-
- if (!A)
- {
- free_int (dm, -1); free_int (rm, -1);free_int (target, -1);
- dm=rm=target=NULL;
- }
-
-
- if (!rm)rm=seq2ecl_mat(CL);
- if (!dm)dm=declare_int (A->nseq, A->nseq);
- if (!target)target=declare_int (A->nseq*A->nseq, 3);
-
- ns=vcalloc (2, sizeof (int));
- ls=declare_int (2,A->nseq);
-
-
- for (a=0; a<A->nseq-1; a++)
- for (b=a+1; b<A->nseq; b++)
- {
- ns[0]=2;
- ls[0][0]=a;
- ls[0][1]=b;
- score=sub_aln2ecl_raw_score (A, CL, ns[0], ls[0]);
- dm[a][b]=dm[b][a]=MAX(0,(rm[a][b]-score));
- }
- for (n=0,a=0; a<A->nseq; a++)
- {
- for (b=a; b<A->nseq; b++, n++)
- {
-
- target[n][0]=a;
- target[n][1]=b;
- for ( c=0; c<A->nseq; c++)
- {
- if (c!=a && c!=b)target[n][2]+=dm[a][c]+dm[b][c];
- }
- }
- }
- sort_int_inv (target,3, 2, 0, n-1);
-
- for (a=0; a<A->nseq; a++)
- {
- if (target[a][0]==target[a][1])
- {
- ns[0]=1;
- ls[0][0]=target[a][0];
- }
- else
- {
- ns[0]=2;
- ls[0][0]=target[a][0]; ls[0][1]=target[a][1];
- }
-
- for (ns[1]=0,b=0; b<A->nseq; b++)
- {
- if (b!=target[a][0] && b!=target[a][1])ls[1][ns[1]++]=b;
- }
-
- ungap_sub_aln (A, ns[0], ls[0]);
- ungap_sub_aln (A, ns[1], ls[1]);
-
- A->score_aln=pair_wise (A, ns, ls,CL);
- fprintf ( stderr, "\nSEQ: %d %d SCORE=%d\n",target[a][0],target[a][1], aln2ecl_raw_score(A, CL));
- }
- return A;
-}
-
-int get_best_group ( int **used, Constraint_list *CL);
-int seq_aln_thr1(Alignment *A, int **used, int threshold, Constraint_list *CL);
-int seq_aln_thr2( Alignment*A, int **used, int threshold, int g, Constraint_list *CL);
-
-int get_best_group ( int **used, Constraint_list *CL)
-{
- int a,b,c,d,n, tot,stot, best_tot, best_seq, nseq;
- int ns[2];
- int *ls[2];
-
- best_seq=0;
- nseq=((CL->S)->nseq);
- tot=best_tot=0;
- for (a=0; a<nseq; a++)
- {
- if ( used[a][0]==-1)continue;
- for ( tot=0,b=0; b< nseq; b++)
- {
- if ( a==b) continue;
- if ( used[b][0]==-1)continue;
- ns[0]=used[a][1];
- ls[0]=used[a]+2;
- ns[1]=used[b][1];
- ls[1]=used[b]+2;
- for (stot=0, n=0,c=0; c<ns[0]; c++)
- for (d=0; d<ns[1]; d++, n++)
- {
- stot+=(CL->DM)->similarity_matrix[ls[0][c]][ls[1][d]];
- }
- if (n>0)stot/=n;
- tot+=stot;
- }
- if (tot>best_tot)
- {
- best_tot=tot;
- best_seq=a;
- }
- }
- return best_seq;
-}
-
-
-
-Alignment * seq2aln_group (Alignment *A, int N, Constraint_list *CL)
-{
-
- NT_node P;
- int a;
-
- char *list, **list2;
- Alignment *F;
-
-
- fprintf (CL->local_stderr, "\n##### DPA ##### Compute Fast Alignment");
- A=iterative_tree_aln (A,1, CL);
- fprintf (CL->local_stderr, "\n##### DPA ##### Identify Nodes");
- P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- set_node_score (A, P, "idmat_sim");
- fprintf (CL->local_stderr, "\n##### DPA ##### Split Nodes");
- list=split_nodes_nseq (A,P,N, list=vcalloc (P->nseq*200, sizeof (char)));
-
- list2=string2list (list);
- fprintf (CL->local_stderr, "\n##### DPA ##### Save Nodes");
-
- F=A;
- for (a=1; a<atoi (list2[0]); a++)
- {
- A->A=main_read_aln(list2[a], NULL);
- A=A->A;
- }
- fprintf (CL->local_stderr, "\n##### DPA ##### Finished");
- vfree (list); free_char (list2, -1);
-
- A=F;
- while (A)
- {
- A=A->A;
- }
-
- return F;
-}
-
-
-
-
-Alignment * seq_aln ( Alignment*A, int n,Constraint_list *CL)
-{
-
- int **used, a, t,n1, nseq;
-
-
- n1=nseq=(CL->S)->nseq;
- used=declare_int (nseq, nseq+3);
-
-
- for (a=0; a< nseq; a++)
- {
- used[a][1]=1;
- used[a][2]=a;
- }
-
-
- for (t=50; t>=0 && nseq>1; t-=5)
- {
- nseq=seq_aln_thr1 (A, used,t, CL);
- }
-
- vfree (used);
- return A;
-}
-
-int seq_aln_thrX(Alignment *A, int **used, int threshold, Constraint_list *CL)
-{
- int n=0,a;
- seq_aln_thr1(A,used,threshold,CL);
- for ( a=0; a< (CL->S)->nseq; a++)
- n+=(used[a][1]>0)?1:0;
-
- return n;
-}
-int seq_aln_thr1(Alignment *A, int **used, int threshold, Constraint_list *CL)
-{
- int a,g, nseq, n_groups;
- nseq=(CL->S)->nseq;
-
- g=get_best_group(used, CL);
-
- used[g][0]=1;
-
-
-
- while ( seq_aln_thr2 (A, used, threshold,g, CL)!=0)
- {
- g=get_best_group (used, CL);
- used[g][0]=1;
- }
-
- for (n_groups=0,a=0; a< nseq; a++)
- if ( used[a][1]!=0)
- {
- n_groups++;
- used[a][0]=0;
- }
- return n_groups;
-}
-
-
-int seq_aln_thr2( Alignment*A, int **used, int threshold, int g, Constraint_list *CL)
-{
- int a, b,c,d;
- int ns[2], *ls[2];
- int nseq, n_members;
- double sim;
-
- n_members=0;
-
- nseq=((CL->S)->nseq);
- used[g][0]=1;
- ns[0]=used[g][1];
- ls[0]=used[g]+2;
-
- for ( a=0; a< nseq; a++)
- {
- if (used[a][0]!=0);
- else
- {
- ns[1]=used[a][1];
- ls[1]=used[a]+2;
-
- ungap_sub_aln (A, ns[0], ls[0]);
- ungap_sub_aln (A, ns[1], ls[1]);
-
- A->score_aln=pair_wise (A, ns, ls,CL);
-
- for (sim=0,b=0; b<ns[0]; b++)
- {
- for (c=0; c<ns[1]; c++)
- {
- sim+=generic_get_seq_sim (A->seq_al[ls[0][b]], A->seq_al[ls[1][c]], NULL,"idmat_sim2");
- }
- }
- sim/=(double)(ns[0]*ns[1]);
- if (sim>=threshold)
- {
-
- used[g][1]+=ns[1];
- for (d=0; d<ns[1]; d++)
- ls[0][ns[0]++]=ls[1][d];
- used[a][0]=-1;
- used[a][1]=0;
- b=ns[0];c=ns[1];
- n_members++;
- }
-
- }
- }
-
- if (n_members>0)used[g][0]=-1;
- return n_members;
-}
-/****************************************************************************/
-/* */
-/* */
-/* Alignment Methods */
-/* */
-/* */
-/****************************************************************************/
-
-Alignment * tsp_aln (Alignment *A, Constraint_list *CL, Sequence *S)
-{
- int a, b ;
- int ** distances;
- int *ns, **ls;
- int **used;
-
- A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
- ns=vcalloc (2, sizeof (int));
- ls=declare_int (2, (CL->S)->nseq);
- used=declare_int ( A->nseq, 2);
-
-
- CL->DM=cl2distance_matrix (CL, NOALN,NULL,NULL,0);
- distances=declare_int (A->nseq+1, A->nseq+1);
- distances=duplicate_int ( (CL->DM)->score_similarity_matrix, -1, -1);
-
- for (a=0; a< A->nseq; a++)
- {
- used[a][0]=a;
- for (b=0; b< A->nseq; b++)
- {
- used[a][1]+=distances[a][b];
- }
- }
-
- sort_int_inv (used,2,1,0,(CL->S)->nseq-1);
-
- ls[0][ns[0]++]=used[0][0];
- ns[1]=1;
-
- for (a=1; a< S->nseq; a++)
- {
- fprintf ( stderr, "\n%s %d", (CL->S)->name[used[a][0]], used[a][1]);
- ls[1][0]=used[a][0];
- pair_wise ( A,ns,ls, CL);
- ls[0][ns[0]++]=used[a][0];
- }
-
- A->nseq=(CL->S)->nseq;
- return A;
-
-}
-
-Alignment *stack_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
- {
- int a,n;
- int**group_list;
- int *n_groups;
- char dp_mode[100];
-
-
- sprintf ( dp_mode, "%s", CL->dp_mode);
- sprintf (CL->dp_mode, "gotoh_pair_wise");
-
- n_groups=vcalloc ( 2, sizeof (int));
- group_list=declare_int ( 2, A->nseq);
-
- n=A->nseq;
-
- for ( a=0; a<n; a++)ungap(A->seq_al[a]);
- for ( a=1; a<n; a++)
- {
- n_groups[0]=a;
- n_groups[1]=1;
- group_list[0][a-1]=a-1;
- group_list[1][0] =a;
-
- pair_wise ( A,n_groups, group_list, CL);
- fprintf ( stderr, "\n\t[%d]->[%d]", a,(int)strlen ( A->seq_al[0]));
- }
- fprintf (stderr, "\n");
- vfree(n_groups);
- free_int ( group_list, -1);
- sprintf (CL->dp_mode, "%s",dp_mode);
-
- return A;
- }
-Alignment *realign_aln_clust ( Alignment*A, Constraint_list *CL);
-Alignment *realign_aln_random_bipart_n ( Alignment*A, Constraint_list *CL, int n);
-Alignment *iterate_aln ( Alignment*A, int nit, Constraint_list *CL)
-{
- int it;
- int mode=1;
- int score, iscore, delta;
- fprintf ( CL->local_stderr, "Iterated Refinement: %d cycles START: score= %d\n", nit,iscore=aln2ecl_raw_score (A, CL) );
-
-
- if ( nit==-1)nit=A->nseq*2;
- if ( A->len_aln==0)A=very_fast_aln (A, A->nseq, CL);
- A=reorder_aln (A,(CL->S)->name, A->nseq);
-
- for (it=0; it< nit; it++)
- {
- //CL->local_stderr=output_completion (CL->local_stderr,it, nit,1, "");
- if (mode==0)A=realign_aln (A, CL);
- else if (mode ==1)A=realign_aln_random_bipart (A, CL);
- else if (mode ==2)A=realign_aln_clust (A, CL);
- else if (mode ==3)A=realign_aln_random_bipart_n (A, CL,2);
-
-
- score=aln2ecl_raw_score (A, CL);
- delta=iscore-score;
- fprintf (CL->local_stderr, "\n\tIteration Cycle: %d Score=%d Improvement= %d", it+1,score, delta);
- }
- fprintf ( CL->local_stderr, "\nIterated Refinement: Completed Improvement=%d\n", delta);
- return A;
-}
-
-int get_next_best (int seq, int nseq, int *used, int **dm);
-int get_next_best (int seq, int nseq, int *used, int **dm)
-{
- int a,set, d, bd, bseq;
-
- for (set=0,a=0; a< nseq; a++)
- {
- if (used[a] || seq==a)continue;
- d=dm[seq][a];
- if (set==0 || d>bd)
- {
- bseq=a;
- bd=d;
- set=1;
- }
- }
- return bseq;
-}
-Alignment * full_sorted_aln (Alignment *A, Constraint_list *CL)
-{
- int a,b;
- A=sorted_aln_seq (0, A, CL);
- print_aln(A);
- for (a=1; a<A->nseq; a++)
- {
- A=A->A=copy_aln (A, NULL);
- for (b=0; b<A->nseq; b++)ungap(A->seq_al[b]);
- A=sorted_aln_seq (a, A, CL);
- print_aln(A);
- }
- return A;
-}
-Alignment * sorted_aln (Alignment *A, Constraint_list *CL)
-{
- return sorted_aln_seq (-1, A, CL);
-}
-Alignment * sorted_aln_seq (int new_seq, Alignment *A, Constraint_list *CL)
-{
- int a, b=0, nseq;
- int *ns, **ls, **score, *used, **dm;
- int old_seq;
-
- dm=(CL->DM)->score_similarity_matrix;
- nseq=(CL->S)->nseq;
- score=declare_int (nseq, 3);
- used=vcalloc (nseq, sizeof (int));
- ls=declare_int (2, nseq);
- ns=vcalloc (2, sizeof (int));
-
-
- if ( new_seq==-1)
- {
- for (a=0; a<nseq; a++)
- {
- score[a][0]=a;
- score[a][1]=b;
- for ( b=0; b<nseq; b++)
- score[a][2]+=dm[a][b];
- }
- sort_int ( score,3, 2, 0, nseq-1);
- old_seq=new_seq=score[nseq-1][0];
- }
-
- for (a=1; a< nseq; a++)
- {
- used[new_seq]=1;
- ls[0][ns[0]++]=new_seq;
- ns[1]=1;
- ls[1][0]=get_next_best(new_seq,nseq, used,dm);
- old_seq=new_seq;
- new_seq=ls[1][0];
-
- A->score_aln=pair_wise (A, ns, ls,CL);
-
- }
- return A;
-}
-
-Alignment * ungap_aln4tree (Alignment *A);
-Alignment * ungap_aln4tree (Alignment *A)
-{
- int t, n, max_sim, sim;
- Alignment *B;
-
-
-
- n=35;
- max_sim=60;
-
- t=A->len_aln/10;
-
- B=copy_aln (A, NULL);
- B=ungap_aln_n(B, n);
- return B;
-
- sim=aln2sim (B, "idmat");
- while (B->len_aln<t && sim>max_sim && n>0)
- {
- n-=10;
- B=copy_aln (A, B);
- B=ungap_aln_n(B, n);
- sim=aln2sim (B, "idmat");
- }
- if ( B->len_aln<t && sim>max_sim)B=copy_aln (A, B);
- return B;
-}
-
-
-
-Alignment * iterative_tree_aln (Alignment *A,int n, Constraint_list *CL)
-{
- NT_node **T=NULL;
- int a;
-
- T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
- for ( a=0; a< n; a++)
- {
-
- Alignment *B;
-
- B=copy_aln (A, NULL);
- B=ungap_aln_n (B, 20);
- sprintf ( CL->distance_matrix_mode, "aln");
-
- CL->DM=cl2distance_matrix ( CL,B,NULL,NULL, 1);
- free_aln (B);
-
- degap_aln (A);
- T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
-
- tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
- }
- return A;
-}
-
-Alignment *profile_aln (Alignment *A, Constraint_list *CL)
-{
- int a,nseq,nseq2;
- int **ls, *ns;
- nseq=A->nseq;
- nseq2=2*nseq;
- ls=declare_int (2, nseq2);
- ns=vcalloc (2, sizeof (int));
-
- A=realloc_aln2(A,nseq2, A->len_aln);
- for (a=0; a< nseq; a++)
- ls[0][ns[0]++]=a;
- for ( a=0; a<nseq;a++)
- {
- sprintf (A->seq_al[a+nseq], "%s", (CL->S)->seq[a]);
- sprintf (A->name[a+nseq], "%s", (CL->S)->name[a]);
- A->order[a+nseq][0]=a;
- }
-
- ns[1]=1;
- for (a=0; a<nseq; a++)
- {
-
- ls[1][0]=a+nseq;
- A->score_aln=pair_wise (A, ns, ls,CL);
-
- ls[0][ns[0]++]=a+nseq;
- }
- for (a=0; a< nseq; a++)
- {
- sprintf (A->seq_al[a], "%s", A->seq_al[a+nseq]);
- }
- A->nseq=nseq;
- return A;
-}
-
-Alignment * iterative_aln ( Alignment*A, int n,Constraint_list *CL)
-{
- int *ns,**ls, **score, **dm;
- int a,b, nseq, max;
- ls=declare_int (2, A->nseq);
- ns=vcalloc (2, sizeof (int));
- ls[0][ns[0]++]=0;
-
-
-
-
-
- nseq=(CL->S)->nseq;
- score=declare_int (nseq,2);
- dm=(CL->DM)->score_similarity_matrix;
- for (a=0; a<nseq; a++)
- {
- score[a][0]=a;
- for ( b=0; b<nseq; b++)
- score[a][1]+=dm[a][b];
- score[a][1]/=nseq;
- }
- sort_int ( score,2, 1, 0, nseq-1);
-
- max=20;
- for (a=0; a<max; a++)
- {
- ns[0]=nseq-1;
- for (ns[0]=0,b=0; b<nseq; b++)
- if (b!=score[a][0])ls[0][ns[0]++]=b;
-
- fprintf (stderr, "[%s %s %d]",(CL->S)->name[score[a][0]],A->name[score[a][0]], score[a][1]);
- ns[1]=1;
- ls[1][0]=score[a][0];
- ungap_sub_aln ( A, ns[0], ls[0]);
- ungap_sub_aln ( A, ns[1], ls[1]);
- A->score_aln=pair_wise (A, ns, ls,CL);
- ls[0][ns[0]++]=a;
- }
-
- return A;
-}
-Alignment *simple_progressive_aln (Sequence *S, NT_node **T, Constraint_list *CL, char *mat)
-{
- int a;
- Alignment *A;
-
-
- A=seq2aln (S, NULL, RM_GAP);
-
- if ( !CL)
- {
-
- CL=declare_constraint_list (S, NULL, NULL, 0, NULL, NULL);
- sprintf ( CL->dp_mode, "myers_miller_pair_wise");
- sprintf ( CL->tree_mode, "nj");
- sprintf ( CL->distance_matrix_mode, "idscore");
- CL=choose_extension_mode ("matrix", CL);
- CL->gop=-10;
- CL->gep=-1;
- if (mat)CL->M=read_matrice (mat);
- CL->pw_parameters_set=1;
- CL->local_stderr=stderr;
- }
-
- if ( !T)T=make_tree (A, CL, CL->gop, CL->gep,S, NULL,MAXIMISE);
- for ( a=0; a< A->nseq; a++)ungap (A->seq_al[a]);
-
- tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
- A=reorder_aln ( A,A->tree_order,A->nseq);
-
- return A;
-}
-
-Alignment *very_fast_aln ( Alignment*A, int nseq, Constraint_list *CL)
-{
-char command[10000];
-char *tmp_seq;
-char *tmp_aln;
-FILE *fp;
-
-if ( CL && CL->local_stderr)fp=CL->local_stderr;
-else fp=stderr;
-
- fprintf (fp, "\n[Computation of an Approximate MSA...");
- tmp_seq= vtmpnam (NULL);
- tmp_aln= vtmpnam (NULL);
- output_fasta_seq ((tmp_seq=vtmpnam (NULL)), A);
- sprintf ( command, "t_coffee -infile=%s -special_mode quickaln -outfile=%s %s -outorder=input", tmp_seq, tmp_aln, TO_NULL_DEVICE);
- my_system ( command);
- A->nseq=0;
- A=main_read_aln (tmp_aln,A);
- fprintf (fp, "]\n");
- return A;
-}
-
-static NT_node* SNL;
-NT_node* tree_aln ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
-{
- int a;
-
-
-
- A->ibit=0;
- if ( strm ((CL->TC)->use_seqan, "NO"))
- {
- static char *tmp;
- NT_node *T;
- if (!tmp)tmp=vtmpnam(NULL);
- if ( CL && CL->dp_mode && strstr (CL->dp_mode, "collapse"))dump_constraint_list (CL, tmp, "w");
- T=local_tree_aln (LT, RT, A, nseq, CL);
-
- if ( CL && CL->dp_mode && strstr (CL->dp_mode, "collapse"))
- {
- empty_constraint_list (CL);
- undump_constraint_list (CL, tmp);
-
- }
- return T;
- }
- else return seqan_tree_aln (LT, RT, A, nseq, CL);
-
-}
-
-NT_node* seqan_tree_aln ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
- {
-
-
- Alignment *B;
-
-
- char *tree, *lib, *seq, *new_aln;
-
-
- //Output tree
- tree=vtmpnam (NULL);
- print_newick_tree (LT->parent, tree);
-
-
- //Output seq
- main_output_fasta_seq (seq=vtmpnam (NULL),B=seq2aln (CL->S,NULL,RM_GAP), NO_HEADER);
- free_aln (B);
-
- //Output lib
- new_aln=vtmpnam (NULL);
- vfclose (save_constraint_list ( CL, 0, CL->ne,lib=vtmpnam(NULL), NULL, "ascii",CL->S));
-
- fprintf (CL->local_stderr, "\n********* USE EXTERNAL ALIGNER: START:\n\tCOMMAND: %s -lib %s -seq %s -usetree %s -outfile %s\n", (CL->TC)->use_seqan,lib, seq, tree, new_aln);
- printf_system ( "%s -lib %s -seq %s -usetree %s -outfile %s", (CL->TC)->use_seqan,lib, seq, tree, new_aln);
- fprintf (CL->local_stderr, "\n********* USE EXTERNAL ALIGNER: END\n");
-
-
- main_read_aln (new_aln, A);
- return tree2ao (LT,RT, A, A->nseq, CL);
-
-
-
- }
-NT_node rec_local_tree_aln ( NT_node P, Alignment*A, Constraint_list *CL, int print);
-NT_node* local_tree_aln ( NT_node l, NT_node r, Alignment*A,int nseq, Constraint_list *CL)
-{
- int a;
- NT_node P, *NL;
- int **min=NULL;
-
- if (!r && !l) return NULL;
- else if (!r)P=l;
- else if (!l)P=r;
- else P=r->parent;
-
- fprintf ( CL->local_stderr, "\nPROGRESSIVE_ALIGNMENT [Tree Based]\n");
- for ( a=0; a<nseq; a++)fprintf (CL->local_stderr,"Group %4d: %s\n",a+1, A->name[a]);
- //1: make sure the Alignment and the Sequences are labeled the same way
- if (CL->translation)vfree (CL->translation);
- CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
- for ( a=0; a< (CL->S)->nseq; a++)
- CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
- A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
- A->nseq=(CL->S)->nseq;
-
- //2 Make sure the tree is in the same order
- recode_tree (P, (CL->S));
- index_tree_node(P);
- initialize_scoring_scheme (CL);
-
- if ( get_nproc()>1 && strstr (CL->multi_thread, "msa"))
- {
- int max_fork;
-
- max_fork=get_nproc()/2;//number of nodes forked, one node =>two jobs
- tree2nnode (P);
- NL=tree2node_list (P, NULL);
- min=declare_int (P->node+1,3);
- for (a=0; a<=P->node; a++)
- {
- NT_node N;
- N=NL[a];
- min[a][0]=a;
- if (!N);
- else if (N && N->nseq==1)min[a][1]=0;
- else
- {
- min[a][1]=MIN(((N->left)->nseq),((N->right)->nseq))*A->nseq+MAX(((N->left)->nseq),((N->right)->nseq));//sort on min and break ties on max
- min[a][2]=MIN(((N->left)->nseq),((N->right)->nseq));
- }
- }
- sort_int_inv (min,3, 1, 0, P->node);
- for (a=0; a<=P->node && a<max_fork; a++)
- {
- if (min[a][2]>1)(NL[min[a][0]])->fork=1;
- }
- }
- free_int (min, -1);
- rec_local_tree_aln (P, A,CL, 1);
- for (a=0; a<P->nseq; a++)sprintf (A->tree_order[a], "%s", (CL->S)->name[P->lseq[a]]);
- A->len_aln=strlen (A->seq_al[0]);
-
- fprintf ( CL->local_stderr, "\n\n");
-
- return NULL;
-}
-
-NT_node rec_local_tree_aln ( NT_node P, Alignment*A, Constraint_list *CL,int print)
-{
- NT_node R,L;
- int score;
-
-
- if (!P || P->nseq==1) return NULL;
- R=P->right;L=P->left;
-
- if (P->fork )
- {
- int s, pid1, pid2;
- char *tmp1, *tmp2;
- tmp1=vtmpnam (NULL);
- tmp2=vtmpnam (NULL);
-
- pid1=vvfork(NULL);
- if (pid1==0)
- {
- if (print==1)
- if (L->nseq>R->nseq)print=1;
-
- initiate_vtmpnam (NULL);
- rec_local_tree_aln (L, A, CL, print);
- dump_msa (tmp1,A, L->nseq, L->lseq);
- myexit (EXIT_SUCCESS);
- }
- else
- {
- pid2=vvfork(NULL);
- if (pid2==0)
- {
- if (print==1)
- if (L->nseq>R->nseq)print=0;
-
-
- initiate_vtmpnam (NULL);
- rec_local_tree_aln (R, A, CL, print);
- dump_msa (tmp2, A, R->nseq, R->lseq);
- myexit (EXIT_SUCCESS);
- }
- }
- vwaitpid (pid1, &s, 0);
- vwaitpid (pid2, &s, 0);
-
-
- undump_msa (A,tmp1);
- undump_msa (A,tmp2);
- }
- else
- {
- rec_local_tree_aln (L, A, CL, print);
- rec_local_tree_aln (R, A, CL, print);
- }
-
- P->score=A->score_aln=score=profile_pair_wise (A,L->nseq, L->lseq,R->nseq,R->lseq,CL);
- score=node2sub_aln_score (A, CL, CL->evaluate_mode,P);
- A->len_aln=strlen (A->seq_al[P->lseq[0]]);
-
- if (print)fprintf(CL->local_stderr, "\n\tGroup %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->[Len=%5d][PID:%d]%s",P->index,R->index,R->nseq,L->index,L->nseq, A->len_aln,getpid(),(P->fork==1)?"[Forked]":"" );
-
- return P;
-}
-
-
-
-NT_node* tree2ao ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
- {
- int *n_s;
- int ** l_s;
- int a, b;
-
- static int n_groups_done, do_split=0;
- int nseq2align=0;
- int *translation;
-
-
- NT_node P=NULL;
-
-
-
-
- if (n_groups_done==0)
- {
- if (SNL)vfree(SNL);
- SNL=vcalloc ( (CL->S)->nseq, sizeof (NT_node));
-
- if (CL->translation)vfree(CL->translation);
- CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
-
- for ( a=0; a< (CL->S)->nseq; a++)
- CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
-
- n_groups_done=(CL->S)->nseq;
- A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
- A->nseq=nseq;
- }
-
- translation=CL->translation;
- n_s=vcalloc (2, sizeof ( int));
- l_s=declare_int ( 2, nseq);
-
-
- if ( RT->parent !=LT->parent)fprintf ( stderr, "Tree Pb [FATAL:%s]", PROGRAM);
- else P=RT->parent;
-
- if ( LT->leaf==1 && RT->leaf==0)
- tree2ao ( RT->left, RT->right,A, nseq,CL);
-
- else if ( RT->leaf==1 && LT->leaf==0)
- tree2ao ( LT->left, LT->right,A,nseq,CL);
-
- else if (RT->leaf==0 && LT->leaf==0)
- {
- tree2ao ( LT->left, LT->right,A,nseq,CL);
- tree2ao ( RT->left, RT->right,A,nseq,CL);
- }
-
- if ( LT->leaf==1 && RT->leaf==1)
- {
- /*1 Identify the two groups of sequences to align*/
-
- nseq2align=LT->nseq+RT->nseq;
- n_s[0]=LT->nseq;
- for ( a=0; a< LT->nseq; a++)l_s[0][a]=translation[LT->lseq[a]];
- if ( LT->nseq==1)LT->group=l_s[0][0];
-
- n_s[1]=RT->nseq;
- for ( a=0; a< RT->nseq; a++)l_s[1][a]=translation[RT->lseq[a]];
- if ( RT->nseq==1)RT->group=l_s[1][0];
-
-
- P->group=n_groups_done++;
-
- if (nseq2align==nseq)
- {
- for (b=0, a=0; a< n_s[0]; a++, b++)sprintf ( A->tree_order[b],"%s", (CL->S)->name[l_s[0][a]]);
- for (a=0; a< n_s[1] ; a++, b++)sprintf ( A->tree_order[b], "%s",(CL->S)->name[l_s[1][a]]);
- n_groups_done=0;
- }
- }
- if (P->parent)P->leaf=1;
- if ( LT->isseq==0)LT->leaf=0;
- if ( RT->isseq==0)RT->leaf=0;
-
- if (RT->isseq){SNL[translation[RT->lseq[0]]]=RT;RT->score=100;}
- if (LT->isseq){SNL[translation[LT->lseq[0]]]=LT;LT->score=100;}
-
- do_split=split_condition (nseq2align,A->score_aln,CL);
- if (CL->split && do_split)
- {
-
- for (a=0; a< P->nseq; a++)SNL[CL->translation[P->lseq[a]]]=NULL;
- SNL[CL->translation[RT->lseq[0]]]=P;
-
- }
-
- vfree ( n_s);
- free_int ( l_s, 2);
- return SNL;
-
- }
-
-NT_node* tree_realn ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
- {
- int *n_s;
- int ** l_s;
- int a, b;
- int score;
- static int n_groups_done;
- int nseq2align=0;
- int *translation;
-
-
- NT_node P=NULL;
-
-
-
-
- if (n_groups_done==0)
- {
- if (SNL)vfree(SNL);
- SNL=vcalloc ( (CL->S)->nseq, sizeof (NT_node));
-
- if (CL->translation)vfree(CL->translation);
- CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
-
- for ( a=0; a< (CL->S)->nseq; a++)
- CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
- if (nseq>2)fprintf ( CL->local_stderr, "\nPROGRESSIVE_ALIGNMENT [Tree Based]\n");
- else fprintf ( CL->local_stderr, "\nPAIRWISE_ALIGNMENT [No Tree]\n");
- n_groups_done=(CL->S)->nseq;
- A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
- A->nseq=nseq;
- }
-
- translation=CL->translation;
- n_s=vcalloc (2, sizeof ( int));
- l_s=declare_int ( 2, nseq);
-
-
- if ( nseq==2)
- {
- n_s[0]=n_s[1]=1;
- l_s[0][0]=name_is_in_list ((CL->S)->name[0],(CL->S)->name, (CL->S)->nseq, MAXNAMES);
- l_s[1][0]=name_is_in_list ((CL->S)->name[1],(CL->S)->name, (CL->S)->nseq, MAXNAMES);
- A->score_aln=score=pair_wise (A, n_s, l_s,CL);
-
- vfree ( n_s);
- free_int ( l_s, 2);
- return SNL;
- }
- else
- {
- if ( RT->parent !=LT->parent)fprintf ( stderr, "Tree Pb [FATAL:%s]", PROGRAM);
- else P=RT->parent;
-
- if ( LT->leaf==1 && RT->leaf==0)
- tree_realn ( RT->left, RT->right,A, nseq,CL);
-
- else if ( RT->leaf==1 && LT->leaf==0)
- tree_realn ( LT->left, LT->right,A,nseq,CL);
-
- else if (RT->leaf==0 && LT->leaf==0)
- {
- tree_realn ( LT->left, LT->right,A,nseq,CL);
- tree_realn ( RT->left, RT->right,A,nseq,CL);
- }
-
- if ( LT->leaf==1 && RT->leaf==1 && (RT->nseq+LT->nseq)<nseq)
- {
- /*1 Identify the two groups of sequences to align*/
- int *list, s, id1, id2;
- list=vcalloc (nseq, sizeof (int));
- for (a=0; a<LT->nseq; a++)
- {
- s=translation[LT->lseq[a]];
- list[s]=1;
- }
- for (a=0; a<RT->nseq; a++)
- {
- s=translation[RT->lseq[a]];
- list[s]=1;
- }
- for (a=0; a<nseq; a++)
- {
- s=list[a];
- l_s[s][n_s[s]++]=a;
- }
-
- vfree (list);
-
- id1=sub_aln2sim (A, n_s, l_s, "idmat_sim");
-
-
- ungap_sub_aln (A, n_s[0],l_s[0]);
- ungap_sub_aln (A, n_s[1],l_s[1]);
- P->score=A->score_aln=score=pair_wise (A, n_s, l_s,CL);
- id2=sub_aln2sim (A, n_s, l_s, "idmat_sim");
-
-
-
-
- if (nseq2align==nseq)
- {
- for (b=0, a=0; a< n_s[0]; a++, b++)sprintf ( A->tree_order[b],"%s", (CL->S)->name[l_s[0][a]]);
- for (a=0; a< n_s[1] ; a++, b++)sprintf ( A->tree_order[b], "%s",(CL->S)->name[l_s[1][a]]);
- n_groups_done=0;
- }
- }
- if (P->parent)P->leaf=1;
- //Recycle the tree
- if ( LT->isseq==0)LT->leaf=0;
- if ( RT->isseq==0)RT->leaf=0;
-
- if (RT->isseq){SNL[translation[RT->lseq[0]]]=RT;RT->score=100;}
- if (LT->isseq){SNL[translation[LT->lseq[0]]]=LT;LT->score=100;}
-
- vfree ( n_s);
- free_int ( l_s, 2);
- return SNL;
- }
-
-
- }
-
-
-
-Alignment* profile_tree_aln ( NT_node P,Alignment*A,Constraint_list *CL, int threshold)
-{
- int *ns, **ls, a, sim;
- NT_node LT, RT, D, UD;
- Alignment *F;
- static NT_node R;
- static int n_groups_done;
-
-
- //first pass
- //Sequences must be in the same order as the tree sequences
- if (!P->parent)
- {
- R=P;
- n_groups_done=P->nseq+1;
- }
-
- LT=P->left;
- RT=P->right;
-
- if (LT->leaf==0)A=delayed_tree_aln1 (LT, A,CL, threshold);
- if (RT->leaf==0)A=delayed_tree_aln1 (RT, A,CL, threshold);
-
- ns=vcalloc (2, sizeof (int));
- ls=declare_int ( 2,R->nseq);
-
- if ( LT->nseq==1)
- {
- ls[0][ns[0]++]=LT->lseq[0];
- LT->group=ls[0][0]+1;
- }
- else
- node2seq_list (LT,&ns[0], ls[0]);
-
- if ( RT->nseq==1)
- {
- ls[1][ns[1]++]=RT->lseq[0];
- RT->group=ls[1][0]+1;
- }
- else
- node2seq_list (RT,&ns[1], ls[1]);
-
-
- P->group=++n_groups_done;
- fprintf (CL->local_stderr, "\n\tGroup %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,RT->group, ns[1],LT->group, ns[0]);
-
- P->score=A->score_aln=pair_wise (A, ns, ls,CL);
- sim=sub_aln2sim(A, ns, ls, "idmat_sim1");
-
- if ( sim<threshold)
- {
- UD=(ns[0]<=ns[1])?RT:LT;
- D= (ns[0]<=ns[1])?LT:RT;
-
- UD->aligned=1;
- D->aligned=0;
-
- fprintf (CL->local_stderr, "[Delayed (Sim=%4d). Kept Group %4d]",sim,UD->group);
-
-
- ungap_sub_aln (A, ns[0],ls[0]);
- ungap_sub_aln (A, ns[1],ls[1]);
- A->nseq=MAX(ns[0],ns[1]);
-
- F=A;
- while (F->A)F=F->A;
- F->A=main_read_aln (output_fasta_sub_aln (NULL, A, ns[(D==LT)?0:1], ls[(D==LT)?0:1]), NULL);
- if ( P==R)
- {
- F=F->A;
- F->A=main_read_aln (output_fasta_sub_aln (NULL, A, ns[(D==LT)?1:0], ls[(D==LT)?1:0]), NULL);
- }
- if (F->A==NULL)
- {
- printf_exit (EXIT_FAILURE, stderr, "\nError: Empty group");
- }
- }
- else
- {
- LT->aligned=1; RT->aligned=1;
- fprintf (CL->local_stderr, "[Score=%4d][Len=%5d]",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), (int)strlen ( A->seq_al[ls[0][0]]));
- A->nseq=ns[0]+ns[1];
- if (P==R)
- {
- F=A;
- while (F->A)F=F->A;
- F->A=main_read_aln (output_fasta_sub_aln2 (NULL, A, ns, ls), NULL);
- }
- }
- P->nseq=0;
- for (a=0; a<LT->nseq;a++)P->lseq[P->nseq++]=LT->lseq[a];
- for (a=0; a<RT->nseq;a++)P->lseq[P->nseq++]=RT->lseq[a];
-
- P->aligned=1;
-
- vfree ( ns);
- free_int ( ls,-1);
- return A;
-}////////////////////////////////////////////////////////////////////////////////////////
-//
-// Frame Tree Aln
-//
-////////////////////////////////////////////////////////////////////////////////////////
-
-//Alignment *frame_tree_aln (Alignment *A, Constraint_list *CL)
-//{
-
-
-////////////////////////////////////////////////////////////////////////////////////////
-//
-// Delayed Tree Aln
-//
-////////////////////////////////////////////////////////////////////////////////////////
-int delayed_pair_wise (Alignment *A, int *ns, int **ls,Constraint_list *CL);
-NT_node* delayed_tree_aln_mode1 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL);
-NT_node* delayed_tree_aln_mode2 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL);
-int paint_nodes2aligned ( NT_node P,char **list, int n);
-
-int reset_visited_nodes ( NT_node P);
-int reset_visited_nodes2 ( NT_node P);
-Alignment * make_delayed_tree_aln (Alignment *A,int n, Constraint_list *CL)
-{
- NT_node **T=NULL;
- int a;
-
- T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- delayed_tree_aln_mode1 ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
-
- for ( a=0; a< n; a++)
- {
-
- sprintf ( CL->distance_matrix_mode, "aln");
- CL->DM=cl2distance_matrix ( CL,A,NULL,NULL, 1);
- degap_aln (A);
- T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- delayed_tree_aln_mode1 ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
- }
-
- return A;
-}
-NT_node* delayed_tree_aln_mode1 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
-{
- NT_node P;
-
-
-
- P=LT->parent;P->nseq=nseq;
- paint_nodes2aligned (P, NULL, 0);
-
- A=delayed_tree_aln1 (P, A, CL,50);
- A=delayed_tree_aln2 (P, A, CL, 0);
- return NULL;
-}
-
-NT_node* delayed_tree_aln_mode2 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
-{
- NT_node P;
-
- int thr=50;
-
- P=LT->parent;P->nseq=nseq;
-
- A=delayed_tree_aln1 (P, A, CL,thr);
- thr-=10;
- while (thr>=0)
- {
- A=delayed_tree_aln2 (P, A, CL, thr);
- thr-=10;
- }
- return NULL;
-}
-
-Alignment* delayed_tree_aln1 ( NT_node P,Alignment*A,Constraint_list *CL, int threshold)
-{
- int *ns, **ls, a, sim;
- NT_node LT, RT, D, UD;
-
- static NT_node R;
- static int n_groups_done;
-
-
- //first pass
- //Sequences must be in the same order as the tree sequences
- if (!P->parent)
- {
- R=P;
- n_groups_done=P->nseq+1;
- }
-
- LT=P->left;
- RT=P->right;
-
- if (LT->leaf==0)A=delayed_tree_aln1 (LT, A,CL, threshold);
- if (RT->leaf==0)A=delayed_tree_aln1 (RT, A,CL, threshold);
-
- ns=vcalloc (2, sizeof (int));
- ls=declare_int ( 2,R->nseq);
-
-
- node2seq_list (LT,&ns[0], ls[0]);
- if ( LT->nseq==1)LT->group=LT->lseq[0]+1;
-
- node2seq_list (RT,&ns[1], ls[1]);
- if ( RT->nseq==1)RT->group=RT->lseq[0]+1;
-
-
- P->group=++n_groups_done;
-
-
- if ( ns[0]==0 || ns[1]==0)
- {
- fprintf (CL->local_stderr, "\n\tF-Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->Skipped",P->group,RT->group, ns[1],LT->group, ns[0]);
-
- LT->aligned=(ns[0]==0)?0:1;
- RT->aligned=(ns[1]==0)?0:1;
- }
- else
- {
- fprintf (CL->local_stderr, "\n\tF-Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,RT->group, ns[1],LT->group, ns[0]);
- P->score=A->score_aln=pair_wise (A, ns, ls,CL);
- sim=sub_aln2max_sim(A, ns, ls, "idmat_sim1");
-
-
- if ( sim<threshold)
- {
- UD=(ns[0]<=ns[1])?RT:LT;
- D= (ns[0]<=ns[1])?LT:RT;
-
- UD->aligned=1;
- D->aligned=0;
-
- fprintf (CL->local_stderr, "[Delayed (Sim=%4d). Kept Group %4d]",sim,UD->group);
-
- ungap_sub_aln (A, ns[0],ls[0]);
- ungap_sub_aln (A, ns[1],ls[1]);
- A->nseq=MAX(ns[0],ns[1]);
- }
- else
- {
- LT->aligned=1; RT->aligned=1;
- fprintf (CL->local_stderr, "[Score=%4d][Len=%5d]",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), (int)strlen ( A->seq_al[ls[0][0]]));
- A->nseq=ns[0]+ns[1];
- }
- P->nseq=0;
- for (a=0; a<LT->nseq;a++)P->lseq[P->nseq++]=LT->lseq[a];
- for (a=0; a<RT->nseq;a++)P->lseq[P->nseq++]=RT->lseq[a];
-
- P->aligned=1;
- }
- vfree ( ns);
- free_int ( ls,-1);
- return A;
-}
-
-Alignment* delayed_tree_aln2 ( NT_node P,Alignment*A,Constraint_list *CL, int thr)
-{
-
- NT_node LT, RT, D;
-
- static NT_node R;
-
-
- LT=P->left;
- RT=P->right;
- if (!P->parent)
- {
- R=P;
- fprintf (CL->local_stderr, "\n");
- }
- if (!LT->aligned && !RT->aligned)
- {
- printf_exit (EXIT_FAILURE, stderr, "ERROR: Unresolved Node On Groups %d [FATAL:%s]\n", P->group,PROGRAM);
- }
- else if (!LT->aligned || !RT->aligned)
- {
- int *ns, **ls, sim;
- ns=vcalloc (2, sizeof (int));
- ls=declare_int (2, R->nseq);
-
- node2seq_list (R,&ns[0], ls[0]);
-
- D=(!LT->aligned)?LT:RT;
- D->aligned=1;
- node2seq_list (D,&ns[1], ls[1]);
-
- fprintf (CL->local_stderr, "\tS-Delayed Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,D->group, ns[1],R->group, ns[0]);
- P->score=A->score_aln=pair_wise (A, ns, ls,CL);
- sim=sub_aln2max_sim(A, ns, ls, "idmat_sim1");
- if (sim<thr)
- {
- fprintf (CL->local_stderr, " [Further Delayed]\n");
- ungap_sub_aln (A, ns[0],ls[0]);
- ungap_sub_aln (A, ns[1],ls[1]);
- D->aligned=0;
- }
- else
- {
- fprintf (CL->local_stderr, "[Score=%4d][Len=%5d][thr=%d]\n",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), (int)strlen ( A->seq_al[ls[0][0]]), thr);
- D->aligned=1;
- }
- vfree (ns);free_int (ls, -1);
- }
- else
- {
- ;
- }
-
- if (LT->leaf==0)A=delayed_tree_aln2 (LT, A,CL, thr);
- if (RT->leaf==0)A=delayed_tree_aln2 (RT, A,CL, thr);
-
- return A;
-}
-
-int delayed_pair_wise (Alignment *A, int *ns, int **ls,Constraint_list *CL)
-{
- int s,s1, s2, a, b;
- int **sim;
-
-
- pair_wise (A, ns, ls, CL);
-
- sim=fast_aln2sim_list (A, "sim3", ns, ls);
-
- sort_int_inv ( sim,3, 2,0, ns[0]*ns[1]-1);
-
- for (a=0; a< 2; a++)
- for ( b=0; b< ns[a]; b++)
- A->order[ls[a][b]][4]=-1;
-
- for (a=0; a< 10 && sim[a][0]!=-1; a++)
- {
- s1=sim[a][0];
- s2=sim[a][1];
- A->order[s1][4]=0;
- A->order[s2][4]=0;
- }
-
- ungap_sub_aln (A, ns[0],ls[0]);
- ungap_sub_aln (A, ns[1],ls[1]);
-
- s=pair_wise (A, ns, ls, CL);
-
- for (a=0; a< 2; a++)
- for ( b=0; b< ns[a]; b++)
- A->order[ls[a][b]][4]=0;
-
- free_int (sim, -1);
- return s;
-}
-
-int node2seq_list2 (NT_node P, int *ns, int *ls)
-{
-
- if ( !P || P->visited ) return ns[0];
- else P->visited=1;
-
- if ( P->isseq)
- {
- ls[ns[0]++]=P->lseq[0];
- }
-
- if (P->left && (P->left) ->aligned)node2seq_list2 (P->left, ns,ls);
- if (P->right && (P->right)->aligned)node2seq_list2 (P->right,ns,ls);
- if (P->aligned && P->parent)node2seq_list2 (P->parent,ns,ls);
-
-
- return ns[0];
-}
-
-int node2seq_list (NT_node P, int *ns, int *ls)
-{
-
- if ( P->isseq && P->aligned)
- {
- ls[ns[0]++]=P->lseq[0];
- }
- else
- {
- if (P->left && (P->left) ->aligned)node2seq_list (P->left, ns,ls);
- if (P->right && (P->right)->aligned)node2seq_list (P->right,ns,ls);
- }
- return ns[0];
-}
-int paint_nodes2aligned ( NT_node P,char **list, int n)
-{
- int r=0;
- if ( P->leaf)
- {
- if ( list==NULL)
- P->aligned=1;
- else if ( name_is_in_list ( P->name, list, n, 100)!=-1)
- P->aligned=1;
- else
- P->aligned=0;
- return P->aligned;
- }
- else
- {
- r+=paint_nodes2aligned (P->left, list, n);
- r+=paint_nodes2aligned (P->right, list, n);
- }
- return r;
-}
-
-int reset_visited_nodes ( NT_node P)
-{
- while (P->parent)P=P->parent;
- return reset_visited_nodes2 (P);
-}
-int reset_visited_nodes2 ( NT_node P)
-{
- int r=0;
- if (P->left)r+=reset_visited_nodes2(P->left);
- if (P->right)r+=reset_visited_nodes2(P->right);
- r+=P->visited;
- P->visited=0;
- return r;
-}
-
-////////////////////////////////////////////////////////////////////////////////////////
-//
-// DPA_MSA
-//
-////////////////////////////////////////////////////////////////////////////////////////
-
-Alignment* dpa_msa2 ( NT_node P,Alignment*A,Constraint_list *CL);
-Alignment *dpa_align_node (NT_node P,Alignment*A,Constraint_list *CL);
-char *node2profile_list (NT_node P,Alignment*A,Constraint_list *CL, char *list);
-char * output_node_aln (NT_node P, Alignment *A, char *name);
-int node2nleaf ( NT_node P);
-
-Alignment* dpa_aln (Alignment*A,Constraint_list *CL)
-{
- NT_node P;
-
- A=iterative_tree_aln (A,1, CL);
- P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- degap_aln (A);
- while (!P->leaf)
- A=dpa_msa2(P, A, CL);
- return A;
-}
-
-int node2nleaf ( NT_node P)
-{
- int n=0;
- if ( P->leaf) return 1;
- else
- {
- n+=node2nleaf ( P->left);
- n+=node2nleaf ( P->right);
- }
- return n;
-}
-Alignment* dpa_msa2 ( NT_node P,Alignment*A,Constraint_list *CL)
-{
- int maxnseq=20;
- int n, n_l, n_r;
- n=node2nleaf (P);
-
-
- if ( n>maxnseq)
- {
- n_l=node2nleaf (P->left);
- n_r=node2nleaf (P->right);
- if (n_l>n_r)
- {
- return dpa_msa2 (P->left, A, CL);
- }
- else
- {
- return dpa_msa2 (P->right, A, CL);
- }
- }
- A=dpa_align_node (P, A, CL);
- P->leaf=1;
- return A;
-}
-
-Alignment *dpa_align_node (NT_node P,Alignment*A,Constraint_list *CL)
-{
-
- char *list, *tmp_aln;
- int a, b;
- Alignment *B;
-
-
- list=vcalloc ( P->nseq*100, sizeof (char));
- list=node2profile_list (P,A, CL, list);
-
- printf_system ( "t_coffee -profile %s -outfile=%s -dp_mode gotoh_pair_wise_lgp -msa_mode iterative_tree_aln -quiet", list,tmp_aln=vtmpnam (NULL));
- B=main_read_aln (tmp_aln, NULL);
- A=realloc_aln (A, B->len_aln+1);
- for ( a=0; a< B->nseq; a++)
- if ( (b=name_is_in_list (B->name[a], A->name, A->nseq, 100))!=-1)
- sprintf (A->seq_al[b], "%s", B->seq_al[a]);
- A->len_aln=B->len_aln;
- free_aln (B);
- vfree (list);
- return A;
-}
-char *node2profile_list (NT_node P,Alignment*A,Constraint_list *CL, char *list)
-{
- if (!P->leaf)
- {
- list=node2profile_list (P->left, A, CL, list);
- list=node2profile_list (P->right, A, CL, list);
- }
- else
- {
-
- list=strcatf (list," %s", output_node_aln (P, A, NULL));
- if ( !P->isseq)P->leaf=0;
- }
- return list;
-}
-char * output_node_aln (NT_node P, Alignment *A, char *name)
-{
- FILE *fp;
- int a;
- if (name==NULL) name=vtmpnam (NULL);
- fp=vfopen (name, "w");
-
- for (a=0; a< P->nseq; a++)
- fprintf ( fp, ">%s\n%s", A->name[P->lseq[a]], A->seq_al[P->lseq[a]]);
- vfclose (fp);
- return name;
-}
-////////////////////////////////////////////////////////////////////////////////////////
-//
-// NEW_DPA_MSA
-//
-////////////////////////////////////////////////////////////////////////////////////////
-
-Alignment * new_dpa_aln (Alignment *A,Constraint_list *CL)
-{
- NT_node P;
-
- char *tmp_aln;
- char *list;
-
- A=make_delayed_tree_aln (A,1, CL);
- P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
- set_node_score (A, P, "idmat_sim");
-
-
- list=split_nodes_nseq (A,P,15, list=vcalloc (P->nseq*200, sizeof (char)));
- printf_system ( "t_coffee -profile %s -outfile=%s -dp_mode gotoh_pair_wise_lgp -msa_mode iterative_tree_aln", list,tmp_aln=vtmpnam (NULL));
- return main_read_aln (tmp_aln, NULL);
-}
-
-char *split_nodes_nseq (Alignment *A, NT_node P, int nseq, char *list)
-{
- int a,n;
-
- n=P->nseq;
- a=100;
- while ( n>=nseq)
- {
- a--;
- n=count_threshold_nodes (A, P, a);
- }
-
- return split_nodes_idmax (A, P, a,list);
-}
-char *split_nodes_idmax (Alignment *A, NT_node P, int t, char *list)
-{
- if (P->isseq || P->score>=t)
- {
- list=strcatf (list," %s", output_node_aln (P, A, NULL));
- }
- else if ( P->score<t)
- {
- list=split_nodes_idmax (A, P->left,t, list);
- list=split_nodes_idmax (A, P->right,t,list);
- }
- return list;
-}
-int count_threshold_nodes (Alignment *A, NT_node P, int t)
-{
- int s=0;
-
- if (P->isseq || P->score>=t)
- {
- s=1;
- }
- else if ( P->score<t)
- {
- s+=count_threshold_nodes (A, P->left,t);
- s+=count_threshold_nodes (A, P->right,t);
- }
-
- return s;
-}
-int set_node_score (Alignment *A, NT_node P, char *mode)
-{
- int a;
- int ns[2], *ls[2];
-
- if (P->isseq) return 0;
- for (a=0; a<2; a++)
- {
- NT_node N;
- N=(a==0)?P->left:P->right;
- ns[a]=N->nseq;
- ls[a]=N->lseq;
- }
- P->score=sub_aln2max_sim(A, ns, ls,mode);
- set_node_score (A,P->left, mode);
- set_node_score (A,P->right, mode);
- return 1;
-}
-/////////////////////////////////////////////////////////////////////////////////////////
-int split_condition (int nseq, int score, Constraint_list *CL)
-{
- int cond1=1, cond2=1;
-
-
- if ( CL->split_nseq_thres)cond1 =(nseq<=CL->split_nseq_thres)?1:0;
- if ( CL->split_score_thres)cond2=(score>=CL->split_score_thres)?1:0;
-
- return (cond1 && cond2);
-}
-int profile_pair_wise (Alignment *A, int n1, int *l1, int n2, int *l2, Constraint_list *CL)
-{
- static int *ns;
- static int **ls;
-
- if (!ns)
- {
- ns=vcalloc (2, sizeof (int));
- ls=vcalloc (2, sizeof (int*));
- }
- ns[0]=n1;
- ls[0]=l1;
- ns[1]=n2;
- ls[1]=l2;
- return pair_wise (A, ns, ls, CL);
-}
-int pair_wise (Alignment *A, int*ns, int **l_s,Constraint_list *CL )
- {
- /*
- CL->maximise
- CL->gop;
- CL->gep
- CL->TG_MODE;
- */
- int score;
-
- int glocal;
- Pwfunc function;
-
-
-
- /*Make sure evaluation functions update their cache if needed*/
- A=update_aln_random_tag (A);
-
- if (! CL->pw_parameters_set)
- {
- fprintf ( stderr, "\nERROR pw_parameters_set must be set in pair_wise [FATAL]\n" );crash("");
- }
-
-
- function=get_pair_wise_function(CL->pair_wise, CL->dp_mode,&glocal);
- if ( CL->get_dp_cost==NULL)CL->get_dp_cost=get_dp_cost;
-
- if (strlen ( A->seq_al[l_s[0][0]])==0 || strlen ( A->seq_al[l_s[1][0]])==0)
- score=empty_pair_wise ( A, ns, l_s, CL, glocal);
- else
- score=function ( A, ns, l_s, CL);
-
- return score;
- }
-
-int empty_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int glocal)
-{
- int n=0, a, b;
- int *l=NULL;
- char *string;
- int l0, l1, len;
-
- if ( glocal==GLOBAL)
- {
- l0=strlen (A->seq_al[l_s[0][0]]);
- l1=strlen (A->seq_al[l_s[1][0]]);
- len=MAX(l1,l0);
-
- if ( len==0)return 0;
- else if (l0>l1){n=ns[1];l=l_s[1];}
- else if (l0<l1){n=ns[0];l=l_s[0];}
- string=generate_null (len);
- for ( a=0; a< n; a++)
- sprintf ( A->seq_al[l[a]], "%s", string);
- A->score=A->score_aln=0;
- A->len_aln=len;
- vfree ( string);
- return 0;
- }
- else if ( glocal==LALIGN)
- {
- A->A=declare_aln (A->S);
- (A->A)->len_aln=0;
- for ( a=0; a< 2; a++)
- for ( b=0; b<ns[a]; b++)
- A->seq_al[l_s[a][b]][0]='\0';
- (A->A)->score_aln=(A->A)->score=0;
- return 0;
- }
- else return 0;
-}
-
-
-
-
-Pwfunc get_pair_wise_function (Pwfunc pw,char *dp_mode, int *glocal)
- {
- /*Returns a function and a mode (Glogal, Local...)*/
-
-
-
- int a;
- static int npw;
- static Pwfunc *pwl;
- static char **dpl;
- static int *dps;
-
- /*The first time: initialize the list of pairwse functions*/
- if ( npw==0)
- {
- pwl=vcalloc ( 100, sizeof (Pwfunc));
- dpl=declare_char (100, 100);
- dps=vcalloc ( 100, sizeof (int));
-
- pwl[npw]=fasta_cdna_pair_wise;
- sprintf (dpl[npw], "fasta_cdna_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=cfasta_cdna_pair_wise;
- sprintf (dpl[npw], "cfasta_cdna_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=idscore_pair_wise;
- sprintf (dpl[npw], "idscore_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=gotoh_pair_wise;
- sprintf (dpl[npw], "gotoh_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=gotoh_pair_wise_lgp;
- sprintf (dpl[npw], "gotoh_pair_wise_lgp");
- dps[npw]=GLOBAL;
- npw++;
-
-
- pwl[npw]=proba_pair_wise;
- sprintf (dpl[npw], "proba_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=biphasic_pair_wise;
- sprintf (dpl[npw], "biphasic_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=subop1_pair_wise;
- sprintf (dpl[npw], "subop1_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=subop2_pair_wise;
- sprintf (dpl[npw], "subop2_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=myers_miller_pair_wise;
- sprintf (dpl[npw], "myers_miller_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=test_pair_wise;
- sprintf (dpl[npw], "test_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=fasta_gotoh_pair_wise;
- sprintf (dpl[npw], "fasta_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
- pwl[npw]=cfasta_gotoh_pair_wise;
- sprintf (dpl[npw], "cfasta_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=very_fast_gotoh_pair_wise;
- sprintf (dpl[npw], "very_fast_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=gotoh_pair_wise_sw;
- sprintf (dpl[npw], "gotoh_pair_wise_sw");
- dps[npw]=LOCAL;
- npw++;
-
- pwl[npw]=cfasta_gotoh_pair_wise_sw;
- sprintf (dpl[npw], "cfasta_sw_pair_wise");
- dps[npw]=LOCAL;
- npw++;
-
- pwl[npw]=gotoh_pair_wise_lalign;
- sprintf (dpl[npw], "gotoh_pair_wise_lalign");
- dps[npw]=LALIGN;
- npw++;
-
- pwl[npw]=sim_pair_wise_lalign;
- sprintf (dpl[npw], "sim_pair_wise_lalign");
- dps[npw]=LALIGN;
- npw++;
-
- pwl[npw]=domain_pair_wise;
- sprintf (dpl[npw], "domain_pair_wise");
- dps[npw]=MOCCA;
- npw++;
-
- pwl[npw]=gotoh_pair_wise;
- sprintf (dpl[npw], "ssec_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=ktup_pair_wise;
- sprintf (dpl[npw], "ktup_pair_wise");
- dps[npw]=LOCAL;
- npw++;
-
- pwl[npw]=precomputed_pair_wise;
- sprintf (dpl[npw], "precomputed_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=myers_miller_pair_wise;
- sprintf (dpl[npw], "default");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=viterbi_pair_wise;
- sprintf (dpl[npw], "viterbi_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=viterbiL_pair_wise;
- sprintf (dpl[npw], "viterbiL_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=viterbiD_pair_wise;
- sprintf (dpl[npw], "viterbiD_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=seq_viterbi_pair_wise;
- sprintf (dpl[npw], "seq_viterbi_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=pavie_pair_wise;
- sprintf (dpl[npw], "pavie_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=glocal_pair_wise;
- sprintf (dpl[npw], "glocal_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=linked_pair_wise;
- sprintf (dpl[npw], "linked_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
-
- pwl[npw]=linked_pair_wise_collapse;
- sprintf (dpl[npw], "linked_pair_wise_collapse");
- dps[npw]=GLOBAL;
- npw++;
-
-
- /*
- pwl[npw]=viterbiDGL_pair_wise;
- sprintf (dpl[npw], "viterbiDGL_pair_wise");
- dps[npw]=GLOBAL;
- npw++;
- */
- }
-
- for ( a=0; a< npw; a++)
- {
- if ( (dp_mode && strm (dpl[a], dp_mode)) || pwl[a]==pw)
- {
- pw=pwl[a];
- if (dp_mode)sprintf (dp_mode,"%s", dpl[a]);
- glocal[0]=dps[a];
- return pw;
- }
- }
- fprintf ( stderr, "\n[%s] is an unknown mode for dp_mode[FATAL]\n", dp_mode);
- crash ( "\n");
- return NULL;
- }
-
-
-/*******************************************************************************/
-/* */
-/* */
-/* Util Functions */
-/* */
-/* */
-/*******************************************************************************/
-
-char *build_consensus ( char *seq1, char *seq2, char *dp_mode)
- {
- Alignment *A;
- char *buf;
- int a;
- char c1, c2;
- static char *mat;
-
-
- if ( !mat) mat=vcalloc ( STRING, sizeof (char));
-
-
- A=align_two_sequences (seq1, seq2, strcpy(mat,"idmat"), 0, 0,dp_mode);
- buf=vcalloc ( A->len_aln+1, sizeof (char));
-
- for ( a=0; a< A->len_aln; a++)
- {
- c1=A->seq_al[0][a];
- c2=A->seq_al[1][a];
- if (is_gap(c1) && is_gap(c2))buf[a]='-';
- else if (is_gap(c1))buf[a]=c2;
- else if (is_gap(c2))buf[a]=c1;
- else if (c1!=c2){vfree (buf);buf=NULL;free_aln(A);return NULL;}
- else buf[a]=c1;
- }
- buf[a]='\0';
- free_sequence (free_aln (A), -1);
- return buf;
- }
-
-
-
-#ifdef FASTAL
-
-combine_profile () Comobine two profiles into one, using the edit sequence produce by the DP
- edit_sequence () insert the gaps using the
-
-int fastal (int argv, char **arg)
-{
- Sequence *S;
- int a, b;
- SeqHasch *H=NULL;
- int ktup=2;
-
- S=get_fasta_sequence (arg[1], NULL);
-
-
-
- for (a=0; a<S->nseq-1; a++)
- {
- for (b=a+1; b<S->nseq; b++)
- {
- mat[b][a]=mat[a][b]addrand()%100;
- }
- }
-
- int_dist2nj_tree (s, S->name, S->nseq, tree_name);
- T=main_read_tree (BT);
- =fastal_tree_aln (T->L,T->R,S);
-}
-
-
-
-
-NT_node fastal_tree_aln ( NT_node P, Sequence *S)
-{
- int score;
-
-
- if (!P || P->nseq==1) return NULL;
- R=P->right;L=P->left;
-
- fastal_tree_aln (P->left,S);
- fastal_tree_aln (P->right,S);
- fastal_pair_wise (P);
- return P;
-}
-
-
-NT_node fastal_pair_wise (NT_node P)
-{
- //X- 1
- //-X 2
- //XX 3
- //-- 4
-
- tb=fastal_align_profile ((P->right)->prf, (P->left)->prf);
-
- l=strlen (tb);
- for (a=0; a< l; a++)
- {
- pr1=pr2=0;
- if (tb[a]== 1 || tb[a] ==3)pr1=1;
- if (tb[a]== 2 || tb[a] ==3)pr2=1;
-
- for (b=0; b<20; b++)
- P->prf[a][b]=((pr1==1)?(P->right)->prf[ppr1][b]:0) + ((pr2==1)?(P->left)->prf[ppr2][b]:0);
- ppr1+=pr1;
- ppr2+=pr2;
- }
- free_int ((P->left)->prf, -1);
- free_int ((P->right)->prf, -1);
-}
-#endif
-/******************************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