src/jalview/analysis/PaSiMap.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.bin.Console;
  26 import jalview.datamodel.AlignmentView;
  27 import jalview.datamodel.Point;
  28 import jalview.datamodel.SequenceI;
  29 import jalview.datamodel.SequenceGroup;
  30 import jalview.gui.PairwiseAlignPanel;
  31 import jalview.gui.PaSiMapPanel;
  32 import jalview.math.Matrix;
  33 import jalview.math.MatrixI;
  34 import jalview.viewmodel.AlignmentViewport;
  35
  36
  37 import java.io.PrintStream;
  38 import java.util.Hashtable;
  39 import java.util.HashMap;
  40 import java.util.HashSet;
  41 import java.util.Enumeration;
  42
  43 /**
  44  * Performs Principal Component Analysis on given sequences
  45  * @AUTHOR MorellThomas
  46  */
  47 public class PaSiMap implements Runnable
  48 {
  49   /*
  50    * inputs
  51    */
  52   final private AlignmentViewport seqs;
  53
  54   final private ScoreModelI scoreModel;
  55
  56   final private SimilarityParamsI similarityParams;
  57
  58   final private byte dim = 3;
  59
  60   /*
  61    * outputs
  62    */
  63   private PairwiseAlignPanel alignment;
  64
  65   private MatrixI pairwiseScores;
  66
  67   private MatrixI eigenMatrix;
  68
  69   /**
  70    * Constructor given the sequences to compute for, the similarity model to
  71    * use, and a set of parameters for sequence comparison
  72    *
  73    * @param sequences
  74    * @param sm
  75    * @param options
  76    */
  77   public PaSiMap(AlignmentViewport sequences, ScoreModelI sm,
  78           SimilarityParamsI options)
  79   {
  80     this.seqs = sequences;
  81     this.scoreModel = sm;
  82     this.similarityParams = options;
  83   }
  84
  85   /**
  86    * Returns Eigenvalue
  87    *
  88    * @param i
  89    *          Index of diagonal within matrix
  90    *
  91    * @return Returns value of diagonal from matrix
  92    */
  93   public double getEigenvalue(int i)
  94   {
  95     return eigenMatrix.getD()[i];
  96   }
  97
  98   /**
  99    * Returns coordinates for each datapoint
 100    *
 101    * @param l
 102    *          DOCUMENT ME!
 103    * @param n
 104    *          DOCUMENT ME!
 105    * @param mm
 106    *          DOCUMENT ME!
 107    * @param factor ~ is 1
 108    *
 109    * @return DOCUMENT ME!
 110    */
 111   public Point[] getComponents(int l, int n, int mm, float factor)
 112   {
 113     Point[] out = new Point[getHeight()];
 114
 115     for (int i = 0; i < out.length; i++)
 116     {
 117       float x = (float) component(i, l) * factor;
 118       float y = (float) component(i, n) * factor;
 119       float z = (float) component(i, mm) * factor;
 120       out[i] = new Point(x, y, z);
 121     }
 122
 123     return out;
 124   }
 125
 126   /**
 127    * DOCUMENT ME!
 128    *
 129    * @param n
 130    *          DOCUMENT ME!
 131    *
 132    * @return DOCUMENT ME!
 133    */
 134   public double[] component(int n)
 135   {
 136     // n = index of eigenvector
 137     double[] out = new double[getWidth()];
 138
 139     for (int i = 0; i < out.length; i++)
 140     {
 141       out[i] = component(i, n);
 142     }
 143
 144     return out;
 145   }
 146
 147   /**
 148    * DOCUMENT ME!
 149    *
 150    * @param row
 151    *          DOCUMENT ME!
 152    * @param n
 153    *          DOCUMENT ME!
 154    *
 155    * @return DOCUMENT ME!
 156    */
 157   double component(int row, int n)
 158   {
 159     return eigenMatrix.getValue(row, n);
 160   }
 161
 162   /**
 163    * Answers a formatted text report of the PaSiMap calculation results (matrices
 164    * and eigenvalues) suitable for display
 165    *
 166    * @return
 167    */
 168   public String getDetails()
 169   {
 170     StringBuilder sb = new StringBuilder(1024);
 171     sb.append("PaSiMap calculation using ").append(scoreModel.getName())
 172             .append(" sequence similarity matrix\n========\n\n");
 173     PrintStream ps = wrapOutputBuffer(sb);
 174
 175     /*
 176      * coordinates matrix, with D vector
 177      */
 178     sb.append(" --- Pairwise correlation coefficients ---\n");
 179     pairwiseScores.print(ps, "%8.6f ");
 180     ps.println();
 181
 182     sb.append(" --- Eigenvalues ---\n");
 183     eigenMatrix.printD(ps, "%15.4e");
 184     ps.println();
 185
 186     sb.append(" --- Coordinates ---\n");
 187     eigenMatrix.print(ps, "%8.6f ");
 188     ps.println();
 189
 190     return sb.toString();
 191   }
 192
 193   /**
 194    * Performs the PaSiMap calculation
 195    *
 196    * creates a new gui/PairwiseAlignPanel with the input sequences (AlignmentViewport)
 197    * uses analysis/AlignSeq to creatue the pairwise alignments and calculate the AlignmentScores (float for each pair)
 198    * gets all float[][] scores from the gui/PairwiseAlignPanel
 199    * checks the connections for each sequence with AlignmentViewport seqs.calculateConnectivity(float[][] scores, int dim) (from analysis/Connectivity) -- throws an Exception if insufficient
 200    * creates a math/MatrixI pairwiseScores of the float[][] scores
 201    * copys the scores and fills the diagonal to create a symmetric matrix using math/Matrix.fillDiagonal()
 202    * performs the analysis/ccAnalysis with the symmetric matrix
 203    * gets the eigenmatrix and the eigenvalues using math/Matrix.tqli()
 204    */
 205   @Override
 206   public void run()
 207   {
 208     try
 209     {
 210 for (SequenceI see : seqs.getAlignment().getSequencesArray())
 211 {
 212 System.out.println(see.getName());
 213 }
 214
 215       int nSeqs = seqs.getAlignment().getHeight();
 216       float[][] scores = new float[nSeqs][nSeqs];       // rows, cols
 217
 218       int nSplits = 1;
 219       //while ((nSeqs / nSplits) > 300)         // heap full at 341
 220       while (((float) nSeqs / nSplits) > 5f)            // heap full at 341
 221         nSplits++;
 222       int splitSeqs = (int) Math.ceil((float) nSeqs / nSplits);
 223 System.out.println(String.format("%d -> %d splits into %d seqs", nSeqs, nSplits, splitSeqs));
 224
 225       int[] splitIndices = new int[nSplits];
 226       for (int i = 0; i < nSplits; i++)
 227       {
 228         splitIndices[i] = splitSeqs * (i + 1);  //exclusive!!
 229       }
 230
 231       HashMap<int[], Float> valuesForScores = splitCombineAndAlign(seqs.getAlignment().getSequencesArray(), splitIndices, splitSeqs);
 232
 233       for (int[] coords : valuesForScores.keySet())
 234       {
 235         scores[coords[0]][coords[1]] = valuesForScores.get(coords);
 236       }
 237 pairwiseScores = new Matrix(scores);
 238 pairwiseScores.print(System.out, "%1.4f ");
 239
 240 /*
 241       alignment = new PairwiseAlignPanel(seqs, true, 100, 5);
 242       float[][] scores = alignment.getAlignmentScores();        //bigger index first -- eg scores[14][13]
 243
 244       //Hashtable<SequenceI, Integer> connectivity = seqs.calculateConnectivity(scores, dim);
 245
 246       pairwiseScores = new Matrix(scores);
 247 pairwiseScores.print(System.out, "%1.4f ");
 248 /*
 249       pairwiseScores.fillDiagonal();
 250
 251       eigenMatrix = pairwiseScores.copy();
 252
 253       ccAnalysis cc = new ccAnalysis(pairwiseScores, dim);
 254       eigenMatrix = cc.run();
 255 */
 256
 257     } catch (Exception q)
 258     {
 259       Console.error("Error computing PaSiMap:  " + q.getMessage());
 260       q.printStackTrace();
 261     }
 262   }
 263
 264   /**
 265   * aligns sequences in splits
 266   * Splits each split into halves and aligns them agains halves of other splits
 267   *
 268   * @param seqs
 269   * @param i ~ indices of split
 270   * @param s ~ sequences per split
 271   *
 272   * @return  a map of where to put in scores, value ~ scores[n][m] = v
 273   **/
 274   protected HashMap<int[], Float> splitCombineAndAlign(SequenceI[] seqArray, int[] i, int s)
 275   {
 276     HashMap<int[], Float> result = new HashMap<int[], Float>();
 277
 278     int[][] allGroups = new int[i.length][s];
 279     for (int g = 0; g < i.length; g++)  // group g
 280     {
 281       int e = 0;        // index going through allGroups[][e]
 282       for (int j = g * s; j < i[g]; j++)  // goes through all numbers in one group
 283       {
 284         allGroups[g][e++] = j >= seqArray.length ? -1 : j;
 285       }
 286     }
 287
 288     int g = 0;  // group count
 289     for (int[] group : allGroups)
 290     {
 291       HashSet<SequenceI> sg = new HashSet<SequenceI>();
 292       //SequenceGroup sg = new SequenceGroup();
 293       for (int index : group)
 294       {
 295         if (index == -1)
 296           continue;
 297         //sg.addSequence(seqArray[index], false);
 298         sg.add(seqArray[index]);
 299       }
 300       SequenceI[] sgArray = new SequenceI[sg.size()];
 301       int k = 0;
 302       for (SequenceI seq : sg)
 303       {
 304         sgArray[k++] = seq;
 305       }
 306       //seqs.setSelectionGroup(sg);
 307       //PairwiseAlignPanel pap = new PairwiseAlignPanel(seqs, true, 100, 5);
 308       //float[][] scores = pap.getAlignmentScores();    //bigger index first -- eg scores[14][13]
 309       float[][] scores = simulateAlignment(sgArray);
 310       for (int s1 = 0; s1 < scores.length; s1++)        // row
 311       {
 312         result.put(new int[]{s1 + g * s, s1 + g * s}, Float.NaN);       // self score = Float.NaN
 313         for (int s2 = 0; s2 < s1; s2++)                 // col
 314         {
 315           result.put(new int[]{s1 + g * s, s2 + g * s}, scores[s1][s2]);
 316         }
 317       }
 318       g++;
 319     }
 320
 321     int smallS = (int) Math.ceil((float) s/2);
 322     int[][] newGroups = new int[i.length * 2][smallS];
 323
 324     g = 0;
 325     for (int[] group : allGroups)
 326     {
 327       int[] split1 = new int[smallS];
 328       int[] split2 = new int[smallS];
 329       for (int k = 0; k < group.length; k++)
 330       {
 331         if (k < smallS)
 332           split1[k] = group[k];
 333         else
 334           split2[k - smallS] = group[k];
 335       }
 336       newGroups[g++] = split1;
 337       newGroups[g++] = split2;
 338     }
 339
 340     // align each subsplit with subsplits from other split groups
 341     for (int subsplitN = 0; subsplitN < newGroups.length; subsplitN++)
 342     {
 343       int c = 1;                // current split block
 344       while (newGroups[subsplitN][0] > smallS * c)
 345       {
 346         c++;
 347       }
 348       for (int nextSplit = subsplitN + 1; nextSplit < newGroups.length; nextSplit++)
 349       {
 350         if (newGroups[nextSplit][0] >= s * c)   // if next subsplit of next split group -> align seqs
 351         {
 352           HashSet<SequenceI> sg = new HashSet<SequenceI>();
 353           //SequenceGroup sg = new SequenceGroup();
 354           for (int index : newGroups[subsplitN])
 355           {
 356             if (index == -1)
 357               continue;
 358             //sg.addSequence(seqArray[index], false);
 359             sg.add(seqArray[index]);
 360           }
 361           for (int index : newGroups[nextSplit])
 362           {
 363             if (index == -1)
 364               continue;
 365             //sg.addSequence(seqArray[index], false);
 366             sg.add(seqArray[index]);
 367           }
 368           SequenceI[] sgArray = new SequenceI[sg.size()];
 369           int k = 0;
 370           for (SequenceI seq : sg)
 371           {
 372             sgArray[k++] = seq;
 373           }
 374           //seqs.setSelectionGroup(sg);
 375           //PairwiseAlignPanel pap = new PairwiseAlignPanel(seqs, true, 100, 5);
 376           //float[][] scores = pap.getAlignmentScores();        //bigger index first -- eg scores[14][13]
 377           float[][] scores = simulateAlignment(sgArray);
 378           for (int s1 = 0; s1 < scores.length; s1++)    // row
 379           {
 380             for (int s2 = 0; s2 < s1; s2++)                     // col
 381             {
 382               if (s1 >= smallS && s2 < smallS)
 383                 result.put(new int[]{s1 + (nextSplit-1) * smallS, s2 + subsplitN * smallS}, scores[s1][s2]);
 384             }
 385           }
 386         }
 387       }
 388     }
 389
 390     return result;
 391   }
 392
 393   /**
 394   * simulate the alignment of a PairwiseAlignPanel
 395   *
 396   * @param seqs
 397   * @return alignment scores
 398   */
 399   protected float[][] simulateAlignment(SequenceI[] seqs)
 400   {
 401     float[][] result = new float[seqs.length][seqs.length];
 402     for (int i = 1; i < seqs.length; i++)
 403     {
 404       for (int j = 0; j < i; j++)
 405       {
 406         String[] seqStrings = new String[2];
 407         seqStrings[0] = seqs[i].getSequenceAsString();
 408         seqStrings[1] = seqs[j].getSequenceAsString();
 409
 410         AlignSeq as = new AlignSeq(seqs[i], seqStrings[0], seqs[j], seqStrings[1], AlignSeq.PEP);
 411         as.calcScoreMatrix();
 412         as.traceAlignmentWithEndGaps();
 413         as.scoreAlignment();
 414         as.printAlignment(System.out);
 415         result[i][j] = as.getAlignmentScore();
 416       }
 417     }
 418     return result;
 419   }
 420
 421   /**
 422    * Returns a PrintStream that wraps (appends its output to) the given
 423    * StringBuilder
 424    *
 425    * @param sb
 426    * @return
 427    */
 428   protected PrintStream wrapOutputBuffer(StringBuilder sb)
 429   {
 430     PrintStream ps = new PrintStream(System.out)
 431     {
 432       @Override
 433       public void print(String x)
 434       {
 435         sb.append(x);
 436       }
 437
 438       @Override
 439       public void println()
 440       {
 441         sb.append("\n");
 442       }
 443     };
 444     return ps;
 445   }
 446
 447   /**
 448    * Answers the N dimensions of the NxM PaSiMap matrix. This is the number of
 449    * sequences involved in the pairwise score calculation.
 450    *
 451    * @return
 452    */
 453   public int getHeight()
 454   {
 455     // TODO can any of seqs[] be null?
 456     return eigenMatrix.height();// seqs.getSequences().length;
 457   }
 458
 459   /**
 460    * Answers the M dimensions of the NxM PaSiMap matrix. This is the number of
 461    * sequences involved in the pairwise score calculation.
 462    *
 463    * @return
 464    */
 465   public int getWidth()
 466   {
 467     // TODO can any of seqs[] be null?
 468     return eigenMatrix.width();// seqs.getSequences().length;
 469   }
 470
 471   /**
 472    * Answers the sequence pairwise similarity scores which were the first step
 473    * of the PaSiMap calculation
 474    *
 475    * @return
 476    */
 477   public MatrixI getPairwiseScores()
 478   {
 479     return pairwiseScores;
 480   }
 481
 482   public void setPairwiseScores(MatrixI m)
 483   {
 484     pairwiseScores = m;
 485   }
 486
 487   public MatrixI getEigenmatrix()
 488   {
 489     return eigenMatrix;
 490   }
 491
 492   public void setEigenmatrix(MatrixI m)
 493   {
 494     eigenMatrix = m;
 495   }
 496
 497   public PairwiseAlignPanel getAlignments()
 498   {
 499     return alignment;
 500   }
 501
 502   public String getAlignmentOutput()
 503   {
 504     return alignment.getAlignmentOutput();
 505   }
 506
 507   public byte getDim()
 508   {
 509     return dim;
 510   }
 511 }