+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdarg.h>
-#include <ctype.h>
-#include <string.h>
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "define_header.h"
-#include "dp_lib_header.h"
-void print_atom ( Atom*A);
-
-float **** quantile_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL,Pdb_param *PP, FILE *fp);
-float **** irmsdmin_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL,Pdb_param *PP, FILE *fp);
-int apdb ( int argc, char *argv[])
- {
-
- Constraint_list *CL=NULL;
- Sequence *S=NULL;
- Alignment *A=NULL;
- Alignment *EA=NULL;
- Pdb_param *pdb_param;
-
- Fname *F=NULL;
- char *file_name;
- int a,c;
-
- int n_pdb;
-
-/*PARAMETERS VARIABLES*/
- int garbage;
- char *parameters;
- FILE *fp_parameters;
-
- int quiet;
- char *se_name;
- FILE *le=NULL;
-
- char **list_file;
- int n_list;
- char **struc_to_use;
- int n_struc_to_use;
-
- char *aln;
- char *repeat_seq;
- char *repeat_pdb;
-
- char *color_mode;
- char *comparison_io;
-
- int n_excluded_nb;
-
- float maximum_distance;
- float similarity_threshold;
- float md_threshold;
-
-
- int print_rapdb;
-
- char *outfile;
- char *run_name;
-
- char *apdb_outfile;
- char *cache;
-
- char **out_aln_format;
- int n_out_aln_format;
-
- char *output_res_num;
- char *local_mode;
- float filter;
- int filter_aln;
- int irmsd_graph;
- int nirmsd_graph;
- int n_template_file;
- char **template_file_list;
- char *mode;
- int prot_min_sim;
- int prot_max_sim;
- int prot_min_cov;
- int pdb_min_sim;
- int pdb_max_sim;
- int pdb_min_cov;
- int gapped;
-
-
- char *prot_blast_server;
- char *pdb_blast_server;
-
-
- char *pdb_db;
- char *prot_db;
- int min_ncol;
-
- argv=standard_initialisation (argv, &argc);
-
-/*PARAMETER PROTOTYPE: READ PARAMETER FILE */
- declare_name (parameters);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-parameters" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "R_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Read the files in the parameter file" ,\
- /*Parameter*/ ¶meters ,\
- /*Def 1*/ "NULL" ,\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- if ( parameters && parameters[0])
- {
- argv[argc]=vcalloc ( VERY_LONG_STRING, sizeof(char));
- a=0;
- fp_parameters=vfopen (parameters, "r");
- while ((c=fgetc (fp_parameters))!=EOF)argv[1][a++]=c;
- vfclose (fp_parameters);
- argv[argc][a]='\0';
- argc++;
- argv=break_list ( argv, &argc, "=:;, \n");
- }
-/*PARAMETER PROTOTYPE*/
- declare_name (se_name);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-quiet" ,\
- /*Flag*/ &quiet ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &se_name ,\
- /*Def 1*/ "stderr" ,\
- /*Def 2*/ "/dev/null" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
- le=vfopen ( se_name, "w");
- fprintf ( le, "\nPROGRAM: %s\n",argv[0]);
-
-/*PARAMETER PROTOTYPE: IN */
- list_file=declare_char ( 200, STRING);
- n_list=get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-in" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 200 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ list_file ,\
- /*Def 1*/ "",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: IN */
- struc_to_use=declare_char ( 200, STRING);
- n_struc_to_use=get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-struc_to_use" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 200 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ struc_to_use ,\
- /*Def 1*/ "",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: COMPARISON IO */
- declare_name (comparison_io);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-io_format" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 200 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &comparison_io,\
- /*Def 1*/ "hsgd0123456",\
- /*Def 2*/ "" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: ALN */
- declare_name (aln);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-aln" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &aln,\
- /*Def 1*/ "",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: ALN */
-
- declare_name (repeat_seq);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-repeat_seq" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &repeat_seq,\
- /*Def 1*/ "",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: ALN */
- declare_name (repeat_pdb);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-repeat_pdb" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &repeat_pdb,\
- /*Def 1*/ "",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: Nb to exclude */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-n_excluded_nb" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Exclude the N Nb on each side of the central residue. -1 triggers an automatic setting equal to the window size corresponding to the sphere" ,\
- /*Parameter*/ &n_excluded_nb ,\
- /*Def 1*/ "-1" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: diatances to count */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-similarity_threshold" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &similarity_threshold,\
- /*Def 1*/ "70" ,\
- /*Def 2*/ "70" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: diatances to count */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-filter" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Filter by only keeping the best quantile" ,\
- /*Parameter*/ &filter,\
- /*Def 1*/ "1.00" ,\
- /*Def 2*/ "1.00" ,\
- /*Min_value*/ "-1.00" ,\
- /*Max Value*/ "1.00" \
- );
-/*PARAMETER PROTOTYPE: diatances to count */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-filter_aln" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Lower Case For Residues Filtered Out" ,\
- /*Parameter*/ &filter_aln,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "0" ,\
- /*Max Value*/ "1" \
- );
-/*PARAMETER PROTOTYPE: diatances to count */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-irmsd_graph" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Outputs the irmsd, position/position" ,\
- /*Parameter*/ &irmsd_graph,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "0" ,\
- /*Max Value*/ "1" \
- );
-/*PARAMETER PROTOTYPE: diatances to count */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-nirmsd_graph" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Outputs the NIRMSD VS N Removed Residues Curve" ,\
- /*Parameter*/ &nirmsd_graph,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "0" ,\
- /*Max Value*/ "1" \
- );
-/*PARAMETER PROTOTYPE: -rmsd_threshold */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-md_threshold" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &md_threshold ,\
- /*Def 1*/ "1" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: -maximum distances */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-maximum_distance" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &maximum_distance ,\
- /*Def 1*/ "10" ,\
- /*Def 2*/ "10" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-
-/*PARAMETER PROTOTYPE: -print_rapdb */
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-print_rapdb" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Prints the neighborhood of each pair of aligned residues, along with the associated local score" ,\
- /*Parameter*/ &print_rapdb ,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: RUN_NAME */
- declare_name (run_name);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-run_name" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &run_name ,\
- /*Def 1*/ "default" ,\
- /*Def 2*/ "" ,\
- /*Min_value*/ "default" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: OUTFILE */
-/*PARAMETER PROTOTYPE: OUTFILE */
- declare_name ( outfile);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-outfile" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &outfile ,\
- /*Def 1*/ "no" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "default" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: OUTFILE */
- declare_name ( apdb_outfile);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-apdb_outfile" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &apdb_outfile ,\
- /*Def 1*/ "stdout" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: OUTPUT_FORMAT */
- out_aln_format=declare_char ( 200, STRING);
- n_out_aln_format=get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-output" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 200 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ out_aln_format,\
- /*Def 1*/ "score_html" ,\
- /*Def 2*/ "" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-
-
-/*PARAMETER PROTOTYPE: INFILE */
- declare_name (color_mode);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-color_mode" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &color_mode ,\
- /*Def 1*/ "apdb" ,\
- /*Def 2*/ "irmsd" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-/*PARAMETER PROTOTYPE: INFILE */
- declare_name (output_res_num);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-seqnos" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/ &output_res_num ,\
- /*Def 1*/ "off" ,\
- /*Def 2*/ "on" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- declare_name (cache);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-cache" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "use,ignore,update,local, directory name" ,\
- /*Parameter*/ &cache ,\
- /*Def 1*/ "use" ,\
- /*Def 2*/ "update" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
- declare_name (local_mode);
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-local_mode" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Mode for choosing the Neighborhood (bubble or window)\nWhen selecting window, maximum distance becomes the window 1/2 size, in residues\nWhen using sphere, maximum_distance is the sphere radius in Angstrom" ,\
- /*Parameter*/ &local_mode ,\
- /*Def 1*/ "sphere" ,\
- /*Def 2*/ "window" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-/*PARAMETER PROTOTYPE: IN */
- template_file_list=declare_char (100, STRING);
- n_template_file=get_cl_param( \
- /*argc*/ argc , \
- /*argv*/ argv , \
- /*output*/ &le ,\
- /*Name*/ "-template_file" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1000 ,\
- /*DOC*/ "List of templates file for the sequences",\
- /*Parameter*/ template_file_list , \
- /*Def 1*/ "_SELF_P_",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- /*PARAMETER PROTOTYPE: MODE */
- declare_name (mode);
- get_cl_param( \
- /*argc*/ argc , \
- /*argv*/ argv , \
- /*output*/ &le ,\
- /*Name*/ "-mode" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "S" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Mode: irmsd, ",\
- /*Parameter*/ &mode , \
- /*Def 1*/ "irmsd",\
- /*Def 2*/ "stdin" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
-
-
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-prot_min_sim" ,\
- /*Flag*/ &prot_min_sim ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Minimum similarity between a sequence and its PDB target" ,\
- /*Parameter*/ &prot_min_sim ,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "20" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- set_int_variable ("prot_min_sim", prot_min_sim);
-
-get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-prot_max_sim" ,\
- /*Flag*/ &prot_max_sim ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Maximum similarity between a sequence and its BLAST relatives" ,\
- /*Parameter*/ &prot_max_sim ,\
- /*Def 1*/ "90" ,\
- /*Def 2*/ "100" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- set_int_variable ("prot_max_sim", prot_max_sim);
-
-get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-prot_min_cov" ,\
- /*Flag*/ &prot_min_cov ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Minimum coverage of a sequence by its BLAST relatives" ,\
- /*Parameter*/ &prot_min_cov ,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "0" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-set_int_variable ("prot_min_cov", prot_min_cov);
-
-get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-pdb_min_sim" ,\
- /*Flag*/ &pdb_min_sim ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Minimum similarity between a sequence and its PDB target" ,\
- /*Parameter*/ &pdb_min_sim ,\
- /*Def 1*/ "35" ,\
- /*Def 2*/ "35" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
- set_int_variable ("pdb_min_sim", pdb_min_sim);
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-pdb_max_sim" ,\
- /*Flag*/ &pdb_max_sim ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Maximum similarity between a sequence and its PDB target" ,\
- /*Parameter*/ &pdb_max_sim ,\
- /*Def 1*/ "100" ,\
- /*Def 2*/ "0" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- set_int_variable ("pdb_max_sim", pdb_max_sim);
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-pdb_min_cov" ,\
- /*Flag*/ &pdb_min_cov ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Minimum coverage of a sequence by its PDB target" ,\
- /*Parameter*/ &pdb_min_cov ,\
- /*Def 1*/ "50" ,\
- /*Def 2*/ "25" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-set_int_variable ("pdb_min_cov", pdb_min_cov);
-
-
-
-declare_name (pdb_blast_server);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-pdb_blast_server" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/&pdb_blast_server ,\
- /*Def 1*/ "EBI" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-declare_name (prot_blast_server);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-blast" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/&prot_blast_server ,\
- /*Def 1*/ "" ,\
- /*Def 2*/ "" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- //make sure that -blast and -blast_server are both supported blast>blast_server
- if ( !prot_blast_server[0])
- {
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-blast_server" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/&prot_blast_server ,\
- /*Def 1*/ "EBI" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- }
- set_string_variable ("blast_server", prot_blast_server);
-
-
-
- declare_name (pdb_db);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-pdb_db" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "Non Redundant PDB database" ,\
- /*Parameter*/&pdb_db ,\
- /*Def 1*/ "pdb" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- set_string_variable ("pdb_db", pdb_db);
-
-
-declare_name (prot_db);
- get_cl_param(\
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-protein_db" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "W_F" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/&prot_db ,\
- /*Def 1*/ "uniprot" ,\
- /*Def 2*/ "default" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
-
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-gapped" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "ND" ,\
- /*Parameter*/&gapped ,\
- /*Def 1*/ "0" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- get_cl_param( \
- /*argc*/ argc ,\
- /*argv*/ argv ,\
- /*output*/ &le ,\
- /*Name*/ "-min_ncol" ,\
- /*Flag*/ &garbage ,\
- /*TYPE*/ "D" ,\
- /*OPTIONAL?*/ OPTIONAL ,\
- /*MAX Nval*/ 1 ,\
- /*DOC*/ "minimum number of columns (negative: fraction)" ,\
- /*Parameter*/&min_ncol ,\
- /*Def 1*/ "4" ,\
- /*Def 2*/ "1" ,\
- /*Min_value*/ "any" ,\
- /*Max Value*/ "any" \
- );
- // set the correct mode:
- if ( strm (argv[0], "trmsd"))sprintf (mode, "trmsd");
-
- set_string_variable ("prot_db", prot_db);
-
-
- if (argc==1){myexit (EXIT_SUCCESS);}
-
- if ( strm (outfile,"no"))n_out_aln_format=0;
-
- get_cl_param( argc, argv,&le, NULL,NULL,NULL,0,0,NULL);
- prepare_cache (cache);
-
-
- if (strm ( aln, ""))
- sprintf ( aln, "%s", argv[1]);
-
- if (!is_aln (aln))
- {
- printf_exit (EXIT_FAILURE, stderr, "\n\n---- ERROR: File %s must be a valid alignment [FATAL:%s-%s]\n\n",aln,argv[0], PROGRAM);
- }
-
- pdb_param=vcalloc ( 1, sizeof(Pdb_param));
-
- pdb_param->similarity_threshold=similarity_threshold;
-
- pdb_param->md_threshold=md_threshold;
- pdb_param->maximum_distance=maximum_distance;
-
- if ( n_excluded_nb>0)
- pdb_param->n_excluded_nb=n_excluded_nb;
- else if ( n_excluded_nb==-1)
- pdb_param->n_excluded_nb=(int)((float)maximum_distance/(float)1.57);
- /* Exclude all the nb within the bubble at +1, +2, +n*/
- pdb_param->print_rapdb=print_rapdb;
- pdb_param->comparison_io=comparison_io;
-
- pdb_param->local_mode=local_mode;
- pdb_param->color_mode=lower_string (color_mode);
- pdb_param->filter=filter;
- pdb_param->filter_aln=filter_aln;
- pdb_param->irmsd_graph=irmsd_graph;
- pdb_param->nirmsd_graph=nirmsd_graph;
-
- sprintf ( list_file[n_list++], "S%s", aln);
-
-
- if (!strm (repeat_seq, ""))
- {
-
- sprintf ( template_file_list[0], "%s", process_repeat (list_file[0], repeat_seq, repeat_pdb));
- fprintf ( le, "\n##Turn a repeat List into a Template File\n");
- le=display_file_content (le,template_file_list[0]);
- fprintf ( le, "\n\n");
- }
- S=read_seq_in_n_list (list_file, n_list, NULL, NULL);
-
- le=display_sequences_names ( S,le,0, 0);
-
- if ( n_template_file)
- {
- fprintf ( le, "\nLooking For Sequence Templates:\n");
- for ( a=0; a< n_template_file; a++)
- {
- fprintf ( le, "\n\tTemplate Type: [%s] Mode Or File: [%s] [Start", template_type2type_name(template_file_list[a]), template_file_list[a]);
- S=seq2template_seq(S, template_file_list[a], F);
- fprintf ( le, "]");
- }
- }
-
- if ( !strm (run_name, "default"))
- {
- F=parse_fname(run_name);
- sprintf (F->name, "%s", F->full);
- }
- else
- {
- F=parse_fname (aln);
- }
-
- for ( a=0; a< S->nseq; a++)
- {
- char *p;
-
- p=seq2T_value (S, a, "template_file", "_P_");
-
- if (p)sprintf (S->file[a], "%s",p);
- }
-
- CL=declare_constraint_list ( S,NULL, NULL, 0,NULL, NULL);
- CL->T=vcalloc (S->nseq,sizeof (Ca_trace*));
-
-
- for ( n_pdb=0,a=0; a<S->nseq; a++)
- {
- if ( !is_pdb_file ( S->file[a])){CL->T[a]=NULL;continue;}
- CL->T[a]=read_ca_trace (S->file[a], "ATOM");
- CL->T[a]=trim_ca_trace (CL->T[a], S->seq[a]);
- (CL->T[a])->pdb_param=pdb_param;
- n_pdb++;
- }
-
- A=declare_aln (S);
-
-
- A->residue_case=KEEP_CASE;
- A=main_read_aln(aln, A);
- EA=copy_aln (A, EA);
- A->CL=CL;
-
- if ( strm (apdb_outfile, "default"))
- sprintf ( apdb_outfile, "%s.apdb_result", F->name);
-
- if ( n_pdb<2)
- {
- FILE *fp;
- fp=vfopen (apdb_outfile, "w");
- fprintf (fp, "\nYour Alignment Does Not Contain Enough Sequences With a known Structure\n");
- fprintf (fp, "To Use APDB, your alignment must include at least TWO sequences with a known structure.\n");
- fprintf (fp, "These sequences must be named according to their PDB identifier, followed by the chain index (if any) ex: 1fnkA\n");
- fprintf (fp, "[FATAL:%s]\n", PROGRAM);
- vfclose (fp);
- }
- else if ( strm (mode, "irmsd"))
- {
- EA=analyse_pdb ( A, EA, apdb_outfile);
- }
- else if ( strm (mode, "msa2tree") || strm (mode, "trmsd"))
- {
- EA=msa2struc_dist ( A, EA,F->name, gapped, min_ncol);
- }
- le=display_output_filename ( le, "APDB_RESULT", "APDB_RESULT_FORMAT_01", apdb_outfile, CHECK);
-
- if ( n_pdb>=2)
- {
- declare_name (file_name);
- for ( a=0; a< n_out_aln_format; a++)
- {
- if ( strm2( outfile, "stdout", "stderr"))sprintf (file_name, "%s", outfile);
- else if ( strm (outfile, "default"))
- sprintf (file_name, "%s.%s",F->name, out_aln_format[a]);
- else
- sprintf (file_name, "%s.%s",outfile,out_aln_format[a]);
-
- output_format_aln (out_aln_format[a],A,EA,file_name);
- le=display_output_filename ( le, "MSA", out_aln_format[a], file_name, CHECK);
- }
- }
- return EXIT_SUCCESS;
- }
-
-
-
-Constraint_list * set_constraint_list4align_pdb (Constraint_list *CL,int seq, char *dp_mode, char *local_mode, char *param_file)
-{
- static Constraint_list *PWCL;
- static Pdb_param *pdb_param;
- char **x;
- int n;
-
- if ( !CL)
- {
- free_constraint_list (PWCL);
- return NULL;
- }
- else if ( !PWCL)
- {
- PWCL=declare_constraint_list ( CL->S,NULL, NULL, 0,NULL, NULL);
-
- pdb_param=vcalloc ( 1, sizeof(Pdb_param));
- pdb_param->N_ca=0;
- pdb_param->max_delta=2.0;
- pdb_param->maximum_distance=14;
- declare_name (pdb_param->local_mode);
- sprintf (pdb_param->local_mode, "%s", local_mode);
- pdb_param->scale=50;
-
- PWCL->pw_parameters_set=1;
- PWCL->S=CL->S;
- PWCL->lalign_n_top=10;
- PWCL->sw_min_dist=10;
-
- PWCL->T=vcalloc ( (PWCL->S)->nseq, sizeof (Ca_trace*));
-
- PWCL->extend_jit=0;
- PWCL->maximise=1;
- /*PWCL->gop=-40;*/
- PWCL->gop=-50;
- PWCL->gep=-20;
- sprintf (CL->matrix_for_aa_group, "vasiliky");
- PWCL->use_fragments=0;
- PWCL->ktup=0;
- PWCL->TG_MODE=1;
- }
-
-
- if ( param_file && check_file_exists ( param_file) )
- {
- if ( (x=get_parameter ( "-nca", &n, param_file))!=NULL){pdb_param->N_ca=atoi(x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-max_delta", &n, param_file))!=NULL){pdb_param->max_delta=atof(x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-maximum_distance", &n, param_file))!=NULL){pdb_param->maximum_distance=atoi(x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-local_mode", &n, param_file))!=NULL){sprintf (pdb_param->local_mode, "%s",x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-scale", &n, param_file))!=NULL){pdb_param->scale=atoi(x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-gapopen", &n, param_file))!=NULL){PWCL->gop=atoi(x[0]);free_char (x, -1);}
- if ( (x=get_parameter ( "-gapext" , &n, param_file))!=NULL){PWCL->gep=atof(x[0]);free_char (x, -1);}
-
- }
-
-
-
-
- sprintf ( PWCL->dp_mode, "%s", dp_mode);
-
- if (strm (PWCL->dp_mode, "lalign"))sprintf (PWCL->dp_mode,"sim_pair_wise_lalign");
- else if (strm (PWCL->dp_mode, "sw"))sprintf (PWCL->dp_mode,"gotoh_pair_wise_sw");
-
- local_mode=pdb_param->local_mode;
- if ( strm ( local_mode, "hasch_ca_trace_nb")) PWCL->evaluate_residue_pair=evaluate_ca_trace_nb;
- else if ( strm ( local_mode, "hasch_ca_trace_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble;
- else if ( strm ( local_mode, "hasch_ca_trace_sap1_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_sap1_bubble;
- else if ( strm ( local_mode, "hasch_ca_trace_sap2_bubble")) PWCL->evaluate_residue_pair=evaluate_ca_trace_sap2_bubble;
-
- else if ( strm ( local_mode, "hasch_ca_trace_transversal")) PWCL->evaluate_residue_pair=evaluate_ca_trace_transversal;
- else if ( strm ( local_mode, "hasch_ca_trace_bubble_2")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble_2;
- else if ( strm ( local_mode, "hasch_ca_trace_bubble_3")) PWCL->evaluate_residue_pair=evaluate_ca_trace_bubble_3;
- else if ( strm ( local_mode, "custom_pair_score_function1")) PWCL->evaluate_residue_pair=custom_pair_score_function1;
- else if ( strm ( local_mode, "custom_pair_score_function2")) PWCL->evaluate_residue_pair=custom_pair_score_function2;
- else if ( strm ( local_mode, "custom_pair_score_function3")) PWCL->evaluate_residue_pair=custom_pair_score_function3;
- else if ( strm ( local_mode, "custom_pair_score_function4")) PWCL->evaluate_residue_pair=custom_pair_score_function4;
- else if ( strm ( local_mode, "custom_pair_score_function5")) PWCL->evaluate_residue_pair=custom_pair_score_function5;
- else if ( strm ( local_mode, "custom_pair_score_function6")) PWCL->evaluate_residue_pair=custom_pair_score_function6;
- else if ( strm ( local_mode, "custom_pair_score_function7")) PWCL->evaluate_residue_pair=custom_pair_score_function7;
- else if ( strm ( local_mode, "custom_pair_score_function8")) PWCL->evaluate_residue_pair=custom_pair_score_function8;
- else if ( strm ( local_mode, "custom_pair_score_function9")) PWCL->evaluate_residue_pair=custom_pair_score_function9;
- else if ( strm ( local_mode, "custom_pair_score_function10")) PWCL->evaluate_residue_pair=custom_pair_score_function10;
-
-
- else
- {
- fprintf ( stderr, "\n%s is an unknown hasch mode, [FATAL]\n", local_mode);
- myexit (EXIT_FAILURE);
- }
-
- if ( PWCL->T[seq]);
- else
- {
- PWCL->T[seq]=read_ca_trace (is_pdb_struc((CL->S)->name[seq]), "ATOM");
- (PWCL->T[seq])->pdb_param=pdb_param;
- PWCL->T[seq]=trim_ca_trace (PWCL->T[seq], (CL->S)->seq[seq]);
- PWCL->T[seq]=hasch_ca_trace(PWCL->T[seq]);
-
- }
-
-
- return PWCL;
-}
-
-
-
-int evaluate_ca_trace_nb (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
-
- return (int)(neighborhood_match(CL, s1,r1, s2, r2, (CL->T[s1])->Chain,(CL->T[s2])->Chain ));
- }
-int evaluate_ca_trace_sap2_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
-
-
-
- return sap2_neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble );
-
- }
-int evaluate_ca_trace_sap1_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
- /*
- Function documentation: start
-
- int evaluate_ca_trace_sap1_bubble (Constraint_list *CL, int s1, int s2, int r1, int r2)
- This function evaluates the cost for matching two residues:
-
- a1 is the cost for matching the two neighborood ( bubble type), using sap
- a1: [0,+100], +100 is the best possible match.
- a2 is the residue type weight:
- min=worst substitution value
- best=best of r1/r1, r2/r2-min
-
- a2=(r1/r2 -min)/best --> a1:[0, 100]
-
- score=a1*a2-->[-inf, +10000];
- */
-
-
-
- float a1;
-
-
- a1=(int) sap1_neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble );
-
- return (int)a1;
-
-
- }
-int evaluate_ca_trace_bubble (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
- /*
- Function documentation: start
-
- int evaluate_ca_trace_bubble (Constraint_list *CL, int s1, int s2, int r1, int r2)
- This function evaluates the cost for matching two residues:
-
- a1 is the cost for matching the two neighborood ( bubble type)
- a1: [-inf,+100-scale], +100-scale is the best possible match.
-
- */
-
-
-
- float a1;
-
-
-
- a1=(int) neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble )-((CL->T[s1])->pdb_param)->scale;
-
- return a1;
-
-
- }
-int evaluate_ca_trace_transversal (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
- return (int)(transversal_match (CL, s1, r1, s2, r2, (CL->T[s1])->Transversal,(CL->T[s2])->Transversal ));
- }
-
-int evaluate_ca_trace_bubble_3 (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
- /*This Mode evaluates :
-
- 1-The Bubble
- 2-The Match of the transversal residues
- */
-
- int a1, l1;
- int a2, l2;
- int a;
-
- l1=MAX(( (CL->T[s1])->Chain )->nb[r1][0] ,((CL->T[s2])->Chain )->nb[r2][0]);
- l2=MAX(( (CL->T[s1])->Bubble)->nb[r1][0], ((CL->T[s2])->Bubble)->nb[r2][0]);
-
- a1=(int)(neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Bubble,(CL->T[s2])->Bubble ));
- a2=(int)(transversal_match (CL, s1, r1, s2, r2, (CL->T[s1])->Transversal,(CL->T[s2])->Transversal ));
-
- if ( !l1 && !l2)return 0;
- a=(a1+a2)/2;
- return a;
- }
-int evaluate_ca_trace_bubble_2 (Constraint_list *CL, int s1, int r1, int s2, int r2)
- {
- /*This Mode evaluates :
- 1-The Ca neighborhood
- 2-The Bubble
- */
-
-
- return (int)((neighborhood_match (CL, s1, r1, s2, r2, (CL->T[s1])->Chain,(CL->T[s2])->Chain )));
- }
-
-
-/*********************************************************************************************/
-/* */
-/* FUNCTIONS FOR COMPARING TWO NEIGHBORHOODS:START */
-/* */
-/*********************************************************************************************/
-float matrix_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
-
- {
- /*
- Function documentation: start
-
- float matrix_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- This function evaluates the matrix for matching two residues:
-
- min=worst substitution value
- best=best of r1/r1, r2/r2-min
-
- a2=(r1/r2 -min)/best --> a1:[0, 100]
-
- score=a1*a2-->[-inf, +10000];
- */
-
-
-
- float a2;
- float m1, m2, m;
- static float min=0;
- int a, b;
-
- if ( !CL->M)
- {
- CL->M=read_matrice ( "pam250mt");
- min=CL->M[0][0];
- for ( a=0; a< 26; a++)
- for ( b=0; b< 26; b++)min=MIN(min, CL->M[a][b]);
- }
-
- if ( r1<=0 || r2<=0)return 0;
- m1=CL->M[(CL->S)->seq[s1][r1-1]-'A'][(CL->S)->seq[s1][r1-1]-'A']-min;
- m2=CL->M[(CL->S)->seq[s2][r2-1]-'A'][(CL->S)->seq[s2][r2-1]-'A']-min;
- m=MAX(m1, m2);
- a2=(CL->M[(CL->S)->seq[s1][r1-1]-'A'][(CL->S)->seq[s2][r2-1]-'A']-min)/m;
-
- return a2;
- }
-
-
-float transversal_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- {
- int a, l1, l2;
- float score=0;
- float delta, max_delta;
- float max;
- Pdb_param*PP;
-
- PP=(CL->T[s1])->pdb_param;
- max_delta=PP->max_delta;
-
- l1=nbs1->nb[r1][0];
- l2=nbs2->nb[r2][0];
-
- if ( l1!=l2 || l1<(PP->N_ca)) return 0;
-
-
- max=MAX(l1, l2)*max_delta;
- for ( delta=0,a=0; a< l2 ; a++)
- {
-
- delta+=max_delta-FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][a]));
- }
- score=(delta*100)/max;
-
-
-
- return score;
- }
-
-float neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- {
- static float **table;
- static int table_size;
- int a, b, l1, l2;
- float score=0;
- float ins, del, sub;
- float delta, max_delta;
- float max;
- Pdb_param*PP;
-
-
- PP=(CL->T[s1])->pdb_param;
- max_delta=PP->max_delta;
-
-
- if ( r1> 0 && r2 >0) {r1--; r2--;}
- else return 0;
-
- l1=nbs1->nb[r1][0];
- l2=nbs2->nb[r2][0];
-
- if (table_size< (MAX(l1, l2)+1))
- {
- table_size=MAX(l1, l2)+1;
- if ( table)free_float (table, -1);
- table=NULL;
- }
- if ( !table) table=declare_float (table_size, table_size);
-
-
- max=MAX(l1, l2)*max_delta;
- if ( max==0)return 0;
-
-
- table[0][0]=0;
- for ( b=1; b<=l2; b++)
- {
- table[0][b]=0;
- }
- for ( a=1; a<=l1; a++)
- {
- table[a][0]=0;
- for ( b=1; b<=l2 ; b++)
- {
-
- delta=max_delta-FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][b]));
-
- del=table[a-1][b];
- ins=table[a][b-1];
- sub= table[a-1][b-1]+delta;
-
- if ( del >= ins && del >= sub)score=del;
- else if ( ins >= del && ins >= sub) score=ins;
- else score=sub;
- table[a][b]=score;
- }
- }
-
-
- score=((((score)*100)/max));
-
-
- return score;
- }
-
-float sap1_neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- {
- /*
- Function documentation: start
-
- float sap1_neighborhood_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- This function is adapted from Taylor, Orengo, Protein Structure Alignment JMB 1989, (208)1-22
- It is the first function where
- score= A/(|dra-drb|+b)
-
- Function documentation: end
- */
-
- static float **table;
- static int table_size;
- int a, b, l1, l2;
- float score=0;
- float ins, del, sub;
- float delta;
- float max;
-
- int A=50;
- int B=5;
-
-
-
-
-
-
- if ( r1> 0 && r2 >0) {r1--; r2--;}
- else return 0;
-
- l1=nbs1->nb[r1][0];
- l2=nbs2->nb[r2][0];
-
- if (table_size< (MAX(l1, l2)+1))
- {
- table_size=MAX(l1, l2)+1;
- if ( table)free_float (table, -1);
- table=NULL;
- }
- if ( !table) table=declare_float (table_size, table_size);
-
-
- max=MAX(l1, l2)*(A/B);
- if ( max==0)return 0;
-
-
- table[0][0]=0;
- for ( b=1; b<=l2; b++)
- {
- table[0][b]=0;
- }
- for ( a=1; a<=l1; a++)
- {
- table[a][0]=0;
- for ( b=1; b<=l2 ; b++)
- {
-
- delta=A/(FABS((nbs1->d_nb[r1][a]-nbs2->d_nb[r2][b]))+B);
-
- del=table[a-1][b];
- ins=table[a][b-1];
- sub= table[a-1][b-1]+delta;
- if ( del >= ins && del >= sub)score=del;
- else if ( ins >= del && ins >= sub) score=ins;
- else score=sub;
- table[a][b]=score;
- }
- }
-
-
- score=((score*100))/(max);
-
-
- return score;
- }
-
-float sap2_neighborhood_match (Constraint_list *CL, int s1, int r1, int s2, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- {
- /*
- Function documentation: start
-
- float sap1_neighborhood_match (Constraint_list *CL, int s1, int s2, int r1, int r2, Struct_nb *nbs1, Struct_nb *nbs2)
- This function is adapted from Taylor, Orengo, Protein Structure Alignment JMB 1989, (208)1-22
- It is the first function where
- score= A/(|dra-drb|+b)
-
- Function documentation: end
- */
-
- static float **table;
- static int table_size;
- int a, b, l1, l2;
- float score=0;
- float ins, del, sub;
- float delta;
- float max;
-
- Amino_acid **pep1;
- Amino_acid **pep2;
- static Atom *vX_1, *vY_1, *vZ_1;
- static Atom *vX_2, *vY_2, *vZ_2;
- static Atom *ca1, *ca2;
- float val;
-
- int A=50;
- int B=2;
-
-
-
-
- if ( r1> 0 && r2 >0) {r1--; r2--;}
- else return 0;
-
- /*Make up the referencial*/
- pep1=(CL->T[s1])->peptide_chain;
- pep2=(CL->T[s2])->peptide_chain;
-
- /*Get Referencial for CA1*/
- if ( (pep1[r1])->C)vX_1 =diff_atom(pep1[r1]->C,pep1[r1]->CA, vX_1);
- if ( (pep1[r1])->N)vY_1 =diff_atom(pep1[r1]->N,pep1[r1]->CA, vY_1);
- if ( (pep1[r1])->CB)vZ_1=diff_atom(pep1[r1]->CB,(pep1[r1])->CA,vZ_1);
- else vZ_1=add_atom (vX_1, vY_1, vZ_1);
-
-
-
-
-
- /*Get Referencial for CA2*/
- if ( (pep2[r2])->C)vX_2 =diff_atom((pep2[r2])->C,(pep2[r2])->CA, vX_2);
- if ( (pep2[r2])->N)vY_2 =diff_atom((pep2[r2])->N,(pep2[r2])->CA, vY_2);
- if ( (pep2[r2])->CB)vZ_2=diff_atom((pep2[r2])->CB,(pep2[r2])->CA, vZ_2);
- else vZ_2=add_atom (vX_2, vY_2, vZ_2);
-
-
-
-
- /*END OF GETTING REFERENCIAL*/
-
- /*Test
- if ( r1>1 && r2>1)
- {
- fprintf (stdout,"\n\t*******");
-
- fprintf (stdout, "RESIDUE %d %c", r1, (CL->S)->seq[s1][r1]);
- if ( (pep1[r1])->CA)fprintf (stdout,"\n\tCA ");print_atom (pep1[r1]->CA );
- if ( (pep1[r1])->C)fprintf (stdout,"\n\tC ");print_atom (pep1[r1]->C );
- if ( (pep1[r1])->N)fprintf (stdout,"\n\tN ");print_atom (pep1[r1]->N );
- if ( (pep1[r1])->CB)fprintf (stdout,"\n\tCB ");print_atom (pep1[r1]->CB );
- fprintf (stdout,"\n\t*******");
- fprintf (stdout,"\n\tvX ");print_atom ( vX_1);
- fprintf (stdout,"\n\tvY ");print_atom ( vY_1);
- fprintf (stdout,"\n\tvZ ");print_atom ( vZ_1);
-
- ca1= copy_atom ((pep1[r1-1])->CA, ca1);
- ca1 =diff_atom(ca1,(pep1[r1])->CA, ca1);
- fprintf (stdout,"\n\tca ");print_atom ( ca1);
- fprintf ( stdout, "\n\tSQ1=%d ", (int)square_atom(ca1));
- ca1=reframe_atom(vX_1, vY_1, vZ_1, ca1, ca1);
- fprintf ( stdout, "\n\tSQ2=%d ", (int)square_atom(ca1));
- fprintf (stdout,"\n\tca ");print_atom ( ca1);
- fprintf (stdout,"\n\n");
- }
- */
-
- l1=nbs1->nb[r1][0];
- l2=nbs2->nb[r2][0];
-
- if (table_size< (MAX(l1, l2)+1))
- {
- table_size=MAX(l1, l2)+1;
- if ( table)free_float (table, -1);
- table=NULL;
- }
- if ( !table) table=declare_float (table_size, table_size);
-
-
- max=MAX(l1, l2)*(A/B);
-
- if ( max==0)return 0;
-
-
- table[0][0]=0;
- for ( b=1; b<=l2; b++)
- {
- table[0][b]=0;
- }
-
- for ( a=1; a<=l1; a++)
- {
- ca1=copy_atom ((CL->T[s1])->structure[nbs1->nb[r1][a]], ca1);
- ca1=diff_atom(ca1,(pep1[r1])->CA, ca1);
- ca1=reframe_atom(vX_1, vY_1, vZ_1, ca1, ca1);
-
- table[a][0]=0;
- for ( b=1; b<=l2 ; b++)
- {
- ca2 =copy_atom((CL->T[s2])->structure[nbs2->nb[r2][b]], ca2);
- ca2 =diff_atom(ca2,(pep2[r2])->CA, ca2);
- ca2 =reframe_atom(vX_2, vY_2, vZ_2, ca2, ca2);
-
- ca2=diff_atom(ca2,ca1,ca2);
- val=square_atom (ca2);
-
- val=(float)sqrt ((double)val);
-
- delta=A/(val+B);
-
-
- del=table[a-1][b];
- ins=table[a][b-1];
- sub= table[a-1][b-1]+delta;
-
- if ( del >= ins && del >= sub)score=del;
- else if ( ins >= del && ins >= sub) score=ins;
- else score=sub;
- table[a][b]=score;
- }
- }
-
-
- score=(((score*100))/(max)-50);
-
-
- return score;
- }
-
-/*********************************************************************************************/
-/* */
-/* APDB */
-/* */
-/*********************************************************************************************/
-float **** irmsdmin_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL, Pdb_param *PP, FILE *fp)
-{
- int s1, s2, a,col1, n,n2=0, t,flag;
- int **pos, **list;
- float nirmsd, min_nirmsd,max_nirmsd,ref_sum, sum, sum2;
- float **normalized_len;
-
- normalized_len=declare_float (A->nseq+1, A->nseq+1);
- for (s1=0; s1<A->nseq; s1++)
- {
- int l1, l2, r1, r2, p;
- for (s2=0; s2<A->nseq; s2++)
- {
- for ( l1=l2=p=0; p< A->len_aln; p++)
- {
- r1=A->seq_al[s1][p];
- r2=A->seq_al[s2][p];
- if (!is_gap(r1) && isupper(r1))l1++;
- if (!is_gap(r2) && isupper(r2))l2++;
- }
- normalized_len[s1][s2]=MIN(l1,l2);
- }
- }
-
- pos=aln2pos_simple (A, A->nseq);
- for ( s1=0; s1< A->nseq; s1++)
- for ( s2=0; s2<A->nseq; s2++)
- {
- if ( s1==s2) continue;
- else if (!(CL->T[A->order[s1][0]]) || !(CL->T[A->order[s2][0]]))continue;
-
- list=declare_int (A->len_aln, 2);
-
- for ( sum=0,n=0,col1=0; col1< A->len_aln; col1++)
- {
- if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
- else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
- else if ( residues[s1][s2][pos[s1][col1]-1][0]==0)continue;
-
- list[n][0]=pos[s1][col1]-1;
- list[n][1]=(int)100000*residues[s1][s2][pos[s1][col1]-1][4];
- sum2+=residues[s1][s2][pos[s1][col1]-1][4];
- n++;
- }
-
- if (n==0)return residues;
-
- sort_int_inv (list, 2, 1,0, n-1);
- for (sum=0,a=0; a<n; a++)
- {
- sum+=list[a][1];
- }
- ref_sum=sum;
- nirmsd=min_nirmsd=max_nirmsd=sum/(n*n);
- t=0;
-
-
- /*1 Find the maximum*/
- sum=ref_sum;
- for (flag=0,a=0; a< n-1; a++)
- {
- sum-=list[a][1];
- nirmsd=sum/((n-(a+1))*(n-(a+1)));
- if (nirmsd<max_nirmsd)flag=1;
- if ((nirmsd>=max_nirmsd) && flag==1)break;
- n2=a;
- }
-
- sum=ref_sum;
- for (a=0; a<n2-1; a++)
- {
- sum-=list[a][1];
- nirmsd=sum/((n-(a+1))*(n-(a+1)));
-
-
- if ( nirmsd<min_nirmsd)
- {
- min_nirmsd=nirmsd;
- t=a;
- if ( PP->nirmsd_graph)
- {
- fprintf ( stdout, "\n_NIRMSD_GRAPH %s %s POS: %4d Removed: %4d NiRMSD: %.2f", A->name[s1], A->name[s2], list[a][0],a,(nirmsd/100000)*normalized_len[s1][s2]);
- }
- }
- }
-
- if ( PP->print_rapdb)
- {
- for ( a=0; a<n; a++)
- {
- if ( list[a][1]>0 && a<=t)fprintf ( stdout, "\nRAPDB QUANTILE REMOVE S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
- else if ( list[a][1]>0 && a>t)fprintf ( stdout, "\nRAPDB QUANTILE KEEP S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
- }
- }
-
- fprintf ( stdout, "\n# MINIMISATION FILTER ON: NiRMSD minimsation resulted in the removal of %d [out of %d] Columns On the alignment %s Vs %s\n", t, n, A->name[s1], A->name[s2]);
- for ( a=0; a<=t; a++)
- {
-
- residues[s1][s2][list[a][0]][0]=0;
- residues[s1][s2][list[a][0]][1]=0;
- residues[s1][s2][list[a][0]][2]=0;
- residues[s1][s2][list[a][0]][3]=0;
- residues[s1][s2][list[a][0]][4]=-1;
-
- }
-
- free_int (list, -1);
- }
- free_float (normalized_len, -1);
- return residues;
-}
-float **** quantile_apdb_filtration ( Alignment *A, float ****residues, Constraint_list *CL, Pdb_param *PP,FILE *fp)
-{
- int s1, s2, a,col1, n, t;
- int **pos, **list;
-
- pos=aln2pos_simple (A, A->nseq);
- for ( s1=0; s1< A->nseq; s1++)
- for ( s2=0; s2<A->nseq; s2++)
- {
- if ( s1==s2) continue;
- else if (!(CL->T[A->order[s1][0]]) || !(CL->T[A->order[s2][0]]))continue;
-
- list=declare_int (A->len_aln, 2);
-
- for ( n=0,col1=0; col1< A->len_aln; col1++)
- {
- if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
- else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
-
- list[n][0]=pos[s1][col1]-1;
- list[n][1]=(int)100*residues[s1][s2][pos[s1][col1]-1][4];
- n++;
-
- }
-
- sort_int_inv (list, 2, 1,0, n-1);
-
- t=quantile_rank ( list,1, n,PP->filter);
-
- if ( PP->print_rapdb)
- {
- for ( a=0; a<n; a++)
- {
- if ( list[a][1]>0 && a<t)fprintf ( stdout, "\nRAPDB QUANTILE REMOVE S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
- else if ( list[a][1]>0 && a>t)fprintf ( stdout, "\nRAPDB QUANTILE KEEP S1: %3d S2: %3d COL: %3d SCORE*100: %d", s1, s2, list[a][0], list[a][1]);
- }
- }
-
- for ( a=0; a<t; a++)
- {
-
- residues[s1][s2][list[a][0]][0]=0;
- residues[s1][s2][list[a][0]][1]=0;
- residues[s1][s2][list[a][0]][2]=0;
- residues[s1][s2][list[a][0]][3]=0;
- residues[s1][s2][list[a][0]][4]=-1;
-
- }
-
- free_int (list, -1);
- }
-
- return residues;
-}
-Alignment * analyse_pdb ( Alignment *A, Alignment *ST, char *results)
- {
- int s1,s2,r1, r2,b, p;
- int **pos;
- float **normalize_len;
- float m2, m4;
- float pair_tot=0, pair_m1, pair_m2, pair_m3, pair_m4, pair_m5, pair_len=0;
- float seq_tot, seq_m1, seq_m2, seq_m3, seq_m4, seq_m5,seq_len;
- float msa_tot, msa_m1, msa_m2, msa_m3, msa_m4, msa_m5, msa_len;
- float iRMSD_unit, iRMSD_max, iRMSD_min;
- float ****residues;
- Pdb_param *PP=NULL;
- Constraint_list *CL;
- char *average_file, *pairwise_file, *total_file, *irmsd_file=0;
- FILE *fp, *average,*pairwise, *total, *irmsd_graph=0;
-
-
- fp =vfopen ( results, "w");
- pairwise=vfopen ((pairwise_file=vtmpnam (NULL)),"w");
- average =vfopen ((average_file =vtmpnam (NULL)),"w");
- total =vfopen ((total_file =vtmpnam (NULL)),"w");
-
-
- CL=A->CL;
-
- for ( s1=0; s1< (A->S)->nseq; s1++)
- if ( CL->T[s1]){PP=(CL->T[s1])->pdb_param;break;}
-
- if (PP->irmsd_graph)irmsd_graph =vfopen ((irmsd_file =vtmpnam (NULL)),"w");
-
- fprintf ( fp, "\nAPDB_RESULT_FORMAT_02\n");
- residues=analyse_pdb_residues ( A, A->CL,PP);
- if ( PP->filter>=0)residues=quantile_apdb_filtration (A, residues, A->CL,PP, fp);
- else if ( PP->filter<0)residues=irmsdmin_apdb_filtration (A, residues, A->CL,PP, fp);
-
- pos=aln2pos_simple (A, A->nseq);
-
-
-
-
-
- /*Compute the alignment length for normalization*/
- normalize_len=declare_float (A->nseq+1, A->nseq+1);
- for (s1=0; s1<A->nseq; s1++)
- {
- int l1, l2, r1, r2;
- for (s2=0; s2<A->nseq; s2++)
- {
- for ( l1=l2=p=0; p< A->len_aln; p++)
- {
- r1=A->seq_al[s1][p];
- r2=A->seq_al[s2][p];
- if (!is_gap(r1) && isupper(r1))l1++;
- if (!is_gap(r2) && isupper(r2))l2++;
- }
- normalize_len[s1][s2]=MIN(l1,l2);
- }
- }
-
- msa_len=msa_tot=msa_m1=msa_m2=msa_m3=msa_m4=msa_m5=0;
-
- for ( s1=0; s1< A->nseq; s1++)
- {
- if ( !(CL->T[A->order[s1][0]]))continue;
- seq_len=seq_tot=seq_m1=seq_m2=seq_m3=seq_m4=seq_m5=0;
- for ( s2=0; s2< A->nseq; s2++)
- {
- if ( s1==s2)continue;
- if ( !(CL->T[A->order[s2][0]]))continue;
- pair_tot=pair_m1=pair_m2=pair_m3=pair_m4=pair_m5=0;
- for ( p=0; p< A->len_aln; p++)
- {
- r1=A->seq_al[s1][p];
- r2=A->seq_al[s2][p];
- b=pos[s1][p]-1;
-
-
- if (PP->filter_aln)
- {
- if (is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0)
- {
- A->seq_al[s1][p]=tolower(r1);
- A->seq_al[s2][p]=tolower(r2);
- }
- else
- {
- A->seq_al[s1][p]=toupper(r1);
- A->seq_al[s2][p]=toupper(r2);
- }
-
- }
-
- if ( PP->irmsd_graph && ( is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0))
- {
-
- fprintf ( irmsd_graph, "\n_IRMSD_GRAPH %10s %10s ALN: %c%c iRMSD: -1.00", A->name[s1], A->name[s2],A->seq_al[s1][p], A->seq_al[s2][p]);
- }
-
- if (is_gap(r1) || is_gap(r2) || residues[s1][s2][b][0]==0)continue;
- pair_tot++;
-
- /*APDB*/
- m2=(residues[s1][s2][b][2]*100)/residues[s1][s2][b][0];
- if (m2>PP->similarity_threshold){pair_m3++;}
-
- /*iRMSD*/
-
- m4=residues[s1][s2][b][4];
-
- if ( PP->irmsd_graph )
- {
- fprintf ( irmsd_graph, "\nIRMSD_GRAPH %10s %10s ALN: %c%c iRMSD: %.2f", A->name[s1], A->name[s2],A->seq_al[s1][p], A->seq_al[s2][p], m4);
- }
- pair_m4+=m4;
- }
- pair_len=normalize_len[s1][s2];
- if ( s1>s2)
- {
-
- fprintf ( pairwise, "\n\n#PAIRWISE: %s Vs %s",A->name[s1], A->name[s2]);
- fprintf ( pairwise, "\n\tPAIRWISE EVALUATED: %6.2f %% [%s Vs %s] ", (pair_len==0)?-1:(pair_tot*100)/pair_len,A->name[s1], A->name[s2]);
- fprintf ( pairwise, "\n\tPAIRWISE APDB: %6.2f %% [%s Vs %s] ", (pair_tot==0)?-1:(pair_m3*100)/pair_tot,A->name[s1], A->name[s2]);
- fprintf ( pairwise, "\n\tPAIRWISE iRMSD: %6.2f Angs [%s Vs %s]", (pair_tot==0)?-1:pair_m4/pair_tot,A->name[s1], A->name[s2]);
- fprintf ( pairwise, "\n\tPAIRWISE NiRMSD: %6.2f Angs [%s Vs %s] [%d pos]", (pair_tot==0)?-1:(pair_m4*pair_len)/(pair_tot*pair_tot), A->name[s1], A->name[s2], (int)pair_tot);
- fprintf ( pairwise, "\n\tRAPDB PAIRS PAIRWISE N_NONEMPTY_PAIRS %d N_MAXIMUM_PAIRS %d",(int) pair_tot, (int)pair_len);
- msa_m3+=pair_m3;
- msa_m4+=pair_m4;
- msa_tot+=pair_tot;
- msa_len+=pair_len;
- }
- seq_m3+=pair_m3;
- seq_m4+=pair_m4;
- seq_tot+=pair_tot;
- seq_len+=pair_len;
-
- }
-
- fprintf ( average, "\n\n#AVERAGE For Sequence %s", A->name[s1]);
- fprintf ( average, "\n\tAVERAGE EVALUATED: %6.2f %% [%s]", (seq_len==0)?-1:(seq_tot*100)/seq_len, A->name[s1]);
- fprintf ( average, "\n\tAVERAGE APDB: %6.2f %% [%s]", (seq_tot==0)?-1:(seq_m3*100)/seq_tot, A->name[s1]);
- fprintf ( average, "\n\tAVERAGE iRMSD: %6.2f Angs [%s]", (seq_tot==0)?-1:seq_m4/seq_tot, A->name[s1]);
- fprintf ( average, "\n\tAVERAGE NiRMS: %6.2f Angs [%s]", (seq_tot==0)?-1:(seq_m4*seq_len)/(seq_tot*seq_tot), A->name[s1]);
- if ( strm (PP->color_mode, "apdb"))ST->score_seq[s1]=(seq_tot==0)?-1:(seq_m3*100)/pair_tot;
- if (PP->print_rapdb)fprintf (average, "\n\tRAPDB PAIRS AVERAGE N_NONEMPTY_PAIRS %d N_MAXIMUM_PAIRS %d", (int)pair_tot, (int)pair_len);
-
- if ( strm (PP->color_mode, "irmsd"))ST->score_seq[s1]=(seq_tot==0)?-1:10*((seq_m4*pair_len)/(seq_tot*seq_tot));
-
- }
- fprintf ( total, "\n\n#TOTAL for the Full MSA");
- fprintf ( total, "\n\tTOTAL EVALUATED: %6.2f %% ", (msa_len==0)?-1:(msa_tot*100)/msa_len);
- fprintf ( total, "\n\tTOTAL APDB: %6.2f %% ", (msa_tot==0)?-1:(msa_m3*100)/msa_tot);
- fprintf ( total, "\n\tTOTAL iRMSD: %6.2f Angs", (msa_tot==0)?-1:msa_m4/msa_tot);
- fprintf ( total, "\n\tTOTAL NiRMSD: %6.2f Angs", (msa_tot==0)?-1:(msa_m4*msa_len)/(msa_tot*msa_tot));
- if (PP->print_rapdb)fprintf (total, "\n\tRAPDB PAIRS TOTAL N_NONEMPTY_PAIRS: %d N_MAXIMUM_PAIRS %d", (int)msa_tot, (int)msa_len);
-
- if ( strm (PP->color_mode, "apdb")) ST->score_aln=ST->score=A->score_aln=A->score=(msa_tot==0)?-1:(msa_m3*100)/msa_tot;
- if ( strm (PP->color_mode, "irmsd"))ST->score_aln=ST->score=A->score_aln=A->score=(msa_tot==0)?-1:10*((msa_m4*msa_len)/(msa_tot*msa_tot));
-
- vfclose (average);vfclose (total); vfclose (pairwise);if (PP->irmsd_graph)vfclose (irmsd_graph);
- fp=display_file_content (fp, pairwise_file);
- fp=display_file_content (fp, average_file);
- fp=display_file_content (fp, total_file);
- if ( PP->irmsd_graph)fp=display_file_content (fp, irmsd_file);
-
- fprintf ( fp, "\n\n# EVALUATED: Fraction of Pairwise Columns Evaluated\n");
- fprintf ( fp, "# APDB: Fraction of Correct Columns according to APDB\n");
- fprintf ( fp, "# iRMDS: Average iRMSD over all evaluated columns\n");
- fprintf ( fp, "# NiRMDS: iRMSD*MIN(L1,L2)/Number Evaluated Columns\n");
- fprintf ( fp, "# Main Parameter: -maximum_distance %.2f Angstrom\n", PP->maximum_distance);
-
- fprintf ( fp, "# Undefined values are set to -1 and indicate LOW Alignment Quality\n");
- fp=print_program_information (fp, NULL);
-
-
-
-
- /*Color Output*/
- for (iRMSD_max=0,iRMSD_min=10000,s1=0; s1<A->nseq; s1++)
- for ( s2=0; s2< A->nseq; s2++)
- for (p=0; p<A->len_aln; p++)
- {
- if ( residues[s1][s2][p][4]>0)
- {
- iRMSD_max=MAX(iRMSD_max, residues[s1][s2][p][4]);
- iRMSD_min=MAX(iRMSD_min, residues[s1][s2][p][4]);
- }
-
- }
- iRMSD_unit=iRMSD_max/8;
-
- for (p=0; p< A->len_aln; p++)
- for ( s1=0; s1< A->nseq; s1++)
- {
-
- for ( p=0; p< A->len_aln; p++)
- {
- r1=A->seq_al[s1][p];
- b=pos[s1][p]-1;
- if ( is_gap(r1) || !(CL->T[A->order[s1][0]]))
- ST->seq_al[s1][p]=NO_COLOR_RESIDUE;
- else
- {
- float tot_m2=0, tot_m4=0, v=0;
- seq_m2=seq_m4=0;
-
- for (s2=0; s2< A->nseq; s2++)
- {
- r2=A->seq_al[s1][p];
- if ( s1==s2) continue;
- if (is_gap(r2) || !(CL->T[A->order[s1][0]]) || residues[s1][s2][b][0]==0)continue;
-
- seq_m2+=m2=(residues[s1][s2][b][2]*100)/residues[s1][s2][b][0];
- tot_m2++;
-
- m4=residues[s1][s2][b][4];
- if (m4>=0)
- {
- seq_m4+=m4;
- tot_m4++;
- }
- }
-
- if (strm ( PP->color_mode, "apdb"))
- {
- if (tot_m2==0)v=NO_COLOR_RESIDUE;
- else v=MIN((seq_m2/(10*tot_m2)),9);
- }
- else if ( strm (PP->color_mode, "irmsd"))
- {
- if ( tot_m4==0)v=NO_COLOR_RESIDUE;
- else v=(8-(int)((seq_m4/(iRMSD_unit*tot_m4))))+1;
- }
- ST->seq_al[s1][p]=v;
-
- }
- }
- }
- for ( p=0; p<A->len_aln; p++) ST->seq_al[A->nseq][p]=NO_COLOR_RESIDUE;
-
-
- ST->generic_comment=vcalloc ( 100, sizeof (int));
- if ( strm (PP->color_mode, "apdb"))
- {
- sprintf ( ST->generic_comment, "# APDB Evaluation: Color Range Blue-[0 %% -- 100 %%]-Red\n# Sequence Score: APDB\n# Local Score: APDB\n\n");
- }
- else if ( strm (PP->color_mode, "irmsd"))
- {
- sprintf ( ST->generic_comment, "\n# iRMSD Evaluation:\n# Sequence score: NiRMSD (Angstrom*10)\n# Local Score: iRMSD, Blue-[%.2f Ang. -- 0.00 Ang.]-Red \n", iRMSD_max);
- }
-
- fprintf ( fp, "\n");
- vfclose (fp);
- free_int (pos, -1);
- return ST;
- }
-float **** analyse_pdb_residues ( Alignment *A, Constraint_list *CL, Pdb_param *pdb_param)
- {
-
- int **pos;
- int s1, s2, rs1, rs2;
- int col1, col2;
- float ****distances;
-
- /*Distances[Nseq][len_aln][4]
- distances...[0]: Number of residues within the bubble
- distances...[1]: Absolute difference of distance of residues within Bubble
- distances...[2]: Number of residues within the bubble with Delta dist < md_threshold
- distances ..[3]: Sum of squared difference of distances
- distances ..[4]: iRMSD
- */
- float d1, d2,delta;
- int wd1, wd2;
- int in_bubble=0;
- int real_res1_col1=0;
- int real_res1_col2;
- int real_res2_col1;
- int real_res2_col2;
- Pdb_param *PP;
- int print_rapdb;
- float nrapdb, rapdb;
- Alignment *BA=NULL;
-
- PP=pdb_param;
- print_rapdb=PP->print_rapdb;
-
- distances=declare_arrayN(4, sizeof (float), A->nseq, A->nseq, 0, 0);
-
- /*Pre-computation of the internal distances----> T[seq]->ca_dist[len][len]*/
- /*Can be avoided if distance_on_request set to 1 */
-
- for ( s1=0; s1< A->nseq; s1++)
- {
- rs1=A->order[s1][0];
- if (CL->T[rs1] && !(CL->T[rs1])->ca_dist)(CL->T[rs1])->ca_dist=measure_ca_distances(CL->T[rs1]);
- for ( s2=0; s2< A->nseq; s2++)
- {
- distances[s1][s2]=declare_float ( A->len_aln, 6);
- }
- }
- pos=aln2pos_simple (A, A->nseq);
-
- for ( s1=0; s1< A->nseq; s1++)
- for ( col1=0; col1< A->len_aln; col1++)
- for ( s2=0; s2<A->nseq; s2++)
- {
- rs1=A->order[s1][0];
- rs2=A->order[s2][0];
- rapdb=0;
- nrapdb=0;
- if ( s1==s2) continue;
- else if (!(CL->T[rs1]) || !(CL->T[rs2]))continue;
- else if ( islower (A->seq_al[s1][col1]) || islower ( A->seq_al[s2][col1]))continue;
- else if ( pos[s1][col1]<=0 || pos[s2][col1]<=0 ) continue;
-
- if ( print_rapdb && s2>s1)
- {
-
- fprintf ( stdout, "RAPDB S1: %s S2: %s POS %d %d %c %d %c ", A->name[s1], A->name[s2], col1+1, pos[s1][col1],A->seq_al[s1][col1], pos[s2][col1],A->seq_al[s2][col1]);
- BA=copy_aln (A, BA);
- lower_string (BA->seq_al[s1]);
- lower_string (BA->seq_al[s2]);
- BA->seq_al[s1][col1]=toupper (BA->seq_al[s1][col1]);
- BA->seq_al[s2][col1]=toupper (BA->seq_al[s2][col1]);
- }
-
- for ( col2=0; col2<A->len_aln; col2++)
- {
-
- if (pos[s1][col2]<=0 || pos[s2][col2]<=0 )continue;
- else if ( FABS((pos[s1][col2]-pos[s1][col1]))<=PP->n_excluded_nb)continue;
- else if ( FABS((pos[s2][col2]-pos[s2][col1]))<=PP->n_excluded_nb)continue;
- else if ( islower (A->seq_al[s1][col2]) || islower ( A->seq_al[s2][col2]))continue;
-
- real_res1_col1=pos[s1][col1]-1;
- real_res1_col2=pos[s1][col2]-1;
-
- real_res2_col1=pos[s2][col1]-1;
- real_res2_col2=pos[s2][col2]-1;
-
- d1=(CL->T[rs1])->ca_dist[real_res1_col1][real_res1_col2];
- d2=(CL->T[rs2])->ca_dist[real_res2_col1][real_res2_col2];
-
- if ( d1==UNDEFINED || d2 == UNDEFINED) continue;
-
-
-
- if ( strm ( PP->local_mode, "sphere"))
- {
- in_bubble= (d1<PP->maximum_distance && d2<PP->maximum_distance)?1:0; ;
- }
- else if ( strm ( PP->local_mode, "window"))
- {
- wd1=FABS((pos[s1][col2]-pos[s1][col1]));
- wd2=FABS((pos[s2][col2]-pos[s2][col1]));
- in_bubble= (wd1<PP->maximum_distance && wd2<PP->maximum_distance)?1:0; ;
- }
-
- if (in_bubble)
- {
- if ( print_rapdb && s2 >s1)
- {
- fprintf ( stdout, "NB %d %d %c %d %c ", col2, pos[s1][col2], A->seq_al[s1][col2], pos[s2][col2], A->seq_al[s2][col2]);
- BA->seq_al[s1][col2]=toupper (BA->seq_al[s1][col2]);
- BA->seq_al[s2][col2]=toupper (BA->seq_al[s2][col2]);
- }
- delta=FABS((d1-d2));
- if (delta<PP->md_threshold)
- distances[s1][s2][real_res1_col1][2]++;
- distances[s1][s2][real_res1_col1][1]+=delta;
- distances[s1][s2][real_res1_col1][0]++;
- distances[s1][s2][real_res1_col1][3]+=delta*delta;
- nrapdb++;
- rapdb+=delta*delta;
- }
- }
-
- if ( nrapdb==0)distances[s1][s2][real_res1_col1][4]=-1;
- else distances[s1][s2][real_res1_col1][4]=(float)sqrt((double)(rapdb/nrapdb));
-
- if ( print_rapdb && s2>s1)
- {
- if (nrapdb==0)
- {
- fprintf ( stdout, "APDB: UNDEFINED\n");
- }
- else
- {
-
- fprintf ( stdout, " APDB: %.2f ",(float)sqrt((double)(rapdb/nrapdb)));
- BA->residue_case=KEEP_CASE;unalign_residues (BA, s1, s2);
- fprintf ( stdout,"SEQ1: %s %s SEQ2: %s %s\n", BA->name[s1], BA->seq_al[s1], BA->name[s2], BA->seq_al[s2]);
- }
- }
-
- }
-
- free_aln (BA);
- free_int (pos, -1);
- return distances;
- }
-int pair_res_suitable4trmsd (int s1,int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
-int aln_column_contains_gap (Alignment *A, int c);
-float aln2ncol4trmsd(Alignment *A, int **pos, Constraint_list *CL, int **lc);
-int pair_columns_suitable4trmsd(int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
-int column_is_suitable4trmsd(int col1,Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s);
-
-
-
-NT_node trmsdmat2tree (float **dm, int **count,Alignment *A);
-Alignment * msa2struc_dist ( Alignment *A, Alignment *ST, char *results, int gapped, int min_ncol4trmsd)
- {
-
- int **pos, c;
- FILE *tl;
- int s1, s2, rs1, rs2;
- int col1, col2;
- float ****distances;
- float **dm,**tdm;
- int **count,**tcount;
- int print_subtrees=0;
- float min, max;
-
- /*Distances[Nseq][len_aln][4]
- distances...[0]: Number of residues within the bubble
- distances...[1]: Absolute difference of distance of residues within Bubble
- distances...[2]: Number of residues within the bubble with Delta dist < md_threshold
- distances ..[3]: Sum of squared difference of distances
- distances ..[4]: iRMSD
- */
- Pdb_param *pdb_param;
- Constraint_list *CL;
- int a, b, ncol, npos,n;
- float d1, d2,delta;
- int wd1, wd2;
- int in_bubble=0;
- int real_res1_col1=0;
- int real_res1_col2;
- int real_res2_col1;
- int real_res2_col2;
- Pdb_param *PP;
- int print_rapdb;
- float nrapdb, rapdb;
- Alignment *BA=NULL;
- NT_node *T0,*T1,*T2,*PT, *POS;
- NT_node BT0, BT10,BT50, BT100=NULL,RBT;
- char **pair_pos_list;
-
- int ntree=0, ntree2;
-
- Alignment *B;
- char *pos_list;
- char *tot_pos_list;
-
- char *struc_tree10;
- char *struc_tree100;
- char *struc_tree50;
- char *struc_tree0;
- char *consense_file;
-
- char *color_struc_tree;
- int **score;
- int proceed=1;
- int **lc;
- int used;
-
- if (min_ncol4trmsd<0)
- {
- min_ncol4trmsd*=-1;
- min_ncol4trmsd=(min_ncol4trmsd*A->len_aln)/100;
- }
- else if ( min_ncol4trmsd>=A->len_aln)
- {
- min_ncol4trmsd=A->len_aln-1;
- }
-
- lc=declare_int (A->nseq, 2);
- for (a=0; a<A->nseq; a++)lc[a][0]=a;
-
- declare_name(tot_pos_list);
- sprintf ( tot_pos_list, "%s.struc_tree.list", results);
-
- declare_name(consense_file);
- sprintf (consense_file, "%s.struc_tree.consense_output", results);
-
- declare_name(pos_list);
- sprintf ( pos_list, "%s.pos_list", results);
-
- declare_name(struc_tree0);
- sprintf ( struc_tree0, "%s.struc_tree.consensus",results);
-
- declare_name(struc_tree10);
- sprintf ( struc_tree10, "%s.struc_tree10",results);
-
- declare_name(struc_tree100);
- sprintf ( struc_tree100, "%s.struc_tree100",results);
-
- declare_name(struc_tree50);
- sprintf ( struc_tree50, "%s.struc_tree50",results);
-
- declare_name(color_struc_tree);
- sprintf ( color_struc_tree, "%s.struc_tree.html", results);
-
- pair_pos_list=declare_char (A->len_aln*A->len_aln+1, 100);
- T1=vcalloc (A->len_aln*A->len_aln+1, sizeof (NT_node));
- T2=vcalloc (A->len_aln+1, sizeof (NT_node));
-
- PT=vcalloc (A->len_aln*A->len_aln+1, sizeof (NT_node));
- POS=vcalloc (A->len_aln+1, sizeof (NT_node));
-
- CL=A->CL;
-
- //Check all sequences have a PDB structure
-
- for (used=0,a=0; a<A->nseq; a++)
- {
- if ( ! seq2P_template_file(A->S,a))
- {
- add_warning (stderr, "Sequence %s removed from the dataset [No Usable Structure]", A->name[a]);
- }
- else
- {
- if (used!=a)
- {
- sprintf (A->name[used], "%s", A->name[a]);
- sprintf (A->seq_al[used], "%s", A->seq_al[a]);
- for (b=0; b<4; b++)A->order[used][b]=A->order[a][b];
- }
- used++;
- }
- }
-
- A->nseq=used;
-
- if (A->nseq<2)myexit (fprintf_error(stderr, "Two sequences at least must have a known structure"));
-
- for ( s1=0; s1< (A->S)->nseq; s1++)
- if ( CL->T[s1]){PP=(CL->T[s1])->pdb_param;break;}
-
- for ( s1=0; s1< A->nseq; s1++)
- {
- rs1=A->order[s1][0];
- if (CL->T[rs1] && !(CL->T[rs1])->ca_dist)(CL->T[rs1])->ca_dist=measure_ca_distances(CL->T[rs1]);
- }
- pos=aln2pos_simple (A, A->nseq);
-
- dm=declare_float (A->nseq, A->nseq);
- count=declare_int (A->nseq, A->nseq);
- tdm=declare_float (A->nseq, A->nseq);
- tcount=declare_int (A->nseq, A->nseq);
-
- PP->maximum_distance=1000;
- sprintf ( PP->local_mode, "sphere");
-
- while ((npos=aln2ncol4trmsd(A,pos,CL,lc))<min_ncol4trmsd && A->nseq>1)
- {
-
- sort_int_inv (lc,2, 1, 0,A->nseq-1);
- add_information (stderr, "Remove Sequence [%s] that contains %d un-suitable positions", A->name[lc[0][0]], lc[0][1]);
- A=remove_seq_from_aln (A, A->name[lc[0][0]]);
- ungap_aln (A);
- pos=aln2pos_simple (A, A->nseq);
- }
- if (!A->nseq)
- {
- myexit (fprintf_error(stderr,"No suitable pair of column supporting a tree"));
- }
- else
- fprintf ( stderr, "\n---- Number of usable positions: %d [%.2f %%]\n", npos, ((float)npos*100)/(float)A->len_aln);
-
- tl=vfopen (tot_pos_list, "w");
- for (ncol=0,ntree=0, col1=0; col1< A->len_aln; col1++)
- {
- int w,tree, cont;
- //output_completion (stderr, col1, A->len_aln,1, "Sample Columns");
- if (!gapped && aln_column_contains_gap (A, col1))continue;
- for ( cont=1,ntree2=0,col2=0; col2<A->len_aln; col2++)
- {
- for (s1=0; s1< A->nseq-1; s1++)
- {
- rs1=A->order[s1][0];
- if (!pair_res_suitable4trmsd (s1,col1, col2, A, pos, PP, CL, &w))continue;
- for ( s2=s1+1; s2<A->nseq; s2++)
- {
- if (!pair_res_suitable4trmsd (s2,col1, col2, A, pos, PP, CL, &w))continue;
-
- rs2=A->order[s2][0];
- real_res1_col1=pos[s1][col1]-1;
- real_res1_col2=pos[s1][col2]-1;
- real_res2_col1=pos[s2][col1]-1;
- real_res2_col2=pos[s2][col2]-1;
-
- d1=(CL->T[rs1])->ca_dist[real_res1_col1][real_res1_col2];
- d2=(CL->T[rs2])->ca_dist[real_res2_col1][real_res2_col2];
-
- delta=FABS((d1-d2));
- dm[s1][s2]=dm[s2][s1]+=delta;
- tdm[s1][s2]=tdm[s2][s1]+=delta;
- tcount[s1][s2]++;
- tcount[s2][s1]++;
-
- count[s1][s2]++;
- count[s2][s1]++;
- }
- }
- }
-
-
-
- if ((POS[col1]=trmsdmat2tree (dm, count, A)))
- {
- T1[ntree]=POS[col1];
- fprintf (tl, "\n>Tree_%d Column\n", col1+1);
- print_tree (T1[ntree], "newick", tl);
- ntree++;
- }
- }
- vfclose (tl);
-
- if (!ntree){fprintf ( stderr, "\nERROR: No suitable pair of column supporting a tree [FATAL]\n"); exit (EXIT_SUCCESS);}
-
- score=treelist2avg_treecmp (T1, NULL);
- display_output_filename( stderr,"TreeList","newick",tot_pos_list, CHECK);
-
- if (treelist_file2consense (tot_pos_list, NULL, consense_file))
- {
- display_output_filename( stderr,"ConsenseTree","phylip",consense_file, CHECK);
- }
- else
- {
- fprintf ( stderr, "\nPhylip is not installed: the program could not produce the consense output. This is not mandatory but useful");
- }
-
- //consensus tree
-
- if ((BT100=treelist2filtered_bootstrap (T1, NULL,score, 1.0)))
- {
- vfclose (print_tree (BT100,"newick", vfopen (struc_tree0, "w")));
- display_output_filename( stderr,"Tree","newick",struc_tree0, CHECK);
- }
- if (print_subtrees)
- {
-
- if ( (BT0=trmsdmat2tree (tdm, tcount, A)))
- {
- vfclose (print_tree (BT0,"newick", vfopen (struc_tree0, "w")));
- display_output_filename( stderr,"Tree","newick",struc_tree0, CHECK);
- }
- if ((BT10=treelist2filtered_bootstrap (T1, NULL,score, 0.1)))
- {
- vfclose (print_tree (BT10,"newick", vfopen (struc_tree10, "w")));
- display_output_filename( stderr,"Tree","newick",struc_tree10, CHECK);
- }
-
- if ((BT50=treelist2filtered_bootstrap (T1, NULL, score,0.5)))
- {
- vfclose (print_tree (BT50,"newick", vfopen (struc_tree50, "w")));
- display_output_filename( stderr,"Tree","newick",struc_tree50, CHECK);
- }
- }
-
-
- if (!BT100)BT100=treelist2filtered_bootstrap (T1, NULL,score, 1.0);
-
- RBT=BT100;
- if (RBT)
- {
- B=copy_aln (A, NULL);
- for (a=0; a<A->len_aln; a++)
- {
- int score;
- Tree_sim *S=NULL;
-
- if (POS[a])
- {
- S=tree_cmp (POS[a], RBT);
- score=S->uw/10;
- }
- else
- {
- score=NO_COLOR_RESIDUE;
- }
-
- for (b=0; b<B->nseq; b++)
- {
- if ( is_gap (B->seq_al[b][a]) || score == NO_COLOR_RESIDUE)
- {
- B->seq_al[b][a]=NO_COLOR_RESIDUE;
- }
- else
- {
- B->seq_al[b][a]=S->uw/10;
- }
- }
- if (S)vfree (S);
- }
-
- output_format_aln ("score_html", A,B,color_struc_tree);
- display_output_filename( stderr,"Colored MSA","score_html",color_struc_tree, CHECK);
- free_aln (BA);
- fprintf ( stderr, "\n");
- }
- fprintf ( stderr, "\n");
- free_int (pos, -1);
- exit (EXIT_SUCCESS);
- return NULL;
- }
-NT_node trmsdmat2tree (float **dm, int **count,Alignment *A)
-{
- float min, max;
- int s1, s2;
- NT_node T;
- int ns;
- int **dm_int;
-
- ns=A->nseq;
- for (s1=0; s1<ns-1; s1++)
- for (s2=s1+1; s2<ns; s2++)
- {
- if ( count [s1][s2])dm[s1][s2]=dm[s2][s1]=dm[s1][s2]/(float)count[s1][s2];
- else
- {
- return NULL;
- }
- if (s1==0 && s2==1)min=max=dm[s1][s2];
- min=MIN(dm[s1][s2], min);
- max=MAX(dm[s1][s2], max);
- }
- dm_int=declare_int (ns, ns);
- for (s1=0; s1<A->nseq-1; s1++)
- for (s2=s1+1; s2<A->nseq; s2++)
- {
- dm_int[s1][s2]=dm_int[s2][s1]=((dm[s1][s2])/(max))*100;
- }
- T=compute_std_tree_2(A, dm_int, "_TMODE_upgma");
- free_int (dm_int, -1);
- for (s1=0; s1<ns; s1++)for ( s2=0; s2<ns; s2++){dm[s1][s2]=count[s1][s2]=0;}
- return T;
-}
-
-int pair_res_suitable4trmsd (int s1,int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
-{
- int rs;
- rs=A->order[s1][0];
- if ( !(CL->T[rs])){s[0]=s1; return 0;}
- else if (is_gap (A->seq_al[s1][col1])){s[0]=s1;return 0;}
- else if (is_gap (A->seq_al[s1][col2])){s[0]=s1;return 0;}
-
- else if (islower(A->seq_al[s1][col1])){s[0]=s1; return 0;}
- else if (islower(A->seq_al[s1][col2])){s[0]=s1; return 0;}
-
- else if ( FABS(((pos[s1][col2])-(pos[s1][col1])))<=PP->n_excluded_nb){s[0]=s1;return 0;}
- else if ((CL->T[rs])->ca_dist[pos[s1][col1]-1][pos[s1][col2]-1]==UNDEFINED){s[0]=s1;return 0;}
- return 1;
-}
-int pair_columns_suitable4trmsd(int col1, int col2, Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
-{
- int s1;
- if (!column_is_suitable4trmsd (col1, A, pos, PP, CL,s))return 0;
- if (!column_is_suitable4trmsd (col2, A, pos, PP, CL,s))return 0;
- for (s1=0; s1<A->nseq; s1++)
- {
- int rs, rr1, rr2;
-
- rs=A->order[s1][0];
- if ( FABS(((pos[s1][col2])-(pos[s1][col1])))<=PP->n_excluded_nb){s[0]=s1;return 0;}
- if ((CL->T[rs])->ca_dist[pos[s1][col1]-1][pos[s1][col2]-1]==UNDEFINED){s[0]=s1;return 0;}
- rr1=pos[s1][col1]-1;
- rr2=pos[s1][col2]-1;
- if ((CL->T[rs])->ca_dist[rr1][rr2]>PP->maximum_distance){s[0]=s1;return 0;}
- }
- return 1;
-}
-int column_is_suitable4trmsd(int col1,Alignment *A, int **pos, Pdb_param *PP, Constraint_list *CL,int *s)
-{
- int s1;
- for ( s1=0; s1<A->nseq; s1++)
- {
- int rs;
- rs=A->order[s1][0];
- if ( !(CL->T[rs])){s[0]=s1; return 0;}
- else if (is_gap (A->seq_al[s1][col1])){s[0]=s1;return 0;}
- else if (islower(A->seq_al[s1][col1])){s[0]=s1; return 0;}
- }
- return 1;
-}
-int aln_column_contains_gap (Alignment *A, int c)
-{
- int a, b;
- if ( !A || c>=A->len_aln || c<0)
- {
- printf ( "\nERROR: values out of range in aln_column_contains_gap [FATL:%s]\n", PROGRAM);
- exit (EXIT_FAILURE);
- }
- for ( a=0; a<A->nseq; a++) if ( is_gap(A->seq_al[a][c]))return 1;
- return 0;
-}
-
-
-float aln2ncol4trmsd(Alignment *A, int **pos, Constraint_list *CL, int **lc)
-{
- //This function estimates the number of columns suitable for constructing a trmsd
- int col1, s1, ncol, n, rs1, real_res1_col1;
-
- for (s1=0; s1<A->nseq; s1++){lc[s1][0]=s1; lc[s1][1]=0;}
- for (ncol=0,col1=0; col1< A->len_aln; col1++)
- {
- for (n=0,s1=0; s1<A->nseq; s1++)
- {
- real_res1_col1=pos[s1][col1]-1;
- rs1=A->order[s1][0];
-
- if (real_res1_col1<0)lc[s1][1]++;
- else if (!((CL->T[A->order[s1][0]])->ca[real_res1_col1]))lc[s1][1]++;
- else n++;
- }
- if (n==A->nseq)
- {
- ncol++;
- }
- }
- return ncol;
-}
-
-float square_atom ( Atom *X)
-{
-
- return X->x*X->x + X->y*X->y + X->z*X->z;
-}
-Atom* reframe_atom ( Atom *X, Atom*Y, Atom *Z, Atom *IN, Atom *R)
- {
- float new_x, new_y, new_z;
-
- if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
-
-
- new_x= X->x*IN->x + Y->x*IN->y +Z->x*IN->z;
- new_y= X->y*IN->x + Y->y*IN->y +Z->y*IN->z;
- new_z= X->z*IN->x + Y->z*IN->y +Z->z*IN->z;
-
- R->x=new_x;
- R->y=new_y;
- R->z=new_z;
- return R;
- }
-
-Atom* add_atom ( Atom *A, Atom*B, Atom *R)
-{
- if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
-
- R->x=A->x+B->x;
- R->y=A->y+B->y;
- R->z=A->z+B->z;
-
- return R;
-}
-Atom* diff_atom ( Atom *A, Atom*B, Atom *R)
-{
- if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
-
- R->x=A->x-B->x;
- R->y=A->y-B->y;
- R->z=A->z-B->z;
-
- return R;
-}
-
-Atom * copy_atom ( Atom *A, Atom*R)
-{
- if ( R==NULL)R=vcalloc ( 1, sizeof (Atom));
- R->num=A->num;
- R->res_num=A->res_num;
- R->x=A->x;
- R->y=A->y;
- R->z=A->z;
-
- sprintf( R->type, "%s", A->type);
- return R;
-}
- void print_atom (Atom *A)
-{
- fprintf ( stdout, "%.2f %.2f %.2f", A->x, A->y, A->z);
-}
-/************************************************************************/
-/* */
-/* NUSSINOV */
-/* */
-/************************************************************************/
-
-/*---------prototypes ----------*/
-static void computeBasePairMatrix(int**M,char*S,int l, int T);
-static int backtrack(int a,int b,int**M,char*S,char*P, int T);
-
-
-
-static int basePair(char x, char y)
-{
- static short **mat;
-
- if (!mat)
- {
- char alp[20];
- int a, b, c1, c2, lc1, lc2;
- mat=declare_short (256, 256);
- sprintf ( alp, "AGCTUagctu");
- for (a=0; a<strlen (alp); a++)
- {
- for (b=a; b<strlen (alp)-1; b++)
- {
- c1=alp[a];c2=alp[b];
- lc1=tolower(c1); lc2=tolower(c2);
- if ( lc1=='g' && lc2=='c')
- mat[c1][c2]=1;
- else if ( lc1=='a' && lc2=='u')
- mat[c1][c2]=1;
- else if ( lc1=='u' && lc2=='g')
- mat [c1][c2]=1;
- mat[c2][c1]=mat[c1][c2];
- }
- }
- }
- return (int)mat[(int)x][(int)y];
-}
-
-
-
-/* ------------------------------------------------------------ */
-
-char *nussinov(char *S, int THRESHOLD)
-{
- char *paren;
- int i;
-
- /*-------------------------------
- S is RNA sequence
- paren is parenthesis expression for
- optimal RNA secondary structure
- THRESHOLD: Min distance between two paired residues
- -------------------------------*/
-
- int **numBasePairs;
- int n;
-
- /*----- initialization --*/
- n = strlen(S);
- paren=vcalloc (n+1, sizeof (char));
- numBasePairs=declare_int (n,n);
-
- for (i=0;i<n;i++) paren[i]='.';
- paren[n]='\0'; // paren is string of same length as S
- computeBasePairMatrix(numBasePairs,S,n, THRESHOLD);
- backtrack(0,n-1,numBasePairs,S,paren, THRESHOLD);
- free_int (numBasePairs, -1);
- return paren;
-}
-
-static void computeBasePairMatrix(int** numBasePairs,char *S,int n, int THRESHOLD)
-{
- int i,j,d,k,max,val,index;
-
- for (d = THRESHOLD; d < n; d++){
- for(i=0; i < n; i++)
- {
- j=i+d;
- if (j < n){
- max=0;
- index=n;
- /*-------------------------------------
- if index<n at end of for-loop, then this
- means that index and j form a base pair,
- and this is noted by numBasePairs[j][i]=index.
- if index=n at end of for-loop, then this
- means that j is not base paired.
- -------------------------------------*/
-
- if ( numBasePairs[i][j-1]>max ){
- max = numBasePairs[i][j-1];
- index = n;
- // j not basepaired with some k such that i<k<j
- }
-
- val = basePair(S[i],S[j]) + numBasePairs[i+1][j-1];
- if ( j-i<= THRESHOLD && val > max ){
- max = val;
- index=i;
- }
- for(k=i; k<=j-THRESHOLD; k++){
- val = basePair(S[k],S[j]) + numBasePairs[i][k-1]
- + numBasePairs[k+1][j-1];
- if (val > max) {
- max = val;
- index=k;
- }
- }
- numBasePairs[i][j]=max;
- if (index<n)
- numBasePairs[j][i]=index;
- else
- numBasePairs[j][i]=-1;
- }
- }
- }
-
-}
-
-
-
-
-static int backtrack(int i, int j, int **numBasePairs,char *S, char *paren, int THRESHOLD)
-{
- int k;
-
- k = numBasePairs[j][i];
- if(k != -1)
- {
- paren[k] = '(';
- paren[j] = ')';
- if( THRESHOLD <= (j-1)-(k+1) )
- backtrack(k+1,j-1,numBasePairs,S,paren, THRESHOLD);
- if (THRESHOLD <= k-1-i )
- backtrack(i,k-1,numBasePairs,S,paren, THRESHOLD);
- }
- else{
- if( THRESHOLD <= j-1-i )
- {
- backtrack(i,j-1,numBasePairs,S,paren, THRESHOLD);
- }
- else
- return 0;
- }
- return 0;}
-
-int count;
-char * rna_struc2rna_lib ( char *seq_name, char *seq, char *name)
-{
- FILE *fp;
- char *st;
-
-
- st=nussinov (seq, 2);
- if ( name==NULL)name=vtmpnam(NULL);
- fp=vfopen ( name, "w");
- fprintf (fp, "! TC_LIB_FORMAT_01\n");
- fprintf (fp, "1\n%s %d %s\n", seq_name, (int)strlen (seq), seq);
- fprintf (fp, "#1 1\n");
- display_rna_ss (0, seq, st, fp);
- fprintf ( fp, "! SEQ_1_TO_N\n");
- vfclose (fp);
- vfree (st);
- //printf_system ( "cp %s test", name);
- return name;
-}
-int display_rna_ss ( int n, char *seq, char *st, FILE *fp)
-{
- char p;
- char string[100];
- static int thread;
-
- while ((p=st[n])!='\0')
- {
- if ( p=='(')
- {
- thread=count++;
- sprintf (string, "%d",n+1);
- n=display_rna_ss (n+1, seq, st, fp);
- fprintf (fp, "%s %d 100\n", string, n+1);
- }
- else if (p=='.');
- else if (p==')')
- {
- return n;
- }
- n++;
- }
- return n;
-}
-/******************************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