src/jalview/analysis/AlignSeq.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 java.util.Locale;
  24
  25 import jalview.analysis.scoremodels.PIDModel;
  26 import jalview.analysis.scoremodels.ScoreMatrix;
  27 import jalview.analysis.scoremodels.ScoreModels;
  28 import jalview.analysis.scoremodels.SimilarityParams;
  29 import jalview.api.analysis.SimilarityParamsI;
  30 import jalview.datamodel.AlignmentAnnotation;
  31 import jalview.datamodel.AlignmentI;
  32 import jalview.datamodel.Mapping;
  33 import jalview.datamodel.Sequence;
  34 import jalview.datamodel.SequenceI;
  35 import jalview.util.Comparison;
  36 import jalview.util.Format;
  37 import jalview.util.MapList;
  38 import jalview.util.MessageManager;
  39
  40 import java.awt.Color;
  41 import java.awt.Graphics;
  42 import java.io.PrintStream;
  43 import java.util.ArrayList;
  44 import java.util.Arrays;
  45 import java.util.List;
  46 import java.util.StringTokenizer;
  47
  48 /**
  49  *
  50  *
  51  * @author $author$
  52  * @version $Revision$
  53  */
  54 public class AlignSeq
  55 {
  56   private static final int MAX_NAME_LENGTH = 30;
  57
  58   private static final int GAP_OPEN_COST = 120;
  59
  60   private static final int GAP_EXTEND_COST = 20;
  61
  62   private static final int GAP_INDEX = -1;
  63
  64   public static final String PEP = "pep";
  65
  66   public static final String DNA = "dna";
  67
  68   private static final String NEWLINE = System.lineSeparator();
  69
  70   float[][] score;
  71
  72   float[][] E;
  73
  74   float[][] F;
  75
  76   int[][] traceback; // todo is this actually used?
  77
  78   int[] seq1;
  79
  80   int[] seq2;
  81
  82   SequenceI s1;
  83
  84   SequenceI s2;
  85
  86   public String s1str;
  87
  88   public String s2str;
  89
  90   int maxi;
  91
  92   int maxj;
  93
  94   int[] aseq1;
  95
  96   int[] aseq2;
  97
  98   /*
  99    * matches in alignment
 100    */
 101   int match=-1;
 102
 103   public String astr1 = "";
 104
 105   public String astr2 = "";
 106
 107   /** DOCUMENT ME!! */
 108   public int seq1start;
 109
 110   /** DOCUMENT ME!! */
 111   public int seq1end;
 112
 113   /** DOCUMENT ME!! */
 114   public int seq2start;
 115
 116   public int seq2end;
 117
 118   int count;
 119
 120   public float maxscore;
 121
 122   int prev = 0;
 123
 124   StringBuffer output = new StringBuffer();
 125
 126   String type; // AlignSeq.PEP or AlignSeq.DNA
 127
 128   private ScoreMatrix scoreMatrix;
 129
 130   /**
 131    * Creates a new AlignSeq object.
 132    *
 133    * @param s1
 134    *          first sequence for alignment
 135    * @param s2
 136    *          second sequence for alignment
 137    * @param type
 138    *          molecule type, either AlignSeq.PEP or AlignSeq.DNA
 139    */
 140   public AlignSeq(SequenceI s1, SequenceI s2, String type)
 141   {
 142     seqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
 143             type);
 144   }
 145
 146   /**
 147    * Creates a new AlignSeq object.
 148    *
 149    * @param s1
 150    *          DOCUMENT ME!
 151    * @param s2
 152    *          DOCUMENT ME!
 153    * @param type
 154    *          DOCUMENT ME!
 155    */
 156   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
 157           String string2, String type)
 158   {
 159     seqInit(s1, string1.toUpperCase(Locale.ROOT), s2,
 160             string2.toUpperCase(Locale.ROOT), type);
 161   }
 162
 163   /**
 164    * DOCUMENT ME!
 165    *
 166    * @return DOCUMENT ME!
 167    */
 168   public float getMaxScore()
 169   {
 170     return maxscore;
 171   }
 172
 173   /**
 174    * DOCUMENT ME!
 175    *
 176    * @return DOCUMENT ME!
 177    */
 178   public int getSeq2Start()
 179   {
 180     return seq2start;
 181   }
 182
 183   /**
 184    * DOCUMENT ME!
 185    *
 186    * @return DOCUMENT ME!
 187    */
 188   public int getSeq2End()
 189   {
 190     return seq2end;
 191   }
 192
 193   /**
 194    * DOCUMENT ME!
 195    *
 196    * @return DOCUMENT ME!
 197    */
 198   public int getSeq1Start()
 199   {
 200     return seq1start;
 201   }
 202
 203   /**
 204    * DOCUMENT ME!
 205    *
 206    * @return DOCUMENT ME!
 207    */
 208   public int getSeq1End()
 209   {
 210     return seq1end;
 211   }
 212
 213   /**
 214    * DOCUMENT ME!
 215    *
 216    * @return DOCUMENT ME!
 217    */
 218   public String getOutput()
 219   {
 220     return output.toString();
 221   }
 222
 223   /**
 224    * DOCUMENT ME!
 225    *
 226    * @return DOCUMENT ME!
 227    */
 228   public String getAStr1()
 229   {
 230     return astr1;
 231   }
 232
 233   /**
 234    * DOCUMENT ME!
 235    *
 236    * @return DOCUMENT ME!
 237    */
 238   public String getAStr2()
 239   {
 240     return astr2;
 241   }
 242
 243   /**
 244    * DOCUMENT ME!
 245    *
 246    * @return DOCUMENT ME!
 247    */
 248   public int[] getASeq1()
 249   {
 250     return aseq1;
 251   }
 252
 253   /**
 254    * DOCUMENT ME!
 255    *
 256    * @return DOCUMENT ME!
 257    */
 258   public int[] getASeq2()
 259   {
 260     return aseq2;
 261   }
 262
 263   /**
 264    *
 265    * @return aligned instance of Seq 1
 266    */
 267   public SequenceI getAlignedSeq1()
 268   {
 269     SequenceI alSeq1 = new Sequence(s1.getName(), getAStr1());
 270     alSeq1.setStart(s1.getStart() + getSeq1Start() - 1);
 271     alSeq1.setEnd(s1.getStart() + getSeq1End() - 1);
 272     alSeq1.setDatasetSequence(
 273             s1.getDatasetSequence() == null ? s1 : s1.getDatasetSequence());
 274     return alSeq1;
 275   }
 276
 277   /**
 278    *
 279    * @return aligned instance of Seq 2
 280    */
 281   public SequenceI getAlignedSeq2()
 282   {
 283     SequenceI alSeq2 = new Sequence(s2.getName(), getAStr2());
 284     alSeq2.setStart(s2.getStart() + getSeq2Start() - 1);
 285     alSeq2.setEnd(s2.getStart() + getSeq2End() - 1);
 286     alSeq2.setDatasetSequence(
 287             s2.getDatasetSequence() == null ? s2 : s2.getDatasetSequence());
 288     return alSeq2;
 289   }
 290   /**
 291    * fraction of seq2 matched in the alignment
 292    * @return NaN or [0..1]
 293    */
 294   public double getS2Coverage()
 295   {
 296     if (match>=0)
 297     {
 298       return ((double)match)/((double)s2.getEnd()-s2.getStart()+1);
 299     }
 300     return Double.NaN;
 301   }
 302   /**
 303    * fraction of seq1 matched in the alignment
 304    * @return NaN or [0..1]
 305    */
 306   public double getS1Coverage()
 307   {
 308     if (match>=0)
 309     {
 310       return ((double)match)/((double)s1.getEnd()-s1.getStart()+1);
 311     }
 312     return Double.NaN;
 313   }
 314
 315   /**
 316    * Construct score matrix for sequences with standard DNA or PEPTIDE matrix
 317    *
 318    * @param s1
 319    *          - sequence 1
 320    * @param string1
 321    *          - string to use for s1
 322    * @param s2
 323    *          - sequence 2
 324    * @param string2
 325    *          - string to use for s2
 326    * @param type
 327    *          DNA or PEPTIDE
 328    */
 329   public void seqInit(SequenceI s1, String string1, SequenceI s2,
 330           String string2, String type)
 331   {
 332     this.s1 = s1;
 333     this.s2 = s2;
 334     setDefaultParams(type);
 335     seqInit(string1, string2);
 336   }
 337
 338   /**
 339    * construct score matrix for string1 and string2 (after removing any existing
 340    * gaps
 341    *
 342    * @param string1
 343    * @param string2
 344    */
 345   private void seqInit(String string1, String string2)
 346   {
 347     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
 348     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
 349
 350     if (s1str.length() == 0 || s2str.length() == 0)
 351     {
 352       output.append(
 353               "ALL GAPS: " + (s1str.length() == 0 ? s1.getName() : " ")
 354                       + (s2str.length() == 0 ? s2.getName() : ""));
 355       return;
 356     }
 357
 358     score = new float[s1str.length()][s2str.length()];
 359
 360     E = new float[s1str.length()][s2str.length()];
 361
 362     F = new float[s1str.length()][s2str.length()];
 363     traceback = new int[s1str.length()][s2str.length()];
 364
 365     seq1 = indexEncode(s1str);
 366
 367     seq2 = indexEncode(s2str);
 368   }
 369
 370   private void setDefaultParams(String moleculeType)
 371   {
 372     if (!PEP.equals(moleculeType) && !DNA.equals(moleculeType))
 373     {
 374       output.append("Wrong type = dna or pep only");
 375       throw new Error(MessageManager
 376               .formatMessage("error.unknown_type_dna_or_pep", new String[]
 377               { moleculeType }));
 378     }
 379
 380     type = moleculeType;
 381     scoreMatrix = ScoreModels.getInstance()
 382             .getDefaultModel(PEP.equals(type));
 383   }
 384
 385   /**
 386    * DOCUMENT ME!
 387    */
 388   public void traceAlignment()
 389   {
 390     // Find the maximum score along the rhs or bottom row
 391     float max = -Float.MAX_VALUE;
 392
 393     for (int i = 0; i < seq1.length; i++)
 394     {
 395       if (score[i][seq2.length - 1] > max)
 396       {
 397         max = score[i][seq2.length - 1];
 398         maxi = i;
 399         maxj = seq2.length - 1;
 400       }
 401     }
 402
 403     for (int j = 0; j < seq2.length; j++)
 404     {
 405       if (score[seq1.length - 1][j] > max)
 406       {
 407         max = score[seq1.length - 1][j];
 408         maxi = seq1.length - 1;
 409         maxj = j;
 410       }
 411     }
 412
 413     int i = maxi;
 414     int j = maxj;
 415     int trace;
 416     maxscore = score[i][j] / 10f;
 417
 418     seq1end = maxi + 1;
 419     seq2end = maxj + 1;
 420
 421     aseq1 = new int[seq1.length + seq2.length];
 422     aseq2 = new int[seq1.length + seq2.length];
 423     match=0;
 424     StringBuilder sb1 = new StringBuilder(aseq1.length);
 425     StringBuilder sb2 = new StringBuilder(aseq2.length);
 426
 427     count = (seq1.length + seq2.length) - 1;
 428
 429     while (i > 0 && j > 0)
 430     {
 431       aseq1[count] = seq1[i];
 432       sb1.append(s1str.charAt(i));
 433       aseq2[count] = seq2[j];
 434       sb2.append(s2str.charAt(j));
 435       trace = findTrace(i, j);
 436
 437       if (trace == 0)
 438       {
 439         match++;
 440         i--;
 441         j--;
 442       }
 443       else if (trace == 1)
 444       {
 445         j--;
 446         aseq1[count] = GAP_INDEX;
 447         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 448       }
 449       else if (trace == -1)
 450       {
 451         i--;
 452         aseq2[count] = GAP_INDEX;
 453         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 454       }
 455
 456       count--;
 457     }
 458
 459     seq1start = i + 1;
 460     seq2start = j + 1;
 461
 462     if (aseq1[count] != GAP_INDEX)
 463     {
 464       aseq1[count] = seq1[i];
 465       sb1.append(s1str.charAt(i));
 466     }
 467
 468     if (aseq2[count] != GAP_INDEX)
 469     {
 470       aseq2[count] = seq2[j];
 471       sb2.append(s2str.charAt(j));
 472       if (aseq1[count]!=GAP_INDEX) {
 473         match++;
 474       }
 475     }
 476     /*
 477      * we built the character strings backwards, so now
 478      * reverse them to convert to sequence strings
 479      */
 480     astr1 = sb1.reverse().toString();
 481     astr2 = sb2.reverse().toString();
 482   }
 483
 484   /**
 485    * DOCUMENT ME!
 486    */
 487   public void printAlignment(PrintStream os)
 488   {
 489     // TODO: Use original sequence characters rather than re-translated
 490     // characters in output
 491     // Find the biggest id length for formatting purposes
 492     String s1id = getAlignedSeq1().getDisplayId(true);
 493     String s2id = getAlignedSeq2().getDisplayId(true);
 494     int nameLength = Math.max(s1id.length(), s2id.length());
 495     if (nameLength > MAX_NAME_LENGTH)
 496     {
 497       int truncateBy = nameLength - MAX_NAME_LENGTH;
 498       nameLength = MAX_NAME_LENGTH;
 499       // JAL-527 - truncate the sequence ids
 500       if (s1id.length() > nameLength)
 501       {
 502         int slashPos = s1id.lastIndexOf('/');
 503         s1id = s1id.substring(0, slashPos - truncateBy)
 504                 + s1id.substring(slashPos);
 505       }
 506       if (s2id.length() > nameLength)
 507       {
 508         int slashPos = s2id.lastIndexOf('/');
 509         s2id = s2id.substring(0, slashPos - truncateBy)
 510                 + s2id.substring(slashPos);
 511       }
 512     }
 513     int len = 72 - nameLength - 1;
 514     int nochunks = ((aseq1.length - count) / len)
 515             + ((aseq1.length - count) % len > 0 ? 1 : 0);
 516     float pid = 0f;
 517
 518     output.append("Score = ").append(score[maxi][maxj]).append(NEWLINE);
 519     output.append("Length of alignment = ")
 520             .append(String.valueOf(aseq1.length - count)).append(NEWLINE);
 521     output.append("Sequence ");
 522     Format nameFormat = new Format("%" + nameLength + "s");
 523     output.append(nameFormat.form(s1id));
 524     output.append(" (Sequence length = ")
 525             .append(String.valueOf(s1str.length())).append(")")
 526             .append(NEWLINE);
 527     output.append("Sequence ");
 528     output.append(nameFormat.form(s2id));
 529     output.append(" (Sequence length = ")
 530             .append(String.valueOf(s2str.length())).append(")")
 531             .append(NEWLINE).append(NEWLINE);
 532
 533     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
 534
 535     for (int j = 0; j < nochunks; j++)
 536     {
 537       // Print the first aligned sequence
 538       output.append(nameFormat.form(s1id)).append(" ");
 539
 540       for (int i = 0; i < len; i++)
 541       {
 542         if ((i + (j * len)) < astr1.length())
 543         {
 544           output.append(astr1.charAt(i + (j * len)));
 545         }
 546       }
 547
 548       output.append(NEWLINE);
 549       output.append(nameFormat.form(" ")).append(" ");
 550
 551       /*
 552        * Print out the match symbols:
 553        * | for exact match (ignoring case)
 554        * . if PAM250 score is positive
 555        * else a space
 556        */
 557       for (int i = 0; i < len; i++)
 558       {
 559         if ((i + (j * len)) < astr1.length())
 560         {
 561           char c1 = astr1.charAt(i + (j * len));
 562           char c2 = astr2.charAt(i + (j * len));
 563           boolean sameChar = Comparison.isSameResidue(c1, c2, false);
 564           if (sameChar && !Comparison.isGap(c1))
 565           {
 566             pid++;
 567             output.append("|");
 568           }
 569           else if (PEP.equals(type))
 570           {
 571             if (pam250.getPairwiseScore(c1, c2) > 0)
 572             {
 573               output.append(".");
 574             }
 575             else
 576             {
 577               output.append(" ");
 578             }
 579           }
 580           else
 581           {
 582             output.append(" ");
 583           }
 584         }
 585       }
 586
 587       // Now print the second aligned sequence
 588       output = output.append(NEWLINE);
 589       output = output.append(nameFormat.form(s2id)).append(" ");
 590
 591       for (int i = 0; i < len; i++)
 592       {
 593         if ((i + (j * len)) < astr2.length())
 594         {
 595           output.append(astr2.charAt(i + (j * len)));
 596         }
 597       }
 598
 599       output.append(NEWLINE).append(NEWLINE);
 600     }
 601
 602     pid = pid / (aseq1.length - count) * 100;
 603     output.append(new Format("Percentage ID = %3.2f\n").form(pid));
 604     output.append(NEWLINE);
 605     try
 606     {
 607       os.print(output.toString());
 608     } catch (Exception ex)
 609     {
 610     }
 611   }
 612
 613   /**
 614    * DOCUMENT ME!
 615    *
 616    * @param i
 617    *          DOCUMENT ME!
 618    * @param j
 619    *          DOCUMENT ME!
 620    *
 621    * @return DOCUMENT ME!
 622    */
 623   public int findTrace(int i, int j)
 624   {
 625     int t = 0;
 626     float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 627             s2str.charAt(j));
 628     float max = score[i - 1][j - 1] + (pairwiseScore * 10);
 629
 630     if (F[i][j] > max)
 631     {
 632       max = F[i][j];
 633       t = -1;
 634     }
 635     else if (F[i][j] == max)
 636     {
 637       if (prev == -1)
 638       {
 639         max = F[i][j];
 640         t = -1;
 641       }
 642     }
 643
 644     if (E[i][j] >= max)
 645     {
 646       max = E[i][j];
 647       t = 1;
 648     }
 649     else if (E[i][j] == max)
 650     {
 651       if (prev == 1)
 652       {
 653         max = E[i][j];
 654         t = 1;
 655       }
 656     }
 657
 658     prev = t;
 659
 660     return t;
 661   }
 662
 663   /**
 664    * DOCUMENT ME!
 665    */
 666   public void calcScoreMatrix()
 667   {
 668     int n = seq1.length;
 669     int m = seq2.length;
 670
 671     // top left hand element
 672     score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 673             s2str.charAt(0)) * 10;
 674     E[0][0] = -GAP_EXTEND_COST;
 675     F[0][0] = 0;
 676
 677     // Calculate the top row first
 678     for (int j = 1; j < m; j++)
 679     {
 680       // What should these values be? 0 maybe
 681       E[0][j] = max(score[0][j - 1] - GAP_OPEN_COST,
 682               E[0][j - 1] - GAP_EXTEND_COST);
 683       F[0][j] = -GAP_EXTEND_COST;
 684
 685       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 686               s2str.charAt(j));
 687       score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
 688               -GAP_EXTEND_COST);
 689
 690       traceback[0][j] = 1;
 691     }
 692
 693     // Now do the left hand column
 694     for (int i = 1; i < n; i++)
 695     {
 696       E[i][0] = -GAP_OPEN_COST;
 697       F[i][0] = max(score[i - 1][0] - GAP_OPEN_COST,
 698               F[i - 1][0] - GAP_EXTEND_COST);
 699
 700       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 701               s2str.charAt(0));
 702       score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
 703       traceback[i][0] = -1;
 704     }
 705
 706     // Now do all the other rows
 707     for (int i = 1; i < n; i++)
 708     {
 709       for (int j = 1; j < m; j++)
 710       {
 711         E[i][j] = max(score[i][j - 1] - GAP_OPEN_COST,
 712                 E[i][j - 1] - GAP_EXTEND_COST);
 713         F[i][j] = max(score[i - 1][j] - GAP_OPEN_COST,
 714                 F[i - 1][j] - GAP_EXTEND_COST);
 715
 716         float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 717                 s2str.charAt(j));
 718         score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
 719                 E[i][j], F[i][j]);
 720         traceback[i][j] = findTrace(i, j);
 721       }
 722     }
 723   }
 724
 725   /**
 726    * Returns the given sequence with all of the given gap characters removed.
 727    *
 728    * @param gapChars
 729    *          a string of characters to be treated as gaps
 730    * @param seq
 731    *          the input sequence
 732    *
 733    * @return
 734    */
 735   public static String extractGaps(String gapChars, String seq)
 736   {
 737     if (gapChars == null || seq == null)
 738     {
 739       return null;
 740     }
 741     StringTokenizer str = new StringTokenizer(seq, gapChars);
 742     StringBuilder newString = new StringBuilder(seq.length());
 743
 744     while (str.hasMoreTokens())
 745     {
 746       newString.append(str.nextToken());
 747     }
 748
 749     return newString.toString();
 750   }
 751
 752   /**
 753    * DOCUMENT ME!
 754    *
 755    * @param f1
 756    *          DOCUMENT ME!
 757    * @param f2
 758    *          DOCUMENT ME!
 759    * @param f3
 760    *          DOCUMENT ME!
 761    *
 762    * @return DOCUMENT ME!
 763    */
 764   private static float max(float f1, float f2, float f3)
 765   {
 766     float max = f1;
 767
 768     if (f2 > f1)
 769     {
 770       max = f2;
 771     }
 772
 773     if (f3 > max)
 774     {
 775       max = f3;
 776     }
 777
 778     return max;
 779   }
 780
 781   /**
 782    * DOCUMENT ME!
 783    *
 784    * @param f1
 785    *          DOCUMENT ME!
 786    * @param f2
 787    *          DOCUMENT ME!
 788    *
 789    * @return DOCUMENT ME!
 790    */
 791   private static float max(float f1, float f2)
 792   {
 793     float max = f1;
 794
 795     if (f2 > f1)
 796     {
 797       max = f2;
 798     }
 799
 800     return max;
 801   }
 802
 803   /**
 804    * Converts the character string to an array of integers which are the
 805    * corresponding indices to the characters in the score matrix
 806    *
 807    * @param s
 808    *
 809    * @return
 810    */
 811   int[] indexEncode(String s)
 812   {
 813     int[] encoded = new int[s.length()];
 814
 815     for (int i = 0; i < s.length(); i++)
 816     {
 817       char c = s.charAt(i);
 818       encoded[i] = scoreMatrix.getMatrixIndex(c);
 819     }
 820
 821     return encoded;
 822   }
 823
 824   /**
 825    * DOCUMENT ME!
 826    *
 827    * @param g
 828    *          DOCUMENT ME!
 829    * @param mat
 830    *          DOCUMENT ME!
 831    * @param n
 832    *          DOCUMENT ME!
 833    * @param m
 834    *          DOCUMENT ME!
 835    * @param psize
 836    *          DOCUMENT ME!
 837    */
 838   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 839           int psize)
 840   {
 841     // TODO method doesn't seem to be referenced anywhere delete??
 842     int max = -1000;
 843     int min = 1000;
 844
 845     for (int i = 0; i < n; i++)
 846     {
 847       for (int j = 0; j < m; j++)
 848       {
 849         if (mat[i][j] >= max)
 850         {
 851           max = mat[i][j];
 852         }
 853
 854         if (mat[i][j] <= min)
 855         {
 856           min = mat[i][j];
 857         }
 858       }
 859     }
 860
 861     jalview.bin.Console.outPrintln(max + " " + min);
 862
 863     for (int i = 0; i < n; i++)
 864     {
 865       for (int j = 0; j < m; j++)
 866       {
 867         int x = psize * i;
 868         int y = psize * j;
 869
 870         // jalview.bin.Console.outPrintln(mat[i][j]);
 871         float score = (float) (mat[i][j] - min) / (float) (max - min);
 872         g.setColor(new Color(score, 0, 0));
 873         g.fillRect(x, y, psize, psize);
 874
 875         // jalview.bin.Console.outPrintln(x + " " + y + " " + score);
 876       }
 877     }
 878   }
 879
 880   /**
 881    * Compute a globally optimal needleman and wunsch alignment between two
 882    * sequences
 883    *
 884    * @param s1
 885    * @param s2
 886    * @param type
 887    *          AlignSeq.DNA or AlignSeq.PEP
 888    */
 889   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
 890           String type)
 891   {
 892     AlignSeq as = new AlignSeq(s1, s2, type);
 893
 894     as.calcScoreMatrix();
 895     as.traceAlignment();
 896     return as;
 897   }
 898
 899   /**
 900    *
 901    * @return mapping from positions in S1 to corresponding positions in S2
 902    */
 903   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
 904   {
 905     ArrayList<Integer> as1 = new ArrayList<Integer>(),
 906             as2 = new ArrayList<Integer>();
 907     int pdbpos = s2.getStart() + getSeq2Start() - 2;
 908     int alignpos = s1.getStart() + getSeq1Start() - 2;
 909     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
 910     boolean lastmatch = false;
 911     // and now trace the alignment onto the atom set.
 912     for (int i = 0; i < astr1.length(); i++)
 913     {
 914       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
 915       if (c1 != '-')
 916       {
 917         alignpos++;
 918       }
 919
 920       if (c2 != '-')
 921       {
 922         pdbpos++;
 923       }
 924
 925       // ignore case differences
 926       if (allowmismatch || (c1 == c2) || (Math.abs(c2 - c1) == ('a' - 'A')))
 927       {
 928         // extend mapping interval
 929         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
 930         {
 931           as1.add(Integer.valueOf(alignpos));
 932           as2.add(Integer.valueOf(pdbpos));
 933         }
 934         lastmatch = true;
 935         lp1 = alignpos;
 936         lp2 = pdbpos;
 937       }
 938       else
 939       {
 940         // extend mapping interval
 941         if (lastmatch)
 942         {
 943           as1.add(Integer.valueOf(lp1));
 944           as2.add(Integer.valueOf(lp2));
 945         }
 946         lastmatch = false;
 947       }
 948     }
 949     // construct range pairs
 950
 951     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
 952             mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
 953     int i = 0;
 954     for (Integer ip : as1)
 955     {
 956       mapseq1[i++] = ip;
 957     }
 958     ;
 959     i = 0;
 960     for (Integer ip : as2)
 961     {
 962       mapseq2[i++] = ip;
 963     }
 964     ;
 965     if (lastmatch)
 966     {
 967       mapseq1[mapseq1.length - 1] = alignpos;
 968       mapseq2[mapseq2.length - 1] = pdbpos;
 969     }
 970     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
 971
 972     jalview.datamodel.Mapping mapping = new Mapping(map);
 973     mapping.setTo(s2);
 974     return mapping;
 975   }
 976
 977   /**
 978    * matches ochains against al and populates seqs with the best match between
 979    * each ochain and the set in al
 980    *
 981    * @param ochains
 982    * @param al
 983    * @param dnaOrProtein
 984    * @param removeOldAnnots
 985    *          when true, old annotation is cleared before new annotation
 986    *          transferred
 987    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
 988    *         List<AlignSeq> alignment between each>
 989    */
 990   public static List<List<? extends Object>> replaceMatchingSeqsWith(
 991           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
 992           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
 993           boolean removeOldAnnots)
 994   {
 995     List<SequenceI> orig = new ArrayList<SequenceI>(),
 996             repl = new ArrayList<SequenceI>();
 997     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
 998     if (al != null && al.getHeight() > 0)
 999     {
1000       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
1001       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
1002
1003       for (SequenceI sq : ochains)
1004       {
1005         SequenceI bestm = null;
1006         AlignSeq bestaseq = null;
1007         float bestscore = 0;
1008         for (SequenceI msq : al.getSequences())
1009         {
1010           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
1011           if (bestm == null || aseq.getMaxScore() > bestscore)
1012           {
1013             bestscore = aseq.getMaxScore();
1014             bestaseq = aseq;
1015             bestm = msq;
1016           }
1017         }
1018         // jalview.bin.Console.outPrintln("Best Score for " + (matches.size() +
1019         // 1) + " :"
1020         // + bestscore);
1021         matches.add(bestm);
1022         aligns.add(bestaseq);
1023         al.deleteSequence(bestm);
1024       }
1025       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
1026       {
1027         SequenceI sq, sp = seqs.get(p);
1028         int q;
1029         if ((q = ochains.indexOf(sp)) > -1)
1030         {
1031           seqs.set(p, sq = matches.get(q));
1032           orig.add(sp);
1033           repl.add(sq);
1034           sq.setName(sp.getName());
1035           sq.setDescription(sp.getDescription());
1036           Mapping sp2sq;
1037           sq.transferAnnotation(sp,
1038                   sp2sq = aligns.get(q).getMappingFromS1(false));
1039           aligs.add(aligns.get(q));
1040           int inspos = -1;
1041           for (int ap = 0; ap < annotations.size();)
1042           {
1043             if (annotations.get(ap).sequenceRef == sp)
1044             {
1045               if (inspos == -1)
1046               {
1047                 inspos = ap;
1048               }
1049               if (removeOldAnnots)
1050               {
1051                 annotations.remove(ap);
1052               }
1053               else
1054               {
1055                 AlignmentAnnotation alan = annotations.remove(ap);
1056                 alan.liftOver(sq, sp2sq);
1057                 alan.setSequenceRef(sq);
1058                 sq.addAlignmentAnnotation(alan);
1059               }
1060             }
1061             else
1062             {
1063               ap++;
1064             }
1065           }
1066           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1067           {
1068             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1069                     Arrays.asList(sq.getAnnotation()));
1070           }
1071         }
1072       }
1073     }
1074     return Arrays.asList(orig, repl, aligs);
1075   }
1076
1077   /**
1078    * compute the PID vector used by the redundancy filter.
1079    *
1080    * @param originalSequences
1081    *          - sequences in alignment that are to filtered
1082    * @param omitHidden
1083    *          - null or strings to be analysed (typically, visible portion of
1084    *          each sequence in alignment)
1085    * @param start
1086    *          - first column in window for calculation
1087    * @param end
1088    *          - last column in window for calculation
1089    * @param ungapped
1090    *          - if true then use ungapped sequence to compute PID
1091    * @return vector containing maximum PID for i-th sequence and any sequences
1092    *         longer than that seuqence
1093    */
1094   public static float[] computeRedundancyMatrix(
1095           SequenceI[] originalSequences, String[] omitHidden, int start,
1096           int end, boolean ungapped)
1097   {
1098     int height = originalSequences.length;
1099     float[] redundancy = new float[height];
1100     int[] lngth = new int[height];
1101     for (int i = 0; i < height; i++)
1102     {
1103       redundancy[i] = 0f;
1104       lngth[i] = -1;
1105     }
1106
1107     // long start = System.currentTimeMillis();
1108
1109     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1110             true);
1111     float pid;
1112     String seqi, seqj;
1113     for (int i = 0; i < height; i++)
1114     {
1115
1116       for (int j = 0; j < i; j++)
1117       {
1118         if (i == j)
1119         {
1120           continue;
1121         }
1122
1123         if (omitHidden == null)
1124         {
1125           seqi = originalSequences[i].getSequenceAsString(start, end);
1126           seqj = originalSequences[j].getSequenceAsString(start, end);
1127         }
1128         else
1129         {
1130           seqi = omitHidden[i];
1131           seqj = omitHidden[j];
1132         }
1133         if (lngth[i] == -1)
1134         {
1135           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1136           lngth[i] = ug.length();
1137           if (ungapped)
1138           {
1139             seqi = ug;
1140           }
1141         }
1142         if (lngth[j] == -1)
1143         {
1144           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1145           lngth[j] = ug.length();
1146           if (ungapped)
1147           {
1148             seqj = ug;
1149           }
1150         }
1151         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1152
1153         // use real sequence length rather than string length
1154         if (lngth[j] < lngth[i])
1155         {
1156           redundancy[j] = Math.max(pid, redundancy[j]);
1157         }
1158         else
1159         {
1160           redundancy[i] = Math.max(pid, redundancy[i]);
1161         }
1162
1163       }
1164     }
1165     return redundancy;
1166   }
1167
1168   public void setScoreMatrix(ScoreMatrix params)
1169   {
1170     if (params!=null)
1171     {
1172       scoreMatrix = params;
1173     }
1174   }
1175 }