binaries/src/ViennaRNA/Cluster/statgeom.c

   1 /*      Quadruple Statistics on an arbitrary alphabet
   2
   3               Vienna RNA Package   ---  Peter F Stadler   1993
   4
   5 */
   6
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <strings.h>
  10 #include <string.h>
  11 #include <ctype.h>
  12 #include "utils.h"
  13 #include "PS3D.h"
  14 #include "distance_matrix.h"
  15
  16 #define PUBLIC
  17 #define PRIVATE    static
  18
  19 #define MIN2(A, B)         ((A) < (B) ? (A) : (B))
  20 #define MIN4(A, B, C, D)   MIN2( (MIN2((A),(B))), (MIN2((C),(D))) )
  21 #define MAX2(A, B)         ((A) > (B) ? (A) : (B))
  22 #define MAX4(A, B, C, D)   MAX2( (MAX2((A),(B))), (MAX2((C),(D))) )
  23
  24 PRIVATE int IBox[16];
  25
  26
  27 PUBLIC float *statgeom(char **seqs, int n_of_seqs);
  28 PUBLIC float *statgeom4(char **ss[4], int nn[4]);
  29 PUBLIC void   printf_stg(float *B);
  30
  31 PRIVATE void SingleBox(char *x1, char *x2, char *x3, char *x4);
  32 PRIVATE void SortSingleBox(void);
  33
  34 /* ----------------------------------------------------------------------- */
  35
  36 PUBLIC float *statgeom(char **seqs, int n_of_seqs)
  37 {
  38    int i,j,k,l;
  39    int i1;
  40
  41    float *B;
  42  float  temp;
  43
  44    if(n_of_seqs < 4) {
  45       fprintf(stderr,"Less than 4 sequences for statistical geometry.\n");
  46       return NULL;
  47    }
  48
  49    B = (float *) space(16*sizeof(float));
  50
  51    for(i=3; i<n_of_seqs; i++) {
  52     for(j=2; j<i; j++) {
  53      for(k=1; k<j; k++) {
  54       for(l=0; l<k; l++) {
  55           SingleBox(seqs[i], seqs[j], seqs[k], seqs[l]);
  56           SortSingleBox();
  57           B[0] = (float) IBox[0];      /* transfer length */
  58           for(i1=1;i1<=15;i1++) B[i1] += ( (float) IBox[i1] );
  59       }
  60      }
  61     }
  62    }
  63
  64    temp = (float) n_of_seqs;
  65    temp = temp*(temp-1.)*(temp-2.)*(temp-3.)/24.;
  66    for(i1=1;i1<=15;i1++) B[i1] /= (temp*B[0]);
  67
  68    return B;
  69 }
  70
  71 /* ----------------------------------------------------------------------- */
  72
  73 PUBLIC float *statgeom4(char **ss[4], int  nn[4])
  74 {
  75    int i,j,k,l;
  76    int i1;
  77
  78    float *B;
  79    float  temp;
  80
  81    B = (float *) space(16*sizeof(float));
  82
  83    for(i=0; i<nn[0]; i++) {
  84     for(j=0; j<nn[1]; j++) {
  85      for(k=0; k<nn[2]; k++) {
  86       for(l=0; l<nn[3]; l++) {
  87           SingleBox(ss[0][i], ss[1][j], ss[2][k], ss[3][l]);
  88           B[0] = (float) IBox[0];      /* transfer length */
  89           for(i1=1;i1<=15;i1++) B[i1] += ( (float) IBox[i1] );
  90       }
  91      }
  92     }
  93    }
  94
  95    for (temp = 1, i=0;i<4;i++) temp *= ((float) nn[i]) ;
  96    for(i1=1;i1<=15;i1++) B[i1] /= (temp*B[0]);
  97
  98    return B;
  99 }
 100 /* ------------------------------------------------------------------------- */
 101
 102 PUBLIC void printf_stg(float *B)
 103 {
 104    printf("> Statistical Geometry.\n");
 105    printf("> %d (sequence length)\n", (int) B[0]);
 106    printf("> AAAA\n");
 107    printf("  %7.5f\n", B[1]);
 108    printf("> BAAA        ABAA        AABA        AAAB\n");
 109    printf("  %7.5f    %7.5f     %7.5f     %7.5f\n", B[2],B[3],B[4],B[5]);
 110    printf("> AABB        ABAB        ABBA\n");
 111    printf("  %7.5f    %7.5f     %7.5f\n", B[6],B[7],B[8]);
 112    printf("> AABC        ABAC        ABCA        BAAC        BACA        BCAA\n");
 113    printf("  %7.5f    %7.5f     %7.5f     %7.5f    %7.5f    %7.5f\n"
 114             ,B[9],B[10],B[11],B[12],B[13],B[14]);
 115    printf("> ABCD\n");
 116    printf("  %7.5f\n",B[15]);
 117 }
 118
 119 /* ------------------------------------------------------------------------- */
 120
 121
 122 PUBLIC void SimplifiedBox(float *B, char *filename)
 123 {
 124    char *tmp, temp[50];
 125    float X,Y,Z;
 126    float T,P;
 127    float t1, t2, t3, t4;
 128    float p1, p2, p3, p4;
 129    float x1, x2, y1, y2,z1,z2;
 130    FILE *fp;
 131    float view[3] = {0.3, 1.5, 0.1};
 132    float axis[3] = {1.0, 0.0, 0.0};
 133
 134    t1 = B[9]+B[10]+B[11]+B[15];
 135    t2 = B[9]+B[12]+B[13]+B[15];
 136    t3 = B[10]+B[12]+B[14]+B[15];
 137    t4 = B[11]+B[13]+B[14]+B[15];
 138
 139    p1 = B[2]; p2 = B[3]; p3 = B[4]; p4 = B[5];
 140
 141    x1 = B[6] + B[14];
 142    x2 = B[6] + B[9];
 143    y1 = B[7] + B[13];
 144    y2 = B[7] + B[10];
 145    z1 = B[8] + B[12];
 146    z2 = B[8] + B[11];
 147
 148    X = (x1+x2)/2.;  Y = (y1+y2)/2.;  Z = (z1+z2)/2;
 149    T = (t1+t2+t3+t4)/4.;
 150    P = (p1+p2+p3+p4)/4.;
 151
 152    printf("> X = %7.5f \n",X);
 153    printf("> Y = %7.5f \n",Y);
 154    printf("> Z = %7.5f \n",Z);
 155    printf("> T = %7.5f \n",T);
 156    printf("> P = %7.5f \n",P);
 157
 158    tmp = get_taxon_label(-1);      /* retrieve the dataset identifier */
 159
 160    temp[0]='\0';
 161    if(tmp) { strcat(temp,tmp);strcat(temp,"_"); free(tmp); }
 162    strcat(temp,filename);
 163
 164    fp = fopen(temp,"w");
 165    if (fp!=NULL) {
 166       ps3d_Preambel(fp, view, axis, "N");
 167       PS_DrawSimplifiedBox(X,Y,Z,T,P, fp);
 168       fclose(fp);
 169    } else fprintf(stderr,"couldn't open %s -- not drawing box\n", temp);
 170 }
 171
 172
 173
 174 /* ------------------------------------------------------------------------- */
 175
 176 PRIVATE void SingleBox(char *x1, char *x2, char *x3, char *x4)
 177 {
 178    int   len,i1,j1,k1,i,M,m;
 179    int   d[4];
 180    char  t[4];
 181
 182    len=strlen(x1);
 183    if(strlen(x2)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
 184    if(strlen(x3)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
 185    if(strlen(x4)!=len) nrerror("Sequences of unequal length in 'SingleBox'");
 186
 187    IBox[0] = len;
 188    for(i=1; i<=15; i++) IBox[i] = 0;
 189
 190    for(i1=0;i1<len;i1++) {
 191       t[0] = x1[i1];
 192       t[1] = x2[i1];
 193       t[2] = x3[i1];
 194       t[3] = x4[i1];
 195       for(j1=0;j1<4;j1++) {
 196          d[j1]=4;
 197          for(k1=0;k1<4;k1++) d[j1] -=(t[j1]!=t[k1]);
 198       }
 199       M=MAX4(d[0],d[1],d[2],d[3]);
 200       switch(M) {
 201        case 4 :     /* Four of a kind */
 202          IBox[1]++;                            /*  A A A A */
 203          break;
 204        case 3 :     /* Three of a kind */
 205          if(d[0] == 1)      IBox[2]++;         /*  B A A A */
 206          else if(d[1] == 1) IBox[3]++;         /*  A B A A */
 207          else if(d[2] == 1) IBox[4]++;         /*  A A B A */
 208          else               IBox[5]++;         /*  A A A B */
 209          break;
 210        case 2 :
 211          m=MIN4(d[0],d[1],d[2],d[3]);
 212          if(m==2){   /* Two Pairs */
 213             if(t[1]==t[0])       IBox[6]++;    /*  A A B B */
 214             else if(t[2]==t[0])  IBox[7]++;    /*  A B A B */
 215             else                 IBox[8]++;    /*  A B B A */
 216          }
 217          else {      /* One Pair */
 218             if(d[0]==2){  /* 0 is in the pair */
 219                if(t[1]==t[0])      IBox[9]++;  /*  A A B C */
 220                else if(t[2]==t[0]) IBox[10]++; /*  A B A C */
 221                else                IBox[11]++ ;/*  A B C A */
 222             }
 223             else if(d[1]==2){  /* 1 is in the pair */
 224                if(t[2]==t[1])      IBox[12]++; /*  B A A C */
 225                else                IBox[13]++; /*  B A C A */
 226             }
 227             else                   IBox[14]++; /*  B C A A */
 228          }
 229          break;
 230        case 1 :      /* No Pair */
 231          IBox[15]++;                           /*  A B C D */
 232          break;
 233        default:
 234          nrerror("This can't happen.");
 235       }
 236    }
 237 }
 238
 239 /* ----------------------------------------------------------------------- */
 240
 241 PRIVATE void SortSingleBox(void)
 242 {
 243    int i;
 244    int M;
 245    int IBB[16];
 246    int s[4];
 247
 248    M = MAX2(MAX2( IBox[6], IBox[7] ), IBox[8] );
 249
 250    if( M== IBox[6] ) {                   /* 12|34       */
 251       IBB[6] = IBox[6];
 252       if( IBox[9] >= IBox[14] ) {        /* 1,2 > 3,4   */
 253          IBB[9]  = IBox[9];
 254          IBB[14] = IBox[14];
 255          if( IBox[7] >= IBox[8] ) {      /*    13|24    */
 256             IBB[7] = IBox[7];
 257             IBB[8] = IBox[8];
 258             if( IBox[10] >= IBox[13] ) { /* 1>2>3>4     */
 259                IBB[10] = IBox[10];
 260                IBB[13] = IBox[13];
 261                IBB[11] = IBox[11];
 262                IBB[12] = IBox[12];
 263                s[0]=0; s[1]=1; s[2]=2; s[3]=3;            /*  1 2 3 4 */
 264             }
 265             else {                       /* 2>1>4>3     */
 266                IBB[10] = IBox[13];
 267                IBB[13] = IBox[10];
 268                IBB[11] = IBox[12];
 269                IBB[12] = IBox[11];
 270                s[0]=1; s[1]=0; s[2]=3; s[3]=2;            /*  2 1 4 3 */
 271             }
 272          }
 273          else {                          /*    14|23    */
 274             IBB[7] = IBox[8];
 275             IBB[8] = IBox[7];
 276             if( IBox[11] >= IBox[12]) {   /* 1>2>4>3     */
 277                IBB[10] = IBox[11];
 278                IBB[13] = IBox[12];
 279                IBB[11] = IBox[10];
 280                IBB[12] = IBox[13];
 281                s[0]=0; s[1]=1; s[2]=3; s[3]=2;            /*  1 2 4 3 */
 282             }
 283             else {                       /* 2>1>3>4     */
 284                IBB[10] = IBox[12];
 285                IBB[13] = IBox[11];
 286                IBB[11] = IBox[13];
 287                IBB[12] = IBox[10];
 288                s[0]=1; s[1]=0; s[2]=2; s[3]=3;            /*  2 1 3 4 */
 289             }
 290          }
 291       }
 292       else {                             /* 3,4 > 1,2   */
 293          IBB[9] = IBox[14];
 294          IBB[14]= IBox[9];
 295          if( IBox[7] >= IBox[8] ) {      /*    31|42    */
 296             IBB[7] = IBox[7];
 297             IBB[8] = IBox[8];
 298             if(IBox[10] >= IBox[13]) {   /* 3>4>1>2     */
 299                IBB[10] = IBox[10];
 300                IBB[13] = IBox[13];
 301                IBB[11] = IBox[12];
 302                IBB[12] = IBox[11];
 303                s[0]=2; s[1]=3; s[2]=0; s[3]=1;            /*  3 4 1 2 */
 304             }
 305             else {                       /*            */
 306                IBB[10] = IBox[13];
 307                IBB[13] = IBox[10];
 308                IBB[11] = IBox[11];
 309                IBB[12] = IBox[12];
 310                s[0]=3; s[1]=2; s[2]=1; s[3]=0;            /*  4 3 2 1 */
 311             }
 312          }
 313          else {                          /*    32|14    */
 314             IBB[7] = IBox[8];
 315             IBB[8] = IBox[7];
 316             if( IBox[11] >= IBox[12]) {   /* 3>4>2>1     */
 317                IBB[10] = IBox[12];
 318                IBB[13] = IBox[11];
 319                IBB[11] = IBox[10];
 320                IBB[12] = IBox[13];
 321                s[0]=2; s[1]=3; s[2]=1; s[3]=0;            /*  3 4 2 1 */
 322             }
 323             else {                       /* 2>1>3>4     */
 324                IBB[10] = IBox[10];
 325                IBB[13] = IBox[13];
 326                IBB[11] = IBox[12];
 327                IBB[12] = IBox[11];
 328                s[0]=2; s[1]=3; s[2]=0; s[3]=1;            /*  3 4 1 2 */
 329             }
 330          }
 331       }
 332    }
 333    else if (M == IBox[7] ) {             /* 13|24       */
 334       IBB[6] = IBox[7];
 335       if( IBox[10] >= IBox[13] ) {       /* 1,3 > 2,4   */
 336          IBB[9]  = IBox[10];
 337          IBB[14] = IBox[13];
 338          if( IBox[6] >= IBox[8]) {       /*    12|34    */
 339             IBB[7] = IBox[6];
 340             IBB[8] = IBox[8];
 341             if( IBox[9] >= IBox[14] ) {  /*    1,2>3,4  */
 342                IBB[10] = IBox[9];
 343                IBB[13] = IBox[14];
 344                IBB[11] = IBox[11];
 345                IBB[12] = IBox[12];
 346                s[0]=0; s[1]=2; s[2]=1; s[3]=3;            /* 1 3 2 4 */
 347             }
 348             else {
 349                IBB[10] = IBox[14];
 350                IBB[13] = IBox[9];
 351                IBB[11] = IBox[12];
 352                IBB[12] = IBox[11];
 353                s[0]=2; s[1]=0; s[2]=3; s[3]=1;            /* 3 1 4 2 */
 354             }
 355          }
 356          else {                          /*    14|23   */
 357             IBB[7] = IBox[8];
 358             IBB[8] = IBox[6];
 359             if( IBox[11] >= IBox[12] ){  /*    1,4 > 2,3 */
 360                IBB[10] = IBox[11];
 361                IBB[13] = IBox[12];
 362                IBB[11] = IBox[9];
 363                IBB[12] = IBox[14];
 364                s[0]=0; s[1]=2; s[2]=3; s[3]=1;            /* 1 3 4 2 */
 365             }
 366             else {
 367                IBB[10] = IBox[12];
 368                IBB[13] = IBox[11];
 369                IBB[11] = IBox[14];
 370                IBB[12] = IBox[9];
 371                s[0]=2; s[1]=0; s[2]=1; s[3]=3;            /* 3 1 2 4 */
 372             }
 373          }
 374       }
 375       else {                             /* 2,4 > 1,3   */
 376          IBB[9]  = IBox[13];
 377          IBB[14] = IBox[10];
 378          if( IBox[6] >= IBox[8]) {       /*    21|43    */
 379             IBB[7] = IBox[6];
 380             IBB[8] = IBox[8];
 381             if( IBox[9] >= IBox[14] ) {  /*    2,1>4,3  */
 382                IBB[10] = IBox[9];
 383                IBB[13] = IBox[14];
 384                IBB[11] = IBox[12];
 385                IBB[12] = IBox[11];
 386                s[0]=1; s[1]=3; s[2]=0; s[3]=2;            /* 2 4 1 3 */
 387             }
 388             else {
 389                IBB[10] = IBox[14];
 390                IBB[13] = IBox[9];
 391                IBB[11] = IBox[11];
 392                IBB[12] = IBox[12];
 393                s[0]=3; s[1]=1; s[2]=2; s[3]=0;            /* 4 2 3 1 */
 394             }
 395          }
 396          else {                          /*    14|23   */
 397             IBB[7] = IBox[8];
 398             IBB[8] = IBox[6];
 399             if( IBox[11] >= IBox[12] ){  /*    1,4 > 2,3 */
 400                IBB[10] = IBox[11];
 401                IBB[13] = IBox[12];
 402                IBB[11] = IBox[14];
 403                IBB[12] = IBox[9];
 404                s[0]=3; s[1]=1; s[2]=0; s[3]=2;            /* 4 2 1 3 */
 405             }
 406             else {
 407                IBB[10] = IBox[12];
 408                IBB[13] = IBox[11];
 409                IBB[11] = IBox[9];
 410                IBB[12] = IBox[14];
 411                s[0]=1; s[1]=3; s[2]=2; s[3]=0;            /* 2 4 3 1 */
 412             }
 413          }
 414       }
 415    }
 416    else {                                /* 14 | 23   */
 417       IBB[6] = IBox[8];
 418       if( IBox[11] >= IBox[12] ) {       /* 1,4 > 2,3 */
 419          IBB[9]  = IBox[11];
 420          IBB[14] = IBox[12];
 421          if(IBox[6] >= IBox[7]) {        /*    12|34  */
 422             IBB[7] = IBox[6];
 423             IBB[8] = IBox[7];
 424             if(IBox[9] >= IBox[14]) {   /*    1,2>3,4 */
 425                IBB[10] = IBox[9];
 426                IBB[13] = IBox[14];
 427                IBB[11] = IBox[10];
 428                IBB[12] = IBox[13];
 429                s[0]=0; s[1]=3; s[2]=1; s[3]=2;          /* 1 4 2 3 */
 430             }
 431             else {                      /*    4,3>2,1 */
 432                IBB[10] = IBox[14];
 433                IBB[13] = IBox[9];
 434                IBB[11] = IBox[13];
 435                IBB[12] = IBox[10];
 436                s[0]=3; s[1]=0; s[2]=2; s[3]=1;          /* 4 1 3 2 */
 437             }
 438          }
 439          else {                         /*    13|24  */
 440             IBB[7] = IBox[7];
 441             IBB[8] = IBox[6];
 442             if( IBox[10] >= IBox[13]) { /*    1,3 > 2,4 */
 443                IBB[10] = IBox[10];
 444                IBB[13] = IBox[13];
 445                IBB[11] = IBox[9];
 446                IBB[12] = IBox[14];
 447                s[0]=0; s[1]=3; s[2]=2; s[3]=1;          /* 1 4 3 2 */
 448             }
 449             else {
 450                IBB[10] = IBox[13];
 451                IBB[13] = IBox[10];
 452                IBB[11] = IBox[14];
 453                IBB[12] = IBox[9];
 454                s[0]=3; s[1]=0; s[2]=1; s[3]=2;          /* 4 1 2 3 */
 455             }
 456          }
 457       }
 458       else {                             /* 2,3 > 1,4 */
 459          IBB[9]  = IBox[12];
 460          IBB[14] = IBox[11];
 461          if(IBox[6] >= IBox[7]) {       /*    34|12  */
 462             IBB[7] = IBox[6];
 463             IBB[8] = IBox[7];
 464             if(IBox[9] >= IBox[14]) {   /*    2,1>4,3 */
 465                IBB[10] = IBox[9];
 466                IBB[13] = IBox[14];
 467                IBB[11] = IBox[13];
 468                IBB[12] = IBox[10];
 469                s[0]=1; s[1]=2; s[2]=0; s[3]=3;          /* 2 3 1 4 */
 470             }
 471             else {                      /*    4,3>2,1 */
 472                IBB[10] = IBox[14];
 473                IBB[13] = IBox[9];
 474                IBB[11] = IBox[10];
 475                IBB[12] = IBox[13];
 476                s[0]=2; s[1]=1; s[2]=3; s[3]=0;          /* 3 2 4 1 */
 477             }
 478          }
 479          else {                         /*    31|42  */
 480             IBB[7] = IBox[7];
 481             IBB[8] = IBox[6];
 482             if( IBox[10] >= IBox[13]) { /*    1,3 > 2,4 */
 483                IBB[10] = IBox[10];
 484                IBB[13] = IBox[13];
 485                IBB[11] = IBox[14];
 486                IBB[12] = IBox[9];
 487                s[0]=2; s[1]=1; s[2]=0; s[3]=4;          /* 3 2 1 4 */
 488             }
 489             else {
 490                IBB[10] = IBox[13];
 491                IBB[13] = IBox[10];
 492                IBB[11] = IBox[9];
 493                IBB[12] = IBox[14];
 494                s[0]=1; s[1]=2; s[2]=3; s[3]=0;          /* 2 3 4 1 */
 495             }
 496          }
 497       }
 498    }
 499    /* HEUREKA */
 500    IBB[0] = IBox[0];
 501    IBB[1] = IBox[1];
 502    IBB[2] = IBox[2+s[0]];
 503    IBB[3] = IBox[2+s[1]];
 504    IBB[4] = IBox[2+s[2]];
 505    IBB[5] = IBox[2+s[3]];
 506    /*
 507    IBB[6...14] see above :-)
 508    */
 509    IBB[15] = IBox[15];
 510
 511    for(i=0;i<=15;i++) IBox[i] = IBB[i];
 512 }
 513
 514 /* ----------------------------------------------------------------------- */
 515