src/jalview/analysis/TreeBuilder.java

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