binaries/src/ViennaRNA/Cluster/treeplot.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include "distance_matrix.h"
   6 #include "utils.h"
   7
   8 #define PUBLIC
   9 #define PRIVATE  static
  10
  11 #define MAXTAXA_FOR_LABELS  200
  12
  13 typedef struct{
  14         int   set1;
  15         int   set2;
  16         float distance;
  17         float distance2;
  18         } Union;
  19
  20 typedef struct _Node_ {
  21                  float  height;
  22                  float  brr;
  23                  float  brl;
  24                  int    whoami;
  25                  int    size;
  26                  struct _Node_  *father;
  27                  struct _Node_  *left;
  28                  struct _Node_  *right;
  29                } Node;
  30
  31 PUBLIC void  PSplot_phylogeny(Union *cluster, char *filename, char *type);
  32
  33 PRIVATE Node *W2Phylo(Union *cluster);
  34 PRIVATE Node *Nj2Phylo(Union *cluster);
  35 PRIVATE void  free_phylo_tree(Node *root);
  36 PRIVATE void fill_height_into_tree(Node *root);
  37 PRIVATE void fill_br_from_height(Node *root);
  38 PRIVATE void plot_branch(Node *root, FILE *fp);
  39 PRIVATE void format_number(float x);
  40
  41 PRIVATE char  str[10];
  42 PRIVATE float threshold;
  43 PRIVATE int   print_labels=1;
  44
  45 /* --------------------------------------------------------------------------*/
  46
  47 PUBLIC void  PSplot_phylogeny(Union *cluster, char *filename, char *type)
  48 {
  49    int   n;
  50    char  outfile[50];
  51    char  *tmp, *tfont;
  52    Node  *root;
  53    Node  *tempnode;
  54    FILE  *fp;
  55    float xsize, ysize, tfontsize, lfontsize, lwidth;
  56
  57    n = cluster[0].set1;
  58    switch(type[0]) {
  59       case 'W' :
  60         root = W2Phylo(cluster);
  61         break;
  62       case 'N' :
  63         root = Nj2Phylo(cluster);
  64         break;
  65       default :
  66         return;
  67    }
  68
  69    outfile[0] = '\0';
  70    tmp = get_taxon_label(-1);   /* retrieve data set identifier */
  71    if(tmp) {strcat(outfile,tmp); strcat(outfile,"_"); free(tmp);}
  72    strcat(outfile,filename);
  73
  74    fp = fopen(outfile,"w");
  75    if (fp==NULL) {
  76       fprintf(stderr,"couldn't open %s -- not doing treeplot\n", outfile);
  77       return;
  78    }
  79    xsize = root->size;
  80    tempnode = root;
  81    while(tempnode->right)
  82       tempnode = tempnode->right;
  83    ysize = tempnode->height;
  84    threshold = ysize*0.06;
  85
  86    lwidth = MIN2(1.5, 8/sqrt((double)n));
  87    tfontsize = MIN2(15,550./(10+n));
  88    lfontsize = 2./3.*tfontsize;
  89    if (n>30) tfont = "Helvetica";
  90    else tfont = "Times-Roman";
  91    if(n > MAXTAXA_FOR_LABELS) print_labels=0;
  92    else print_labels=1;
  93    fprintf(fp,"%%!PS-Adobe-2.0 EPSF-1.2\n");
  94    fprintf(fp,"%%%%Title: TreePlot (%s)\n",type);
  95    fprintf(fp,"%%%%Creator: AnalyseDists\n");
  96    fprintf(fp,"%%%%CreationDate: %s", time_stamp());
  97    /* BoundingBox is only approximate */
  98    fprintf(fp,"%%%%BoundingBox: 35 45 535 640\n");
  99    fprintf(fp,"%%%%Pages: 1\n");
 100    fprintf(fp,"%%%%EndComments: No comment!\n");
 101    fprintf(fp,"288.5 50 translate\n");
 102    fprintf(fp,"%3.1f setlinewidth\n"
 103               "/cmtx matrix currentmatrix def\n", lwidth);
 104    fprintf(fp,"500 %g div 360 %g div scale\n", xsize, ysize);
 105    fprintf(fp,"/rotshow {gsave cmtx setmatrix\n");
 106    fprintf(fp,"          90 rotate 5 %4.1f rmoveto show grestore} def\n",
 107            -tfontsize/4);
 108    fprintf(fp,"/cshow {gsave cmtx setmatrix\n"
 109               "          /Helvetica findfont %4.1f scalefont setfont\n"
 110               "          90 rotate 0 %3.1f rmoveto\n"
 111               "          dup stringwidth pop 2 div neg 0 rmoveto show\n"
 112               "        grestore} def\n", lfontsize, lfontsize/5);
 113    fprintf(fp,"/%s findfont %4.1f scalefont setfont\n", tfont, tfontsize);
 114
 115    fprintf(fp,"0 0 moveto\n");
 116    plot_branch(root,fp);
 117
 118    fprintf(fp,"cmtx setmatrix stroke\nshowpage\n");
 119    free_phylo_tree(root);
 120    fclose(fp);
 121
 122 }
 123
 124 /* --------------------------------------------------------------------------*/
 125
 126 PRIVATE Node *W2Phylo(Union *cluster)
 127 {
 128    int i,n;
 129    float b;
 130    Node **taxa;
 131    Node  *father;
 132
 133    n=cluster[0].set1;
 134    taxa = (Node **) space(sizeof(Node*)*(n+1));
 135
 136    b = sqrt(MAX2(cluster[n-1].distance,0.));
 137    for(i=1;i<=n;i++){
 138       taxa[i] = (Node *) space(sizeof(Node));
 139       taxa[i]->whoami = i;
 140       taxa[i]->size   = 1;
 141       taxa[i]->height = b;
 142    }
 143
 144    for(i=1;i<n;i++) {
 145       father = (Node *) space(sizeof(Node));
 146       father->whoami = 0;
 147       father->left  = taxa[cluster[i].set1];
 148       father->right = taxa[cluster[i].set2];
 149       father->size  = father->left->size+father->right->size;
 150       b = sqrt(MAX2(cluster[n-1].distance-cluster[i].distance,0.));
 151       father->height   = b;
 152       father->left->father   = father;
 153       father->right->father  = father;
 154       taxa[cluster[i].set1] = father;
 155    }
 156    free(taxa);
 157    father->whoami = -1;
 158    fill_br_from_height(father);
 159    return father;
 160 }
 161
 162
 163 /* --------------------------------------------------------------------------*/
 164
 165 PRIVATE void fill_br_from_height(Node *root)
 166 {
 167    root->brr  = root->right->height - root->height;
 168    root->brl  = root->left->height  - root->height;
 169    if(!root->left->whoami)  fill_br_from_height(root->left);
 170    if(!root->right->whoami)  fill_br_from_height(root->right);
 171    return;
 172 }
 173
 174 /* -------------------------------------------------------------------------- */
 175 #define ORDER(X) \
 176 if ((X)->brl>(X)->brr) { \
 177    xnode=(X)->right; \
 178    (X)->right=(X)->left; \
 179    (X)->left=xnode; \
 180    br = (X)->brr; \
 181    (X)->brr=(X)->brl; \
 182    (X)->brl=br; \
 183 }
 184
 185
 186
 187 PRIVATE Node *Nj2Phylo(Union *cluster)
 188 {
 189    int i,n, br;
 190    float h1,h2,maxdist,dist;
 191    Node **taxa;
 192    Node  *father;
 193    Node  *topnode;   /* of the longest path in the tree */
 194    Node  *tempnode, *xnode;
 195    Node  *root;
 196
 197    maxdist = 0.;
 198    n=cluster[0].set1;
 199    taxa = (Node **) space(sizeof(Node*)*(n+1));
 200
 201    for(i=1;i<=n;i++){
 202       taxa[i] = (Node *) space(sizeof(Node));
 203       taxa[i]->whoami = i;
 204       taxa[i]->size   = 1;
 205    }
 206
 207    for(i=1;i<n;i++) {
 208       father = (Node *) space(sizeof(Node));
 209       father->whoami = 0;
 210       h1 = cluster[i].distance +taxa[cluster[i].set1]->height;
 211       h2 = cluster[i].distance2+taxa[cluster[i].set2]->height;
 212       dist = h1 + h2;
 213       if(dist > maxdist) {
 214          maxdist = dist;
 215          topnode = father;
 216       }
 217       if(h1<=h2) {
 218          father->brr   = cluster[i].distance2;
 219          father->brl   = cluster[i].distance;
 220          father->left  = taxa[cluster[i].set1];
 221          father->right = taxa[cluster[i].set2];
 222       }
 223       else {
 224          father->brr   = cluster[i].distance;
 225          father->brl   = cluster[i].distance2;
 226          father->left  = taxa[cluster[i].set2];
 227          father->right = taxa[cluster[i].set1];
 228       }
 229       father->height        = MAX2(h1,h2);
 230       father->size          = father->left->size + father->right->size;
 231       father->left->father   = father;
 232       father->right->father  = father;
 233       taxa[cluster[i].set1] = father;
 234    }
 235
 236    /* new root is inserted in the middle of the longest path */
 237    dist = maxdist/2.;
 238    tempnode = topnode;
 239    while (tempnode->height > dist)
 240       tempnode = tempnode->right;
 241
 242    topnode  = tempnode->father;
 243    /* new root must go between topnode and tempnode */
 244    h1 = dist - tempnode->height;
 245    h2 = topnode->height - dist ;
 246
 247    root = (Node *) space(sizeof(Node));
 248    root->father = NULL;
 249    root->whoami = -1;
 250    root->right  = tempnode;
 251    root->left   = topnode;
 252    root->height = 0.0;
 253    root->brr    = h1;
 254    root->brl    = h2;
 255    root->size   = n;
 256
 257    topnode->size = n - tempnode->size;
 258
 259    tempnode = root;
 260    while(topnode->father) {
 261       topnode->right  = topnode->father;
 262       topnode->brr    = topnode->father->brr;
 263       topnode->father = tempnode;
 264       tempnode = topnode;
 265       topnode  = topnode->right;
 266       topnode->size   = tempnode->size - tempnode->left->size;
 267       ORDER(tempnode);
 268    }
 269    /* topnode is now old root, remove it */
 270    if (tempnode->right==topnode) {
 271       tempnode->right = topnode->left;
 272       tempnode->brr   = topnode->brr + topnode->brl;
 273       tempnode->right->father = tempnode;
 274    } else {
 275       tempnode->left = topnode->left;
 276       tempnode->brl  = topnode->brr + topnode->brl;
 277       tempnode->left->father = tempnode;
 278    }
 279    free(topnode);
 280    ORDER(root);
 281
 282    fill_height_into_tree(root);
 283
 284    return root;
 285 }
 286
 287 /* -------------------------------------------------------------------------- */
 288
 289 PRIVATE void fill_height_into_tree(Node *root)
 290 {
 291    if(root->whoami>0) return;
 292    root->left->height  = root->height + root->brl;
 293    root->right->height = root->height + root->brr;
 294    fill_height_into_tree(root->left);
 295    fill_height_into_tree(root->right);
 296    return;
 297 }
 298
 299 /* -------------------------------------------------------------------------- */
 300
 301 PRIVATE void free_phylo_tree(Node *root) {
 302    if(root->left)  free_phylo_tree(root->left);
 303    if(root->right) free_phylo_tree(root->right);
 304    free(root);
 305 }
 306
 307
 308 /* -------------------------------------------------------------------------- */
 309
 310 PRIVATE void plot_branch(Node *root, FILE *fp)
 311 {
 312    char *label;
 313
 314    fprintf(fp,"currentpoint %g 0  rlineto 0 %g rlineto \n",
 315                      -(float)(root->left->size)/2., root->brr);
 316    if((print_labels)&&(root->brr > threshold)) {
 317       format_number(root->brr);
 318       fprintf(fp,"currentpoint 0 %g rlineto \n", -root->brr/2.);
 319       fprintf(fp,"(%s) cshow moveto\n", str);
 320    }
 321    if(root->right->whoami==0) plot_branch(root->right,fp);
 322    else if(print_labels) {
 323       label = get_taxon_label(root->right->whoami);
 324       fprintf(fp, "(%s) rotshow\n", label);
 325       free(label);
 326    }
 327    fprintf(fp,"moveto\n");
 328    fprintf(fp,"currentpoint %g 0  rlineto 0 %g rlineto \n",
 329                      +(float)(root->right->size)/2., root->brl);
 330    if((print_labels)&&(root->brl > threshold)) {
 331       format_number(root->brl);
 332       fprintf(fp,"currentpoint 0 %g rlineto \n", -root->brl/2.);
 333       fprintf(fp,"(%s) cshow moveto\n",str);
 334    }
 335    if(root->left->whoami==0) plot_branch(root->left,fp);
 336    else if(print_labels) {
 337       label = get_taxon_label(root->left->whoami);
 338       fprintf(fp, "(%s) rotshow\n", label);
 339       free(label);
 340    }
 341    fprintf(fp,"moveto\n");
 342 }
 343
 344 /* -------------------------------------------------------------------------- */
 345
 346 PRIVATE void format_number(float x)
 347 {
 348    if(x >= 1.e5) { sprintf(str, "%4.1g",x ); return;}
 349    if(x >=  100) { sprintf(str, "%.0f" ,x ); return;}
 350    if(x >=   10) { sprintf(str, "%.1f", x ); return;}
 351    if(x >=    1) { sprintf(str, "%.2f", x ); return;}
 352    if(x >=  0.1) { sprintf(str, "%.2f", x ); return;}
 353    if(x >= 0.01) { sprintf(str, "%.3f", x ); return;}
 354    sprintf(str, "%4.1g",x );
 355    return;
 356 }
 357
 358 /* -------------------------------------------------------------------------- */