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