binaries/src/ViennaRNA/Cluster/distance_matrix.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <ctype.h>
   5 #include "utils.h"
   6 #include "StrEdit_CostMatrix.h"
   7
   8 #define  PUBLIC
   9 #define  PRIVATE         static
  10 #define  MAXSEQS         1000
  11
  12 PUBLIC   float **read_distance_matrix(char type[]);
  13 PUBLIC   char  **read_sequence_list(int *n_of_seqs, char *mask);
  14 PUBLIC   float **Hamming_Distance_Matrix(char **seqs, int n_of_seqs);
  15 PUBLIC   float **StrEdit_SimpleDistMatrix(char **seqs, int n_of_seqs);
  16 PUBLIC   float **StrEdit_GotohDistMatrix(char **seqs, int n_of_seqs);
  17 PUBLIC   void    free_distance_matrix(float **x);
  18 PUBLIC   void    printf_distance_matrix(float **x);
  19 PUBLIC   void    printf_taxa_list(void);
  20 PUBLIC   char   *get_taxon_label(int whoami);
  21 PUBLIC   float   StrEdit_SimpleDist(char *str1, char *str2);
  22 PUBLIC   float   StrEdit_GotohDist(char *str1, char *str2);
  23 PUBLIC   void    Set_StrEdit_CostMatrix(char type);
  24 PUBLIC   void    Set_StrEdit_GapCosts(float per_digit, float per_gap);
  25
  26 /* NOTE:   x[0][0] = (float)size_of_matrix;    */
  27
  28 PRIVATE  void    read_taxa_list(void);
  29 PRIVATE  int     string_consists_of(char line[],char *mask);
  30 PRIVATE  float   StrEditCost( int i, int j, char *T1, char *T2);
  31 PRIVATE  int     decode(char id);
  32
  33 PRIVATE  char    Taxa_List[MAXSEQS][50];
  34 PRIVATE  int     Taxa_Numbers[MAXSEQS];
  35 PRIVATE  int     N_of_named_taxa=0;
  36 PRIVATE  char   *file_name;
  37 PRIVATE  char    N_of_infiles=0;
  38 PRIVATE  float **StrEdit_CostMatrix;
  39 PRIVATE  char   *StrEdit_ValidAlphabet;
  40 PRIVATE  float   StrEdit_GapCost   = 1.;
  41 PRIVATE  float   StrEdit_GotohAlpha = 1.;
  42 PRIVATE  float   StrEdit_GotohBeta  = 1.;
  43
  44
  45
  46 PUBLIC float **read_distance_matrix(char type[])
  47 {
  48    char   *line;
  49    float **D;
  50    float   tmp;
  51    int     i,j,size;
  52
  53    while(1) {
  54      type[0]= '\0';
  55      size   =    0;
  56      D      = NULL;
  57      if ((line = get_line(stdin))==NULL) return NULL;
  58      if (*line =='@') return NULL;
  59      if (*line =='*') {
  60        N_of_infiles++;
  61        if(file_name) free(file_name);
  62
  63        if(strlen(line)>1) {
  64          file_name = (char *) space(sizeof(char)*strlen(line));
  65          sscanf(line,"*%s",file_name);
  66        } else {
  67          file_name = (char *) space(10);
  68          sprintf(file_name,"%d",N_of_infiles);
  69        }
  70        read_taxa_list();
  71      }
  72      else if (*line=='>') {
  73        int r;
  74        size = 0;
  75        r = sscanf(line,"> %1s%*[ ] %d", type, &size);
  76        fprintf(stderr, "%d ", r);
  77        if (r==EOF) return NULL;
  78        if((r==2)&&(size>1)) {
  79          D=(float **)space((size+1)*sizeof(float *));
  80          for(i=0; i<=size; i++)
  81            D[i] = (float *)space((size+1)*sizeof(float));
  82          D[0][0] = (float)size;
  83          D[1][1] = 0.0;
  84          for(i=2; i<= size; i++) {
  85            D[i][i] = 0.0;
  86            for(j=1; j<i; j++) {
  87              if (scanf("%f", &tmp)!=1) {
  88                for(i=0;i<=size;i++) free(D[i]);
  89                free(D);
  90                return NULL;
  91              }
  92              D[i][j] = tmp;
  93              D[j][i] = tmp;
  94            }
  95          }
  96          return D;
  97        }
  98        else printf("%s\n",line);
  99      }
 100      else printf(" %s\n", line);
 101      free(line);
 102    }
 103 }
 104
 105 /* ------------------------------------------------------------------------- */
 106
 107 PUBLIC char **read_sequence_list(int *n_of_seqs, char *mask)
 108 {
 109    int     i;
 110    char   *line;
 111    char   *tt[MAXSEQS];
 112    char  **sl;
 113    int     len;
 114
 115    (*n_of_seqs) = 0;
 116    while(1) {
 117       if ((line = get_line(stdin))==NULL) break;
 118
 119       if(line[0]=='\0') break;
 120       if(line[0]=='%')  break;
 121       if(line[0]=='#')  break;
 122       if(line[0]=='@')  break;
 123       if(line[0]=='*') {
 124          N_of_infiles++;
 125          if(file_name) free(file_name);
 126          if(strlen(line)>1) {
 127             file_name = (char *) space(sizeof(char)*strlen(line));
 128             sscanf(line,"*%s",file_name);
 129          } else {
 130             file_name = (char *) space(10);
 131             sprintf(file_name,"%d",N_of_infiles);
 132          }
 133          read_taxa_list();
 134          free(line);
 135          continue;
 136       }
 137
 138       len = strlen(line);
 139       if(string_consists_of(line,mask)){
 140          if(mask[0]=='%') {
 141             for(i=0;i<len;i++){
 142                if(isalpha(line[i])) line[i]=toupper(line[i]);
 143             }
 144             if(mask[0]=='!') {
 145                for(i=0;i<=len;i++){
 146                   switch(toupper(line[i])){
 147                    case 'G' :  case 'A' : case 'X' :
 148                       line[i] = 'R'; break;
 149                     case 'U' : case 'C' : case 'K' : case 'T' :
 150                        line[i] = 'Y'; break;
 151                     default:
 152                       line[i] = '*';
 153                    }
 154                }
 155             }
 156             tt[*n_of_seqs] = (char *)space((len+1)*sizeof(char));
 157             sscanf(line,"%s",tt[*n_of_seqs]);
 158             (*n_of_seqs)++;
 159          }
 160       }
 161       free(line);
 162    }
 163    if(*n_of_seqs == 0) return NULL;
 164    else {
 165      sl = (char **) space((*n_of_seqs)*sizeof(char *));
 166      for(i=0;i<*n_of_seqs; i++) sl[i] = tt[i];
 167      return sl;
 168    }
 169 }
 170
 171 /* -------------------------------------------------------------------------- */
 172
 173 PRIVATE int string_consists_of(char *line, char *mask)
 174 {
 175    int i,j,len;
 176    if((len=strlen(line))==0) return 0; /* empty string ! */
 177    for(j=0;j<len;j++){
 178       for(i=1;i<strlen(mask);i++){     /* mask[0] marks toupper/not_toupper ! */
 179          if(mask[i]==line[j]) break;
 180       }
 181       if(i==(strlen(mask))) break;      /* letter j no in mask */
 182    }
 183    if(j!=len) return 0;       /* loop left before last char -> false */
 184    return 1;                           /* loop left after last char -> true */
 185 }
 186 /* -------------------------------------------------------------------------- */
 187
 188 PUBLIC void free_distance_matrix(float **x)
 189 {
 190    int i,n;
 191    n=(int) x[0][0];
 192    for(i=0;i<=n;i++) free(x[i]);
 193    free(x);
 194    x=NULL;
 195 }
 196
 197 /* -------------------------------------------------------------------------- */
 198
 199 PUBLIC void printf_distance_matrix(float **x)
 200 {
 201    int i,j,n;
 202    n=(int) x[0][0];
 203    printf("> X  %d\n",n);
 204    if(n>1){
 205       for(i=2;i<=n;i++) {
 206          for(j=1;j<i;j++) printf("%g ",x[i][j]);
 207          printf("\n");
 208       }
 209    }
 210 }
 211
 212 /* -------------------------------------------------------------------------- */
 213
 214 PUBLIC float **Hamming_Distance_Matrix(char **seqs, int n_of_seqs)
 215 {
 216    int i,j,k;
 217    float **D;
 218    D = (float **) space((n_of_seqs+1)*sizeof(float *));
 219    for(i=0;i<=n_of_seqs;i++)
 220       D[i] = (float *) space((n_of_seqs+1)*sizeof(float));
 221    D[0][0] = (float) n_of_seqs;
 222
 223    for(i=1; i<n_of_seqs; i++) {
 224       D[i][i] = 0.;
 225       for(j=0;j<i;j++){
 226          if(strlen(seqs[i])!=strlen(seqs[j]))
 227             nrerror("Unequal Seqence Length for Hamming Distance.");
 228          D[i+1][j+1] = 0.0;
 229          for(k=0;k<strlen(seqs[i]);k++)
 230             D[i+1][j+1] += StrEditCost(k+1,k+1,seqs[i],seqs[j]);
 231          /* was :  (float)(seqs[i][k]!=seqs[j][k]); */
 232          D[j+1][i+1] = D[i+1][j+1];
 233       }
 234       D[n_of_seqs][n_of_seqs] = 0.;
 235    }
 236    return D;
 237 }
 238
 239 /* -------------------------------------------------------------------------- */
 240
 241 PUBLIC float **StrEdit_SimpleDistMatrix(char **seqs, int n_of_seqs)
 242 {
 243    int i,j;
 244    float **D;
 245    D = (float **) space((n_of_seqs+1)*sizeof(float *));
 246    for(i=0;i<=n_of_seqs;i++)
 247       D[i] = (float *) space((n_of_seqs+1)*sizeof(float));
 248    D[0][0] = (float) n_of_seqs;
 249
 250    for(i=1; i<n_of_seqs; i++) {
 251       D[i][i] = 0.;
 252       for(j=0;j<i;j++){
 253          D[i+1][j+1] = StrEdit_SimpleDist(seqs[i],seqs[j]);
 254          D[j+1][i+1] = D[i+1][j+1];
 255       }
 256       D[n_of_seqs][n_of_seqs] = 0.;
 257    }
 258    return D;
 259 }
 260
 261 /* -------------------------------------------------------------------------- */
 262
 263 PUBLIC float **StrEdit_GotohDistMatrix(char **seqs, int n_of_seqs)
 264 {
 265    int i,j;
 266    float **D;
 267    D = (float **) space((n_of_seqs+1)*sizeof(float *));
 268    for(i=0;i<=n_of_seqs;i++)
 269       D[i] = (float *) space((n_of_seqs+1)*sizeof(float));
 270    D[0][0] = (float) n_of_seqs;
 271
 272    for(i=1; i<n_of_seqs; i++) {
 273       D[i][i] = 0.;
 274       for(j=0;j<i;j++){
 275          D[i+1][j+1] = StrEdit_GotohDist(seqs[i],seqs[j]);
 276          D[j+1][i+1] = D[i+1][j+1];
 277       }
 278       D[n_of_seqs][n_of_seqs] = 0.;
 279    }
 280    return D;
 281 }
 282
 283 /* -------------------------------------------------------------------------- */
 284
 285 PRIVATE void read_taxa_list(void)
 286 {
 287    char *line;
 288    int i,add_it;
 289
 290    add_it=0;
 291    if ((line = get_line(stdin))==NULL) return;
 292    if(line[0]=='#') {
 293       i=0;
 294       sscanf(line,"#%d", &i);
 295       if(i<=1) N_of_named_taxa = 0;
 296       add_it = i*100000;
 297       free(line);
 298       if ((line = get_line(stdin))==NULL) return;
 299    } else N_of_named_taxa = 0;
 300
 301    do {
 302       if(line[0]=='\0') break;
 303       if(line[0]=='%')  break;
 304       if(line[0]=='#')  break;
 305       if(line[0]=='@')  break;
 306       if(line[0]=='*')  break;
 307       *Taxa_List[N_of_named_taxa]='\0';
 308       sscanf(line,"%d :%49s", &i, Taxa_List[N_of_named_taxa]);
 309       if(*Taxa_List[N_of_named_taxa]) {
 310          Taxa_Numbers[N_of_named_taxa]=i+add_it;
 311          N_of_named_taxa++;
 312       }
 313       free(line);
 314    } while ((line = get_line(stdin))!=NULL);
 315    if (line!=NULL) free(line);
 316    return;
 317 }
 318
 319 /* -------------------------------------------------------------------------- */
 320
 321 PUBLIC void printf_taxa_list(void)
 322 {
 323    int i;
 324    if(N_of_named_taxa>0){
 325       printf("* List of Taxa: %s\n", file_name);
 326       for(i=0;i<N_of_named_taxa;i++)
 327          printf("%3d : %s\n",Taxa_Numbers[i],Taxa_List[i]);
 328       printf("* End of Taxa List\n");
 329    }
 330 }
 331
 332 /* -------------------------------------------------------------------------- */
 333
 334 PUBLIC char *get_taxon_label(int whoami)
 335 {
 336    char *label;
 337    char tmp[20];
 338    int i;
 339
 340    if(whoami<0) {    /* negative arguments return the identifier of the data set */
 341       if(!file_name) return NULL;
 342       label = (char *) space(sizeof(char)*(strlen(file_name)+1));
 343       strcpy(label,file_name);
 344       return label;
 345    }
 346    for(i=0;i<N_of_named_taxa;i++) {
 347       if(whoami==Taxa_Numbers[i]) {
 348          label = (char *) space(sizeof(char)*(strlen(Taxa_List[i])+1));
 349          strcpy(label,Taxa_List[i]);
 350          return label;
 351       }
 352    }
 353    sprintf(tmp,"%d",whoami);
 354
 355    label = (char *) space(sizeof(char)*(strlen(tmp)+1));
 356    strcpy(label,tmp);
 357    return label;
 358 }
 359
 360 /* -------------------------------------------------------------------------- */
 361
 362 PUBLIC float StrEdit_SimpleDist(char *str1, char *str2 )
 363
 364 {
 365    float  **distance;
 366
 367    int           i, j, length1,length2;
 368    float         minus, plus, change, temp;
 369
 370    length1 = strlen(str1);
 371    length2 = strlen(str2);
 372
 373    distance = (float **)  space((length1 +1)*sizeof(float *));
 374    for(i=0; i<= length1; i++)
 375       distance[i] = (float *) space( (length2+1)*sizeof(float));
 376
 377    for(i = 1; i <= length1; i++)
 378       distance[i][0] = distance[i-1][0]+StrEditCost(i,0,str1,str2);
 379    for(j = 1; j <= length2; j++)
 380       distance[0][j] = distance[0][j-1]+StrEditCost(0,j,str1,str2);
 381
 382    for (i = 1; i <= length1; i++) {
 383       for (j = 1; j <= length2 ; j++) {
 384          minus  = distance[i-1][j]  + StrEditCost(i,0,str1,str2);
 385          plus   = distance[i][j-1]  + StrEditCost(0,j,str1,str2);
 386          change = distance[i-1][j-1]+ StrEditCost(i,j,str1,str2);
 387
 388          distance[i][j] = MIN3(minus, plus, change);
 389       }
 390    }
 391    temp = distance[length1][length2];
 392    for(i=0;i<=length1;i++) free(distance[i]);
 393    free(distance);
 394
 395    return temp;
 396 }
 397
 398 /* -------------------------------------------------------------------------- */
 399
 400 PUBLIC float StrEdit_GotohDist(char *str1, char *str2 )
 401 {
 402    float  **D;
 403    float  **E;
 404    float  **F;
 405    int      i, j, length1,length2;
 406    float    temp;
 407
 408    length1 = strlen(str1);
 409    length2 = strlen(str2);
 410
 411    D = space((length1+1)*sizeof(float *));
 412    for(i=0;i<=length1;i++) D[i] = space((length2+1)*sizeof(float));
 413    E = space((length1+1)*sizeof(float *));
 414    for(i=0;i<=length1;i++) E[i] = space((length2+1)*sizeof(float));
 415    F = space((length1+1)*sizeof(float *));
 416    for(i=0;i<=length1;i++) F[i] = space((length2+1)*sizeof(float));
 417
 418    D[0][0] = 0.; E[0][0] = 0.; F[0][0] = 0.;
 419
 420    for(i=1;i<=length1;i++) {
 421       D[i][0] = StrEdit_GotohAlpha + StrEdit_GotohBeta*((float)(i-1));
 422       E[i][0] = StrEdit_GotohAlpha + StrEdit_GotohBeta*((float)(i-1));
 423       F[i][0] = 0.;
 424    }
 425    for(j=1;j<=length2;j++) {
 426       D[0][j] = StrEdit_GotohAlpha + StrEdit_GotohBeta*((float)(j-1));
 427       E[0][j] = 0.;
 428       F[0][j] = StrEdit_GotohAlpha + StrEdit_GotohBeta*((float)(j-1));
 429    }
 430    for(i=1;i<=length1;i++) {
 431       for(j=1;j<=length2;j++) {
 432          E[i][j] = MIN2(  (D[i][j-1]+StrEdit_GotohAlpha),
 433                          (E[i][j-1]+StrEdit_GotohBeta)  );
 434          F[i][j] = MIN2(  (D[i-1][j]+StrEdit_GotohAlpha),
 435                          (F[i-1][j]+StrEdit_GotohBeta)  );
 436          D[i][j] = MIN3(  E[i][j], F[i][j],
 437                          (D[i-1][j-1]+StrEditCost(i,j,str1,str2)) );
 438       }
 439    }
 440    temp = D[length1][length2];
 441    for(i=0;i<=length1;i++) {
 442       free(D[i]); free(E[i]); free(F[i]);
 443    }
 444    free(D); free(E); free(F);
 445
 446    return temp;
 447 }
 448
 449 /* -------------------------------------------------------------------------- */
 450
 451 PRIVATE float StrEditCost(int i, int j, char *T1, char *T2)
 452 {
 453    /* positions i,j from [1..length]; i,j=0 implies Gap */
 454    int i1,j1;
 455    if((i==0)&&(j==0)) nrerror("Edit Cost: Aligned gap characters !!!");
 456    if(i>0) i1 = decode(T1[i-1]); else i1 = 0;
 457    if(j>0) j1 = decode(T2[j-1]); else j1 = 0;
 458    if(StrEdit_CostMatrix==NULL) {
 459       if(i&&j) return (float)(i1!=j1);
 460       else     return (float) StrEdit_GapCost;
 461    }
 462    else return (float) StrEdit_CostMatrix[i1][j1];
 463 }
 464
 465 /* -------------------------------------------------------------------------- */
 466
 467 PRIVATE int decode(char id)
 468 {
 469    int   n,alen;
 470
 471    if (!StrEdit_ValidAlphabet) return (int)id;
 472    alen = strlen(StrEdit_ValidAlphabet);
 473    if(!alen) return (int)id;
 474
 475    for(n=0;n<alen;n++) {
 476       if(id==StrEdit_ValidAlphabet[n]) return n;
 477    }
 478    fprintf(stderr,"Warning: Invalid character in DECODE -> set to ~gap~\n");
 479    return 0;
 480 }
 481
 482 /* -------------------------------------------------------------------------- */
 483
 484 PUBLIC void  Set_StrEdit_CostMatrix(char type)
 485 {
 486    int     i,j;
 487    if(StrEdit_ValidAlphabet) {
 488       free(StrEdit_ValidAlphabet);
 489       StrEdit_ValidAlphabet = NULL;
 490    }
 491    if(StrEdit_CostMatrix) {
 492       free(StrEdit_CostMatrix);
 493       StrEdit_CostMatrix = NULL;
 494    }
 495    switch(type){
 496       case 'D' :
 497         StrEdit_ValidAlphabet = (char*) space((20+2)*sizeof(char));
 498         strcpy(StrEdit_ValidAlphabet,StrEdit_DayhoffA);
 499         StrEdit_CostMatrix    = (float**) space((20+1)*sizeof(float*));
 500         for(i=0;i<=20;i++)
 501            StrEdit_CostMatrix[i] = (float*)space((20+1)*sizeof(float));
 502         for(i=1;i<=20;i++) {
 503            for(j=1;j<=20;j++) {
 504               StrEdit_CostMatrix[i][j] =
 505                  MAX2(StrEdit_DayhoffM[i][i],StrEdit_DayhoffM[j][j])-
 506                                                 StrEdit_DayhoffM[i][j];
 507            }
 508            StrEdit_CostMatrix[i][0] = StrEdit_DayhoffM[i][i];
 509            StrEdit_CostMatrix[0][i] = StrEdit_DayhoffM[i][i];
 510         }
 511         StrEdit_CostMatrix[0][0] = StrEdit_DayhoffM[0][0];
 512         break;
 513       case 'A' :
 514         StrEdit_ValidAlphabet = (char*) space((20+2)*sizeof(char));
 515         strcpy(StrEdit_ValidAlphabet,StrEdit_GLHA);
 516         StrEdit_CostMatrix    = (float**) space((20+1)*sizeof(float*));
 517         for(i=0;i<=20;i++)
 518            StrEdit_CostMatrix[i] = (float*)space((20+1)*sizeof(float));
 519         for(i=1;i<=20;i++) {
 520            for(j=1;j<=20;j++)
 521               StrEdit_CostMatrix[i][j] = StrEdit_GLHM[i][j];
 522            StrEdit_CostMatrix[i][0] = StrEdit_CostMatrix[0][i] =
 523               StrEdit_GapCost;
 524         }
 525         StrEdit_CostMatrix[0][0] = StrEdit_GLHM[0][0];
 526         break;
 527       case 'B' :
 528         StrEdit_ValidAlphabet = space((4+2)*sizeof(char));
 529         strcpy(StrEdit_ValidAlphabet,StrEdit_BinCodeA);
 530         StrEdit_CostMatrix    = (float**) space((4+1)*sizeof(float*));
 531         for(i=0;i<=4;i++)
 532            StrEdit_CostMatrix[i] = (float*)space((4+1)*sizeof(float));
 533         for(i=0;i<=4;i++)
 534            for(j=0;j<=4;j++)
 535               StrEdit_CostMatrix[i][j] = StrEdit_BinCodeM[i][j];
 536         break;
 537       case 'H' :
 538         StrEdit_ValidAlphabet = space((4+2)*sizeof(char));
 539         strcpy(StrEdit_ValidAlphabet,StrEdit_HogewegA);
 540         StrEdit_CostMatrix    = (float**) space((4+1)*sizeof(float*));
 541         for(i=0;i<=4;i++)
 542            StrEdit_CostMatrix[i] = (float*)space((4+1)*sizeof(float));
 543         for(i=0;i<=4;i++)
 544            for(j=0;j<=4;j++)
 545               StrEdit_CostMatrix[i][j] = StrEdit_HogewegM[i][j];
 546         StrEdit_GotohAlpha    = 3.;
 547         StrEdit_GotohBeta     = 0.;
 548         break;
 549    default:
 550         if(!StrEdit_GapCost) StrEdit_GapCost = 1.;  /* This is the simple distance */
 551    }
 552 }
 553
 554 /* -------------------------------------------------------------------------- */
 555
 556 PUBLIC   void    Set_StrEdit_GapCosts(float per_digit, float per_gap)
 557 {
 558    if(per_gap==0.) per_gap = per_digit;
 559    if(per_digit<0) nrerror("Gap Costs invalid.");
 560    if(per_digit>per_gap) nrerror("Gap Costs invalid.");
 561
 562    StrEdit_GapCost     = per_digit;
 563    StrEdit_GotohAlpha  = per_digit;   /* Gotoh gap function g(k) = a + b(k-1) */
 564    StrEdit_GotohBeta   = per_gap;
 565
 566 }
 567
 568 /* -------------------------------------------------------------------------- */