5 #include "distance_matrix.h"
14 #define PRIVATE static
16 PRIVATE char scale1[] = "....,....1....,....2....,....3....,....4";
17 PRIVATE char scale2[] = "....,....5....,....6....,....7....,....8";
19 PRIVATE void usage(void);
21 main(int argc, char *argv[])
33 char DistAlgorithm='H';
35 short Do_Split=0, Do_Wards=0, Do_Stg=1, Do_4_Stg=0, Do_Nj=0, Do_Mat=0;
36 float per_digit, per_gap;
38 mask = space(sizeof(char)*54);
39 strcpy (mask,"%ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz");
40 for (i=1; i<argc; i++) {
41 if (argv[i][0]=='-') {
42 switch ( argv[i][1] ) {
44 if (argv[i][2]=='\0') { Do_Stg = 1 ; break; }
49 for(j=2;j<strlen(argv[i]);j++) {
51 case 's' : Do_Split = 1;
53 case 'w' : Do_Wards = 1;
55 case 'b' : Do_Stg = 1;
59 case 'm' : Do_Mat = 1;
70 if(mask) { free(mask); mask = NULL; }
71 switch (argv[i][2] ) {
76 mask = space(sizeof(char)*54);
78 "%ABCDEFGHIJKLMNOPQRSTUVWabcdefghijklmnopqrstuvwxyz");
79 if(argv[i][3]=='+') mask[0] = '~'; /* make case sensitive */
82 mask = space(sizeof(char)*28);
83 strcpy(mask,"~ABCDEFGHIJKLMNOPQRSTUVW");
86 mask = space(sizeof(char)*28);
87 strcpy(mask,"~abcdefghijklmnopqrstuvwxyz");
90 mask = space(sizeof(char)*12);
91 strcpy(mask,"~1234567890");
94 mask = space(sizeof(char)*64);
96 "%ABCDEFGHIJKLMNOPQRSTUVWabcdefghijklmnopqrstuvwxyz1234567890");
97 if(argv[i][3]=='+') mask[0] = '~'; /* make case sensitive */
100 mask = space(sizeof(char)*10);
101 strcpy(mask,"%GCAUgcau");
102 if(argv[i][3]=='+') mask[0] = '~'; /* make case sensitive */
105 mask = space(sizeof(char)*10);
106 strcpy (mask,"%GCATgcat");
107 if(argv[i][3]=='+') mask[0] = '~'; /* make case sensitive */
109 case 'A' : /* AMINOACIDS */
110 mask = space(sizeof(char)*42);
111 strcpy(mask,"%ACDEFGHIKLMNPQRSTVWYacdefghiklmnpqrstvwy");
112 if(argv[i][3]=='+') mask[0] = '~'; /* make case sensitive */
114 case 'S' : /* SECONDARY STRUCTURES */
115 mask = space(sizeof(char)*6);
116 strcpy(mask,"~().^");
118 case '%' : /* ARBITRARY ALPHABETS */
120 if(argv[i][l] == '+'){
121 mask = space(sizeof(char)*(l-2));
123 for(j=1;j<=l-4;j++) mask[j] = argv[i][j+2];
126 if(argv[i][l] == '!'){
127 mask = space(sizeof(char)*(l-2));
129 for(j=1;j<=l-4;j++) mask[j] = argv[i][j+2];
133 mask = space(sizeof(char)*(l-1));
135 for(j=1;j<=l-3;j++) mask[j] = argv[i][j+2];
143 case 'D' : /* choose algorithm */
153 if(argv[i][3]==',') {
155 sscanf(argv[i],"%[^,],%g",junk,&per_digit);
156 if(per_digit<0.) usage();
158 Set_StrEdit_GapCosts(per_digit,per_gap);
163 if(argv[i][3]==',') {
166 sscanf(argv[i]+4,"%f,%f",&per_digit,&per_gap);
167 if((per_digit<0.)||(per_gap<0)) usage();
168 Set_StrEdit_GapCosts(per_digit,per_gap);
175 case 'd' : /* choose distance matrix */
177 case 'D' : /* Dayhoff Distances */
178 case 'A' : /* Aminoacid Distance (Hofacker & Borstnik) */
179 case 'B' : /* RY distances for nucleotides */
180 case 'H' : /* Hogeweg's Distance for Secondary Structures */
181 case 'S' : /* Simple Distance (superfluous option) */
182 Set_StrEdit_CostMatrix(argv[i][2]);
194 /* END PARSING OF COMMAND LINE */
196 intty = isatty(fileno(stdin));
197 outtty= isatty(fileno(stdout));
201 printf("Input sequences; @ to mark end of input\n");
202 printf("%s%s\n", scale1, scale2);
205 fprintf(stderr,"Input sequences; @ to mark end of input\n");
206 fprintf(stderr,"%s%s\n", scale1, scale2);
210 while ((s=read_sequence_list(&n,mask))!=NULL) {
216 ss[i] = read_sequence_list(&n,mask);
217 if(ss[i]==NULL) nrerror("read_sequences: wrong or insufficient data.");
221 B = statgeom4(ss,nn);
223 SimplifiedBox(B,"box.ps"); /* This is preliminary !!! */
227 for(j=0;j<nn[i];j++) free(ss[i][j]);
230 /* free(ss); */ /* attempt to free a non-heap object */
238 SimplifiedBox(B,"box.ps");
242 if((Do_Split)||(Do_Wards)||(Do_Nj)||(Do_Mat)) {
243 switch(DistAlgorithm) {
245 dm = Hamming_Distance_Matrix(s,n);
246 printf("> %s\n","H (Hamming Distance)");
249 dm = StrEdit_SimpleDistMatrix(s,n);
250 printf("> %s\n","A (Needleman-Wunsch Distance)");
253 dm = StrEdit_GotohDistMatrix(s,n);
254 printf("> %s\n","G (Gotoh Distance)");
257 nrerror("This can't happen.");
261 S = split_decomposition(dm);
267 U = wards_cluster(dm);
268 printf_phylogeny(U,"W");
269 PSplot_phylogeny(U,"wards.ps","Ward's Method");
273 U = neighbour_joining(dm);
274 printf_phylogeny(U,"Nj");
275 PSplot_phylogeny(U,"nj.ps","Neighbor Joining");
278 if(Do_Mat) printf_distance_matrix(dm);
281 if (dm!=NULL) free_distance_matrix(dm);
282 for(i=0;i<n;i++) free(s[i]);
288 /* -----------------------------------------------------------------------*/
290 PRIVATE void usage(void)
292 nrerror("usage: AnalyseSeqs [-X[bswnm]] [-Q] [-M{mask}] \n"
293 " [-D{H|A[,cost]|G[,cost1,cost2]}] [-d{D|B|H|S}]");