src/jalview/analysis/TreeBuilder.java

   1 package jalview.analysis;
   2
   3 import jalview.api.analysis.ScoreModelI;
   4 import jalview.api.analysis.SimilarityParamsI;
   5 import jalview.datamodel.AlignmentView;
   6 import jalview.datamodel.CigarArray;
   7 import jalview.datamodel.SeqCigar;
   8 import jalview.datamodel.SequenceI;
   9 import jalview.datamodel.SequenceNode;
  10 import jalview.math.MatrixI;
  11 import jalview.viewmodel.AlignmentViewport;
  12
  13 import java.util.BitSet;
  14 import java.util.Vector;
  15
  16 public abstract class TreeBuilder
  17 {
  18   public static final String AVERAGE_DISTANCE = "AV";
  19
  20   public static final String NEIGHBOUR_JOINING = "NJ";
  21
  22   protected Vector<BitSet> clusters;
  23
  24   protected SequenceI[] sequences;
  25
  26   public AlignmentView seqData;
  27
  28   protected BitSet done;
  29
  30   protected int noseqs;
  31
  32   int noClus;
  33
  34   protected MatrixI distances;
  35
  36   protected int mini;
  37
  38   protected int minj;
  39
  40   protected double ri;
  41
  42   protected double rj;
  43
  44   SequenceNode maxdist;
  45
  46   SequenceNode top;
  47
  48   double maxDistValue;
  49
  50   double maxheight;
  51
  52   int ycount;
  53
  54   Vector<SequenceNode> node;
  55
  56   private AlignmentView seqStrings;
  57
  58   /**
  59    * Constructor
  60    *
  61    * @param av
  62    * @param sm
  63    * @param scoreParameters
  64    */
  65   public TreeBuilder(AlignmentViewport av, ScoreModelI sm,
  66           SimilarityParamsI scoreParameters)
  67   {
  68     int start, end;
  69     boolean selview = av.getSelectionGroup() != null
  70             && av.getSelectionGroup().getSize() > 1;
  71     seqStrings = av.getAlignmentView(selview);
  72     if (!selview)
  73     {
  74       start = 0;
  75       end = av.getAlignment().getWidth();
  76       this.sequences = av.getAlignment().getSequencesArray();
  77     }
  78     else
  79     {
  80       start = av.getSelectionGroup().getStartRes();
  81       end = av.getSelectionGroup().getEndRes() + 1;
  82       this.sequences = av.getSelectionGroup()
  83               .getSequencesInOrder(av.getAlignment());
  84     }
  85
  86     init(seqStrings, start, end);
  87
  88     computeTree(sm, scoreParameters);
  89   }
  90
  91   public SequenceI[] getSequences()
  92   {
  93     return sequences;
  94   }
  95
  96   /**
  97    * DOCUMENT ME!
  98    *
  99    * @param nd
 100    *          DOCUMENT ME!
 101    *
 102    * @return DOCUMENT ME!
 103    */
 104   double findHeight(SequenceNode nd)
 105   {
 106     if (nd == null)
 107     {
 108       return maxheight;
 109     }
 110
 111     if ((nd.left() == null) && (nd.right() == null))
 112     {
 113       nd.height = ((SequenceNode) nd.parent()).height + nd.dist;
 114
 115       if (nd.height > maxheight)
 116       {
 117         return nd.height;
 118       }
 119       else
 120       {
 121         return maxheight;
 122       }
 123     }
 124     else
 125     {
 126       if (nd.parent() != null)
 127       {
 128         nd.height = ((SequenceNode) nd.parent()).height + nd.dist;
 129       }
 130       else
 131       {
 132         maxheight = 0;
 133         nd.height = (float) 0.0;
 134       }
 135
 136       maxheight = findHeight((SequenceNode) (nd.left()));
 137       maxheight = findHeight((SequenceNode) (nd.right()));
 138     }
 139
 140     return maxheight;
 141   }
 142
 143   /**
 144    * DOCUMENT ME!
 145    *
 146    * @param nd
 147    *          DOCUMENT ME!
 148    */
 149   void reCount(SequenceNode nd)
 150   {
 151     ycount = 0;
 152     // _lycount = 0;
 153     // _lylimit = this.node.size();
 154     _reCount(nd);
 155   }
 156
 157   /**
 158    * DOCUMENT ME!
 159    *
 160    * @param nd
 161    *          DOCUMENT ME!
 162    */
 163   void _reCount(SequenceNode nd)
 164   {
 165     // if (_lycount<_lylimit)
 166     // {
 167     // System.err.println("Warning: depth of _recount greater than number of
 168     // nodes.");
 169     // }
 170     if (nd == null)
 171     {
 172       return;
 173     }
 174     // _lycount++;
 175
 176     if ((nd.left() != null) && (nd.right() != null))
 177     {
 178
 179       _reCount((SequenceNode) nd.left());
 180       _reCount((SequenceNode) nd.right());
 181
 182       SequenceNode l = (SequenceNode) nd.left();
 183       SequenceNode r = (SequenceNode) nd.right();
 184
 185       nd.count = l.count + r.count;
 186       nd.ycount = (l.ycount + r.ycount) / 2;
 187     }
 188     else
 189     {
 190       nd.count = 1;
 191       nd.ycount = ycount++;
 192     }
 193     // _lycount--;
 194   }
 195
 196   /**
 197    * DOCUMENT ME!
 198    *
 199    * @return DOCUMENT ME!
 200    */
 201   public SequenceNode getTopNode()
 202   {
 203     return top;
 204   }
 205
 206   /**
 207    *
 208    * @return true if tree has real distances
 209    */
 210   public boolean hasDistances()
 211   {
 212     return true;
 213   }
 214
 215   /**
 216    *
 217    * @return true if tree has real bootstrap values
 218    */
 219   public boolean hasBootstrap()
 220   {
 221     return false;
 222   }
 223
 224   public boolean hasRootDistance()
 225   {
 226     return true;
 227   }
 228
 229   /**
 230    * Form clusters by grouping sub-clusters, starting from one sequence per
 231    * cluster, and finishing when only two clusters remain
 232    */
 233   void cluster()
 234   {
 235     while (noClus > 2)
 236     {
 237       findMinDistance();
 238
 239       joinClusters(mini, minj);
 240
 241       noClus--;
 242     }
 243
 244     int rightChild = done.nextClearBit(0);
 245     int leftChild = done.nextClearBit(rightChild + 1);
 246
 247     joinClusters(leftChild, rightChild);
 248     top = (node.elementAt(leftChild));
 249
 250     reCount(top);
 251     findHeight(top);
 252     findMaxDist(top);
 253   }
 254
 255   /**
 256    * Returns the minimum distance between two clusters, and also sets the
 257    * indices of the clusters in fields mini and minj
 258    *
 259    * @return
 260    */
 261   protected abstract double findMinDistance();
 262
 263   /**
 264    * Calculates the tree using the given score model and parameters, and the
 265    * configured tree type
 266    * <p>
 267    * If the score model computes pairwise distance scores, then these are used
 268    * directly to derive the tree
 269    * <p>
 270    * If the score model computes similarity scores, then the range of the scores
 271    * is reversed to give a distance measure, and this is used to derive the tree
 272    *
 273    * @param sm
 274    * @param scoreOptions
 275    */
 276   protected void computeTree(ScoreModelI sm, SimilarityParamsI scoreOptions)
 277   {
 278     distances = sm.findDistances(seqData, scoreOptions);
 279
 280     makeLeaves();
 281
 282     noClus = clusters.size();
 283
 284     cluster();
 285   }
 286
 287   /**
 288    * Finds the node, at or below the given node, with the maximum distance, and
 289    * saves the node and the distance value
 290    *
 291    * @param nd
 292    */
 293   void findMaxDist(SequenceNode nd)
 294   {
 295     if (nd == null)
 296     {
 297       return;
 298     }
 299
 300     if ((nd.left() == null) && (nd.right() == null))
 301     {
 302       double dist = nd.dist;
 303
 304       if (dist > maxDistValue)
 305       {
 306         maxdist = nd;
 307         maxDistValue = dist;
 308       }
 309     }
 310     else
 311     {
 312       findMaxDist((SequenceNode) nd.left());
 313       findMaxDist((SequenceNode) nd.right());
 314     }
 315   }
 316
 317   /**
 318    * Calculates and returns r, whatever that is
 319    *
 320    * @param i
 321    * @param j
 322    *
 323    * @return
 324    */
 325   protected double findr(int i, int j)
 326   {
 327     double tmp = 1;
 328
 329     for (int k = 0; k < noseqs; k++)
 330     {
 331       if ((k != i) && (k != j) && (!done.get(k)))
 332       {
 333         tmp = tmp + distances.getValue(i, k);
 334       }
 335     }
 336
 337     if (noClus > 2)
 338     {
 339       tmp = tmp / (noClus - 2);
 340     }
 341
 342     return tmp;
 343   }
 344
 345   protected void init(AlignmentView seqView, int start, int end)
 346   {
 347     this.node = new Vector<SequenceNode>();
 348     if (seqView != null)
 349     {
 350       this.seqData = seqView;
 351     }
 352     else
 353     {
 354       SeqCigar[] seqs = new SeqCigar[sequences.length];
 355       for (int i = 0; i < sequences.length; i++)
 356       {
 357         seqs[i] = new SeqCigar(sequences[i], start, end);
 358       }
 359       CigarArray sdata = new CigarArray(seqs);
 360       sdata.addOperation(CigarArray.M, end - start + 1);
 361       this.seqData = new AlignmentView(sdata, start);
 362     }
 363
 364     /*
 365      * count the non-null sequences
 366      */
 367     noseqs = 0;
 368
 369     done = new BitSet();
 370
 371     for (SequenceI seq : sequences)
 372     {
 373       if (seq != null)
 374       {
 375         noseqs++;
 376       }
 377     }
 378   }
 379
 380   /**
 381    * Merges cluster(j) to cluster(i) and recalculates cluster and node distances
 382    *
 383    * @param i
 384    * @param j
 385    */
 386   void joinClusters(final int i, final int j)
 387   {
 388     double dist = distances.getValue(i, j);
 389
 390     ri = findr(i, j);
 391     rj = findr(j, i);
 392
 393     findClusterDistance(i, j);
 394
 395     SequenceNode sn = new SequenceNode();
 396
 397     sn.setLeft((node.elementAt(i)));
 398     sn.setRight((node.elementAt(j)));
 399
 400     SequenceNode tmpi = (node.elementAt(i));
 401     SequenceNode tmpj = (node.elementAt(j));
 402
 403     findNewDistances(tmpi, tmpj, dist);
 404
 405     tmpi.setParent(sn);
 406     tmpj.setParent(sn);
 407
 408     node.setElementAt(sn, i);
 409
 410     /*
 411      * move the members of cluster(j) to cluster(i)
 412      * and mark cluster j as out of the game
 413      */
 414     clusters.get(i).or(clusters.get(j));
 415     clusters.get(j).clear();
 416     done.set(j);
 417   }
 418
 419   /*
 420    * Computes and stores new distances for nodei and nodej, given the previous
 421    * distance between them
 422    */
 423   protected abstract void findNewDistances(SequenceNode nodei,
 424           SequenceNode nodej, double previousDistance);
 425
 426   /**
 427    * Calculates and saves the distance between the combination of cluster(i) and
 428    * cluster(j) and all other clusters. The form of the calculation depends on
 429    * the tree clustering method being used.
 430    *
 431    * @param i
 432    * @param j
 433    */
 434   protected abstract void findClusterDistance(int i, int j);
 435
 436   /**
 437    * Start by making a cluster for each individual sequence
 438    */
 439   void makeLeaves()
 440   {
 441     clusters = new Vector<BitSet>();
 442
 443     for (int i = 0; i < noseqs; i++)
 444     {
 445       SequenceNode sn = new SequenceNode();
 446
 447       sn.setElement(sequences[i]);
 448       sn.setName(sequences[i].getName());
 449       node.addElement(sn);
 450       BitSet bs = new BitSet();
 451       bs.set(i);
 452       clusters.addElement(bs);
 453     }
 454   }
 455
 456   public AlignmentView getOriginalData()
 457   {
 458     return seqStrings;
 459   }
 460
 461 }