binaries/src/ViennaRNA/Cluster/AS_main.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <unistd.h>
   4 #include <string.h>
   5 #include "distance_matrix.h"
   6 #include "statgeom.h"
   7 #include "split.h"
   8 #include "cluster.h"
   9 #include "treeplot.h"
  10 #include "utils.h"
  11
  12
  13 #define PUBLIC
  14 #define PRIVATE   static
  15
  16 PRIVATE char  scale1[] = "....,....1....,....2....,....3....,....4";
  17 PRIVATE char  scale2[] = "....,....5....,....6....,....7....,....8";
  18
  19 PRIVATE void usage(void);
  20
  21 main(int argc, char *argv[])
  22 {
  23    int      n,i,j,l;
  24    int      intty;
  25    int      outtty;
  26    char    *mask, junk[20];
  27    char   **s;
  28    char   **ss[4];
  29    float   *B;
  30    float  **dm;
  31    Split   *S;
  32    Union   *U;
  33    char     DistAlgorithm='H';
  34    int      nn[4];
  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;
  37
  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] ) {
  43           case 'X':
  44              if (argv[i][2]=='\0') { Do_Stg = 1 ; break; }
  45              Do_Split = 0;
  46              Do_Wards = 0;
  47              Do_Stg   = 0;
  48              Do_Nj    = 0;
  49              for(j=2;j<strlen(argv[i]);j++) {
  50                 switch(argv[i][j]) {
  51                  case 's' :  Do_Split = 1;
  52                     break;
  53                  case 'w' :  Do_Wards = 1;
  54                     break;
  55                  case 'b' :  Do_Stg   = 1;
  56                     break;
  57                  case 'n' :  Do_Nj    = 1;
  58                     break;
  59                  case 'm' :  Do_Mat   = 1;
  60                    break;
  61                  default :
  62                     usage();
  63                 }
  64              }
  65              break;
  66           case 'Q':
  67              Do_4_Stg = 1;
  68              break;
  69           case 'M':
  70              if(mask) { free(mask); mask = NULL; }
  71              switch (argv[i][2] ) {
  72               case '\0' :
  73                 usage();
  74                 break;
  75               case 'a' :
  76                 mask   = space(sizeof(char)*54);
  77                 strcpy(mask,
  78                 "%ABCDEFGHIJKLMNOPQRSTUVWabcdefghijklmnopqrstuvwxyz");
  79                 if(argv[i][3]=='+') mask[0] = '~';  /* make case sensitive */
  80                 break;
  81               case 'u' :
  82                 mask   = space(sizeof(char)*28);
  83                 strcpy(mask,"~ABCDEFGHIJKLMNOPQRSTUVW");
  84                 break;
  85               case 'l' :
  86                 mask   = space(sizeof(char)*28);
  87                 strcpy(mask,"~abcdefghijklmnopqrstuvwxyz");
  88                 break;
  89               case 'c' :
  90                 mask   = space(sizeof(char)*12);
  91                 strcpy(mask,"~1234567890");
  92                 break;
  93               case 'n' :
  94                 mask   = space(sizeof(char)*64);
  95                 strcpy (mask,
  96                 "%ABCDEFGHIJKLMNOPQRSTUVWabcdefghijklmnopqrstuvwxyz1234567890");
  97                 if(argv[i][3]=='+') mask[0] = '~';  /* make case sensitive */
  98                 break;
  99               case 'R' :     /* RNA */
 100                 mask   = space(sizeof(char)*10);
 101                 strcpy(mask,"%GCAUgcau");
 102                 if(argv[i][3]=='+') mask[0] = '~';  /* make case sensitive */
 103                 break;
 104               case 'D' :     /* DNA */
 105                 mask   = space(sizeof(char)*10);
 106                 strcpy (mask,"%GCATgcat");
 107                 if(argv[i][3]=='+') mask[0] = '~';  /* make case sensitive */
 108                 break;
 109               case 'A' :    /*  AMINOACIDS */
 110                 mask   = space(sizeof(char)*42);
 111                 strcpy(mask,"%ACDEFGHIKLMNPQRSTVWYacdefghiklmnpqrstvwy");
 112                 if(argv[i][3]=='+') mask[0] = '~';  /* make case sensitive */
 113                 break;
 114               case 'S' :    /* SECONDARY STRUCTURES */
 115                 mask   = space(sizeof(char)*6);
 116                 strcpy(mask,"~().^");
 117                 break;
 118               case '%' :    /* ARBITRARY ALPHABETS */
 119                 l = strlen(argv[i]);
 120                 if(argv[i][l] == '+'){
 121                    mask =   space(sizeof(char)*(l-2));
 122                    mask[0] = '~';
 123                    for(j=1;j<=l-4;j++) mask[j] = argv[i][j+2];
 124                    mask[l-3]='\0';
 125                 }
 126                 if(argv[i][l] == '!'){
 127                    mask =   space(sizeof(char)*(l-2));
 128                    mask[0] = '!';
 129                    for(j=1;j<=l-4;j++) mask[j] = argv[i][j+2];
 130                    mask[l-3]='\0';
 131                 }
 132                 else {
 133                    mask =   space(sizeof(char)*(l-1));
 134                    mask[0] = '%';
 135                    for(j=1;j<=l-3;j++) mask[j] = argv[i][j+2];
 136                    mask[l-2]='\0';
 137                 }
 138                 break;
 139               default :
 140                 usage();
 141             }
 142             break;
 143           case 'D' :               /* choose algorithm */
 144             switch(argv[i][2]) {
 145               case '\0' :
 146                 usage();
 147                 break;
 148               case 'H' :
 149                 DistAlgorithm = 'H';
 150                 break;
 151               case 'A' :
 152                 DistAlgorithm = 'A';
 153                 if(argv[i][3]==',') {
 154                    per_digit=-1.;
 155                    sscanf(argv[i],"%[^,],%g",junk,&per_digit);
 156                    if(per_digit<0.) usage();
 157                    per_gap = per_digit;
 158                    Set_StrEdit_GapCosts(per_digit,per_gap);
 159                 }
 160                 break;
 161               case 'G' :
 162                 DistAlgorithm = 'G';
 163                 if(argv[i][3]==',') {
 164                    per_digit=-1.;
 165                    per_gap  =-1.;
 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);
 169                 }
 170                 break;
 171               default:
 172                 usage();
 173              }
 174              break;
 175            case 'd' :               /* choose distance matrix */
 176              switch(argv[i][2]) {
 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]);
 183                  break;
 184                default :
 185                  usage();
 186              }
 187              break;
 188           default :
 189              usage();
 190          }
 191       }
 192    }
 193
 194    /* END PARSING OF COMMAND LINE */
 195
 196    intty = isatty(fileno(stdin));
 197    outtty= isatty(fileno(stdout));
 198
 199    if(intty){
 200       if(outtty) {
 201          printf("Input sequences; @ to mark end of input\n");
 202          printf("%s%s\n", scale1, scale2);
 203       }
 204       else {
 205          fprintf(stderr,"Input sequences; @ to mark end of input\n");
 206          fprintf(stderr,"%s%s\n", scale1, scale2);
 207       }
 208    }
 209
 210    while ((s=read_sequence_list(&n,mask))!=NULL) {
 211       dm = NULL;
 212       if(Do_4_Stg) {
 213          ss[0] = s;
 214          nn[0] = n;
 215          for(i=1;i<4;i++) {
 216             ss[i] = read_sequence_list(&n,mask);
 217             if(ss[i]==NULL) nrerror("read_sequences: wrong or insufficient data.");
 218             nn[i] = n;
 219          }
 220          printf_taxa_list();
 221          B = statgeom4(ss,nn);
 222          printf_stg(B);
 223          SimplifiedBox(B,"box.ps");    /* This is preliminary !!! */
 224          free(B);
 225
 226          for(i=0;i<4;i++){
 227             for(j=0;j<nn[i];j++) free(ss[i][j]);
 228             free(ss[i]);
 229          }
 230          /* free(ss); */ /* attempt to free a non-heap object */
 231       }
 232       else {
 233          printf_taxa_list();
 234          if(Do_Stg) {
 235             B = statgeom(s,n);
 236             if (B) {
 237                printf_stg(B);
 238                SimplifiedBox(B,"box.ps");
 239                free(B);
 240             }
 241          }
 242          if((Do_Split)||(Do_Wards)||(Do_Nj)||(Do_Mat)) {
 243             switch(DistAlgorithm) {
 244              case 'H' :
 245                 dm = Hamming_Distance_Matrix(s,n);
 246                 printf("> %s\n","H (Hamming Distance)");
 247                 break;
 248               case 'A' :
 249                 dm = StrEdit_SimpleDistMatrix(s,n);
 250                 printf("> %s\n","A (Needleman-Wunsch Distance)");
 251                 break;
 252               case 'G' :
 253                  dm = StrEdit_GotohDistMatrix(s,n);
 254                 printf("> %s\n","G (Gotoh Distance)");
 255                 break;
 256               default:
 257                 nrerror("This can't happen.");
 258              }
 259          }
 260          if(Do_Split) {
 261             S = split_decomposition(dm);
 262             sort_Split(S);
 263             print_Split(S);
 264             free_Split(S);
 265          }
 266          if(Do_Wards) {
 267             U = wards_cluster(dm);
 268             printf_phylogeny(U,"W");
 269             PSplot_phylogeny(U,"wards.ps","Ward's Method");
 270             free(U);
 271          }
 272          if(Do_Nj) {
 273             U = neighbour_joining(dm);
 274             printf_phylogeny(U,"Nj");
 275             PSplot_phylogeny(U,"nj.ps","Neighbor Joining");
 276             free(U);
 277          }
 278          if(Do_Mat) printf_distance_matrix(dm);
 279       }
 280
 281       if (dm!=NULL) free_distance_matrix(dm);
 282       for(i=0;i<n;i++) free(s[i]);
 283       free(s);
 284    }
 285    return 0;
 286 }
 287
 288 /* -----------------------------------------------------------------------*/
 289
 290 PRIVATE void usage(void)
 291 {
 292    nrerror("usage: AnalyseSeqs [-X[bswnm]] [-Q] [-M{mask}] \n"
 293    "                   [-D{H|A[,cost]|G[,cost1,cost2]}] [-d{D|B|H|S}]");
 294 }