src/jalview/analysis/AlignSeq.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.0b1)
   3  * Copyright (C) 2014 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 of the License, or (at your option) any later version.
  10  *
  11  * Jalview is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14  * PURPOSE.  See the GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  17  * The Jalview Authors are detailed in the 'AUTHORS' file.
  18  */
  19 package jalview.analysis;
  20
  21 import java.util.*;
  22
  23 import java.awt.*;
  24
  25 import jalview.datamodel.*;
  26 import jalview.schemes.*;
  27 import jalview.util.*;
  28
  29 /**
  30  *
  31  *
  32  * @author $author$
  33  * @version $Revision$
  34  */
  35 public class AlignSeq
  36 {
  37   public static final String PEP = "pep";
  38
  39   public static final String DNA = "dna";
  40
  41   static String[] dna =
  42   { "A", "C", "G", "T", "-" };
  43
  44   // "C", "T", "A", "G", "-"};
  45   static String[] pep =
  46   { "A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F",
  47       "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "-" };
  48
  49   int[][] score;
  50
  51   int[][] E;
  52
  53   int[][] F;
  54
  55   int[][] traceback;
  56
  57   int[] seq1;
  58
  59   int[] seq2;
  60
  61   SequenceI s1;
  62
  63   SequenceI s2;
  64
  65   public String s1str;
  66
  67   public String s2str;
  68
  69   int maxi;
  70
  71   int maxj;
  72
  73   int[] aseq1;
  74
  75   int[] aseq2;
  76
  77   public String astr1 = "";
  78
  79   public String astr2 = "";
  80
  81   /** DOCUMENT ME!! */
  82   public int seq1start;
  83
  84   /** DOCUMENT ME!! */
  85   public int seq1end;
  86
  87   /** DOCUMENT ME!! */
  88   public int seq2start;
  89
  90   /** DOCUMENT ME!! */
  91   public int seq2end;
  92
  93   int count;
  94
  95   /** DOCUMENT ME!! */
  96   public int maxscore;
  97
  98   float pid;
  99
 100   int prev = 0;
 101
 102   int gapOpen = 120;
 103
 104   int gapExtend = 20;
 105
 106   int[][] lookup = ResidueProperties.getBLOSUM62();
 107
 108   String[] intToStr = pep;
 109
 110   int defInt = 23;
 111
 112   StringBuffer output = new StringBuffer();
 113
 114   String type;
 115
 116   private int[] charToInt;
 117
 118   /**
 119    * Creates a new AlignSeq object.
 120    *
 121    * @param s1
 122    *          DOCUMENT ME!
 123    * @param s2
 124    *          DOCUMENT ME!
 125    * @param type
 126    *          DOCUMENT ME!
 127    */
 128   public AlignSeq(SequenceI s1, SequenceI s2, String type)
 129   {
 130     SeqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
 131             type);
 132   }
 133
 134   /**
 135    * Creates a new AlignSeq object.
 136    *
 137    * @param s1
 138    *          DOCUMENT ME!
 139    * @param s2
 140    *          DOCUMENT ME!
 141    * @param type
 142    *          DOCUMENT ME!
 143    */
 144   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
 145           String string2, String type)
 146   {
 147     SeqInit(s1, string1.toUpperCase(), s2, string2.toUpperCase(), type);
 148   }
 149
 150   /**
 151    * DOCUMENT ME!
 152    *
 153    * @return DOCUMENT ME!
 154    */
 155   public int getMaxScore()
 156   {
 157     return maxscore;
 158   }
 159
 160   /**
 161    * DOCUMENT ME!
 162    *
 163    * @return DOCUMENT ME!
 164    */
 165   public int getSeq2Start()
 166   {
 167     return seq2start;
 168   }
 169
 170   /**
 171    * DOCUMENT ME!
 172    *
 173    * @return DOCUMENT ME!
 174    */
 175   public int getSeq2End()
 176   {
 177     return seq2end;
 178   }
 179
 180   /**
 181    * DOCUMENT ME!
 182    *
 183    * @return DOCUMENT ME!
 184    */
 185   public int getSeq1Start()
 186   {
 187     return seq1start;
 188   }
 189
 190   /**
 191    * DOCUMENT ME!
 192    *
 193    * @return DOCUMENT ME!
 194    */
 195   public int getSeq1End()
 196   {
 197     return seq1end;
 198   }
 199
 200   /**
 201    * DOCUMENT ME!
 202    *
 203    * @return DOCUMENT ME!
 204    */
 205   public String getOutput()
 206   {
 207     return output.toString();
 208   }
 209
 210   /**
 211    * DOCUMENT ME!
 212    *
 213    * @return DOCUMENT ME!
 214    */
 215   public String getAStr1()
 216   {
 217     return astr1;
 218   }
 219
 220   /**
 221    * DOCUMENT ME!
 222    *
 223    * @return DOCUMENT ME!
 224    */
 225   public String getAStr2()
 226   {
 227     return astr2;
 228   }
 229
 230   /**
 231    * DOCUMENT ME!
 232    *
 233    * @return DOCUMENT ME!
 234    */
 235   public int[] getASeq1()
 236   {
 237     return aseq1;
 238   }
 239
 240   /**
 241    * DOCUMENT ME!
 242    *
 243    * @return DOCUMENT ME!
 244    */
 245   public int[] getASeq2()
 246   {
 247     return aseq2;
 248   }
 249
 250   /**
 251    * DOCUMENT ME!
 252    *
 253    * @return DOCUMENT ME!
 254    */
 255   public SequenceI getS1()
 256   {
 257     return s1;
 258   }
 259
 260   /**
 261    * DOCUMENT ME!
 262    *
 263    * @return DOCUMENT ME!
 264    */
 265   public SequenceI getS2()
 266   {
 267     return s2;
 268   }
 269
 270   /**
 271    * DOCUMENT ME!
 272    *
 273    * @param s1
 274    *          DOCUMENT ME!
 275    * @param string1
 276    *          - string to align for sequence1
 277    * @param s2
 278    *          sequence 2
 279    * @param string2
 280    *          - string to align for sequence2
 281    * @param type
 282    *          DNA or PEPTIDE
 283    */
 284   public void SeqInit(SequenceI s1, String string1, SequenceI s2,
 285           String string2, String type)
 286   {
 287     this.s1 = s1;
 288     this.s2 = s2;
 289     setDefaultParams(type);
 290     SeqInit(string1, string2);
 291   }
 292
 293   public void SeqInit(SequenceI s1, String string1, SequenceI s2,
 294           String string2, ScoreMatrix scoreMatrix)
 295   {
 296     this.s1 = s1;
 297     this.s2 = s2;
 298     setType(scoreMatrix.isDNA() ? AlignSeq.DNA : AlignSeq.PEP);
 299     lookup = scoreMatrix.getMatrix();
 300   }
 301
 302   /**
 303    * construct score matrix for string1 and string2 (after removing any existing
 304    * gaps
 305    *
 306    * @param string1
 307    * @param string2
 308    */
 309   private void SeqInit(String string1, String string2)
 310   {
 311     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
 312     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
 313
 314     if (s1str.length() == 0 || s2str.length() == 0)
 315     {
 316       output.append("ALL GAPS: "
 317               + (s1str.length() == 0 ? s1.getName() : " ")
 318               + (s2str.length() == 0 ? s2.getName() : ""));
 319       return;
 320     }
 321
 322     // System.out.println("lookuip " + rt.freeMemory() + " "+ rt.totalMemory());
 323     seq1 = new int[s1str.length()];
 324
 325     // System.out.println("seq1 " + rt.freeMemory() +" " + rt.totalMemory());
 326     seq2 = new int[s2str.length()];
 327
 328     // System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());
 329     score = new int[s1str.length()][s2str.length()];
 330
 331     // System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());
 332     E = new int[s1str.length()][s2str.length()];
 333
 334     // System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());
 335     F = new int[s1str.length()][s2str.length()];
 336     traceback = new int[s1str.length()][s2str.length()];
 337
 338     // System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());
 339     seq1 = stringToInt(s1str, type);
 340
 341     // System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());
 342     seq2 = stringToInt(s2str, type);
 343
 344     // System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());
 345     // long tstart = System.currentTimeMillis();
 346     // calcScoreMatrix();
 347     // long tend = System.currentTimeMillis();
 348     // System.out.println("Time take to calculate score matrix = " +
 349     // (tend-tstart) + " ms");
 350     // printScoreMatrix(score);
 351     // System.out.println();
 352     // printScoreMatrix(traceback);
 353     // System.out.println();
 354     // printScoreMatrix(E);
 355     // System.out.println();
 356     // /printScoreMatrix(F);
 357     // System.out.println();
 358     // tstart = System.currentTimeMillis();
 359     // traceAlignment();
 360     // tend = System.currentTimeMillis();
 361     // System.out.println("Time take to traceback alignment = " + (tend-tstart)
 362     // + " ms");
 363   }
 364
 365   private void setDefaultParams(String type)
 366   {
 367     setType(type);
 368
 369     if (type.equals(AlignSeq.PEP))
 370     {
 371       lookup = ResidueProperties.getDefaultPeptideMatrix();
 372     }
 373     else if (type.equals(AlignSeq.DNA))
 374     {
 375       lookup = ResidueProperties.getDefaultDnaMatrix();
 376     }
 377   }
 378
 379   private void setType(String type2)
 380   {
 381     this.type = type2;
 382     if (type.equals(AlignSeq.PEP))
 383     {
 384       intToStr = pep;
 385       charToInt = ResidueProperties.aaIndex;
 386       defInt = ResidueProperties.maxProteinIndex;
 387     }
 388     else if (type.equals(AlignSeq.DNA))
 389     {
 390       intToStr = dna;
 391       charToInt = ResidueProperties.nucleotideIndex;
 392       defInt = ResidueProperties.maxNucleotideIndex;
 393     }
 394     else
 395     {
 396       output.append("Wrong type = dna or pep only");
 397       throw new Error("Unknown Type " + type2
 398               + " - dna or pep are the only allowed values.");
 399     }
 400   }
 401
 402   /**
 403    * DOCUMENT ME!
 404    */
 405   public void traceAlignment()
 406   {
 407     // Find the maximum score along the rhs or bottom row
 408     int max = -9999;
 409
 410     for (int i = 0; i < seq1.length; i++)
 411     {
 412       if (score[i][seq2.length - 1] > max)
 413       {
 414         max = score[i][seq2.length - 1];
 415         maxi = i;
 416         maxj = seq2.length - 1;
 417       }
 418     }
 419
 420     for (int j = 0; j < seq2.length; j++)
 421     {
 422       if (score[seq1.length - 1][j] > max)
 423       {
 424         max = score[seq1.length - 1][j];
 425         maxi = seq1.length - 1;
 426         maxj = j;
 427       }
 428     }
 429
 430     // System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);
 431     int i = maxi;
 432     int j = maxj;
 433     int trace;
 434     maxscore = score[i][j] / 10;
 435
 436     seq1end = maxi + 1;
 437     seq2end = maxj + 1;
 438
 439     aseq1 = new int[seq1.length + seq2.length];
 440     aseq2 = new int[seq1.length + seq2.length];
 441
 442     count = (seq1.length + seq2.length) - 1;
 443
 444     while ((i > 0) && (j > 0))
 445     {
 446       if ((aseq1[count] != defInt) && (i >= 0))
 447       {
 448         aseq1[count] = seq1[i];
 449         astr1 = s1str.charAt(i) + astr1;
 450       }
 451
 452       if ((aseq2[count] != defInt) && (j > 0))
 453       {
 454         aseq2[count] = seq2[j];
 455         astr2 = s2str.charAt(j) + astr2;
 456       }
 457
 458       trace = findTrace(i, j);
 459
 460       if (trace == 0)
 461       {
 462         i--;
 463         j--;
 464       }
 465       else if (trace == 1)
 466       {
 467         j--;
 468         aseq1[count] = defInt;
 469         astr1 = "-" + astr1.substring(1);
 470       }
 471       else if (trace == -1)
 472       {
 473         i--;
 474         aseq2[count] = defInt;
 475         astr2 = "-" + astr2.substring(1);
 476       }
 477
 478       count--;
 479     }
 480
 481     seq1start = i + 1;
 482     seq2start = j + 1;
 483
 484     if (aseq1[count] != defInt)
 485     {
 486       aseq1[count] = seq1[i];
 487       astr1 = s1str.charAt(i) + astr1;
 488     }
 489
 490     if (aseq2[count] != defInt)
 491     {
 492       aseq2[count] = seq2[j];
 493       astr2 = s2str.charAt(j) + astr2;
 494     }
 495   }
 496
 497   /**
 498    * DOCUMENT ME!
 499    */
 500   public void printAlignment(java.io.PrintStream os)
 501   {
 502     // TODO: Use original sequence characters rather than re-translated
 503     // characters in output
 504     // Find the biggest id length for formatting purposes
 505     String s1id = s1.getName(), s2id = s2.getName();
 506     int maxid = s1.getName().length();
 507     if (s2.getName().length() > maxid)
 508     {
 509       maxid = s2.getName().length();
 510     }
 511     if (maxid > 30)
 512     {
 513       maxid = 30;
 514       // JAL-527 - truncate the sequence ids
 515       if (s1.getName().length() > maxid)
 516       {
 517         s1id = s1.getName().substring(0, 30);
 518       }
 519       if (s2.getName().length() > maxid)
 520       {
 521         s2id = s2.getName().substring(0, 30);
 522       }
 523     }
 524     int len = 72 - maxid - 1;
 525     int nochunks = ((aseq1.length - count) / len) + 1;
 526     pid = 0;
 527
 528     output.append("Score = " + score[maxi][maxj] + "\n");
 529     output.append("Length of alignment = " + (aseq1.length - count) + "\n");
 530     output.append("Sequence ");
 531     output.append(new Format("%" + maxid + "s").form(s1.getName()));
 532     output.append(" :  " + s1.getStart() + " - " + s1.getEnd()
 533             + " (Sequence length = " + s1str.length() + ")\n");
 534     output.append("Sequence ");
 535     output.append(new Format("%" + maxid + "s").form(s2.getName()));
 536     output.append(" :  " + s2.getStart() + " - " + s2.getEnd()
 537             + " (Sequence length = " + s2str.length() + ")\n\n");
 538
 539     for (int j = 0; j < nochunks; j++)
 540     {
 541       // Print the first aligned sequence
 542       output.append(new Format("%" + (maxid) + "s").form(s1id) + " ");
 543
 544       for (int i = 0; i < len; i++)
 545       {
 546         if ((i + (j * len)) < astr1.length())
 547         {
 548           output.append(astr1.charAt(i + (j * len)));
 549         }
 550       }
 551
 552       output.append("\n");
 553       output.append(new Format("%" + (maxid) + "s").form(" ") + " ");
 554
 555       // Print out the matching chars
 556       for (int i = 0; i < len; i++)
 557       {
 558         if ((i + (j * len)) < astr1.length())
 559         {
 560           if (astr1.charAt(i + (j * len)) == astr2.charAt(i + (j * len))
 561                   && !jalview.util.Comparison.isGap(astr1.charAt(i
 562                           + (j * len))))
 563           {
 564             pid++;
 565             output.append("|");
 566           }
 567           else if (type.equals("pep"))
 568           {
 569             if (ResidueProperties.getPAM250(astr1.charAt(i + (j * len)),
 570                     astr2.charAt(i + (j * len))) > 0)
 571             {
 572               output.append(".");
 573             }
 574             else
 575             {
 576               output.append(" ");
 577             }
 578           }
 579           else
 580           {
 581             output.append(" ");
 582           }
 583         }
 584       }
 585
 586       // Now print the second aligned sequence
 587       output = output.append("\n");
 588       output = output.append(new Format("%" + (maxid) + "s").form(s2id)
 589               + " ");
 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 = output.append("\n\n");
 600     }
 601
 602     pid = pid / (float) (aseq1.length - count) * 100;
 603     output = output.append(new Format("Percentage ID = %2.2f\n\n")
 604             .form(pid));
 605
 606     try
 607     {
 608       os.print(output.toString());
 609     } catch (Exception ex)
 610     {
 611     }
 612   }
 613
 614   /**
 615    * DOCUMENT ME!
 616    *
 617    * @param mat
 618    *          DOCUMENT ME!
 619    */
 620   public void printScoreMatrix(int[][] mat)
 621   {
 622     int n = seq1.length;
 623     int m = seq2.length;
 624
 625     for (int i = 0; i < n; i++)
 626     {
 627       // Print the top sequence
 628       if (i == 0)
 629       {
 630         Format.print(System.out, "%8s", s2str.substring(0, 1));
 631
 632         for (int jj = 1; jj < m; jj++)
 633         {
 634           Format.print(System.out, "%5s", s2str.substring(jj, jj + 1));
 635         }
 636
 637         System.out.println();
 638       }
 639
 640       for (int j = 0; j < m; j++)
 641       {
 642         if (j == 0)
 643         {
 644           Format.print(System.out, "%3s", s1str.substring(i, i + 1));
 645         }
 646
 647         Format.print(System.out, "%3d ", mat[i][j] / 10);
 648       }
 649
 650       System.out.println();
 651     }
 652   }
 653
 654   /**
 655    * DOCUMENT ME!
 656    *
 657    * @param i
 658    *          DOCUMENT ME!
 659    * @param j
 660    *          DOCUMENT ME!
 661    *
 662    * @return DOCUMENT ME!
 663    */
 664   public int findTrace(int i, int j)
 665   {
 666     int t = 0;
 667     int max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
 668
 669     if (F[i][j] > max)
 670     {
 671       max = F[i][j];
 672       t = -1;
 673     }
 674     else if (F[i][j] == max)
 675     {
 676       if (prev == -1)
 677       {
 678         max = F[i][j];
 679         t = -1;
 680       }
 681     }
 682
 683     if (E[i][j] >= max)
 684     {
 685       max = E[i][j];
 686       t = 1;
 687     }
 688     else if (E[i][j] == max)
 689     {
 690       if (prev == 1)
 691       {
 692         max = E[i][j];
 693         t = 1;
 694       }
 695     }
 696
 697     prev = t;
 698
 699     return t;
 700   }
 701
 702   /**
 703    * DOCUMENT ME!
 704    */
 705   public void calcScoreMatrix()
 706   {
 707     int n = seq1.length;
 708     int m = seq2.length;
 709
 710     // top left hand element
 711     score[0][0] = lookup[seq1[0]][seq2[0]] * 10;
 712     E[0][0] = -gapExtend;
 713     F[0][0] = 0;
 714
 715     // Calculate the top row first
 716     for (int j = 1; j < m; j++)
 717     {
 718       // What should these values be? 0 maybe
 719       E[0][j] = max(score[0][j - 1] - gapOpen, E[0][j - 1] - gapExtend);
 720       F[0][j] = -gapExtend;
 721
 722       score[0][j] = max(lookup[seq1[0]][seq2[j]] * 10, -gapOpen, -gapExtend);
 723
 724       traceback[0][j] = 1;
 725     }
 726
 727     // Now do the left hand column
 728     for (int i = 1; i < n; i++)
 729     {
 730       E[i][0] = -gapOpen;
 731       F[i][0] = max(score[i - 1][0] - gapOpen, F[i - 1][0] - gapExtend);
 732
 733       score[i][0] = max(lookup[seq1[i]][seq2[0]] * 10, E[i][0], F[i][0]);
 734       traceback[i][0] = -1;
 735     }
 736
 737     // Now do all the other rows
 738     for (int i = 1; i < n; i++)
 739     {
 740       for (int j = 1; j < m; j++)
 741       {
 742         E[i][j] = max(score[i][j - 1] - gapOpen, E[i][j - 1] - gapExtend);
 743         F[i][j] = max(score[i - 1][j] - gapOpen, F[i - 1][j] - gapExtend);
 744
 745         score[i][j] = max(score[i - 1][j - 1]
 746                 + (lookup[seq1[i]][seq2[j]] * 10), E[i][j], F[i][j]);
 747         traceback[i][j] = findTrace(i, j);
 748       }
 749     }
 750   }
 751
 752   /**
 753    * DOCUMENT ME!
 754    *
 755    * @param gapChar
 756    *          DOCUMENT ME!
 757    * @param seq
 758    *          DOCUMENT ME!
 759    *
 760    * @return DOCUMENT ME!
 761    */
 762   public static String extractGaps(String gapChar, String seq)
 763   {
 764     StringTokenizer str = new StringTokenizer(seq, gapChar);
 765     StringBuffer newString = new StringBuffer();
 766
 767     while (str.hasMoreTokens())
 768     {
 769       newString.append(str.nextToken());
 770     }
 771
 772     return newString.toString();
 773   }
 774
 775   /**
 776    * DOCUMENT ME!
 777    *
 778    * @param i1
 779    *          DOCUMENT ME!
 780    * @param i2
 781    *          DOCUMENT ME!
 782    * @param i3
 783    *          DOCUMENT ME!
 784    *
 785    * @return DOCUMENT ME!
 786    */
 787   public int max(int i1, int i2, int i3)
 788   {
 789     int max = i1;
 790
 791     if (i2 > i1)
 792     {
 793       max = i2;
 794     }
 795
 796     if (i3 > max)
 797     {
 798       max = i3;
 799     }
 800
 801     return max;
 802   }
 803
 804   /**
 805    * DOCUMENT ME!
 806    *
 807    * @param i1
 808    *          DOCUMENT ME!
 809    * @param i2
 810    *          DOCUMENT ME!
 811    *
 812    * @return DOCUMENT ME!
 813    */
 814   public int max(int i1, int i2)
 815   {
 816     int max = i1;
 817
 818     if (i2 > i1)
 819     {
 820       max = i2;
 821     }
 822
 823     return max;
 824   }
 825
 826   /**
 827    * DOCUMENT ME!
 828    *
 829    * @param s
 830    *          DOCUMENT ME!
 831    * @param type
 832    *          DOCUMENT ME!
 833    *
 834    * @return DOCUMENT ME!
 835    */
 836   public int[] stringToInt(String s, String type)
 837   {
 838     int[] seq1 = new int[s.length()];
 839
 840     for (int i = 0; i < s.length(); i++)
 841     {
 842       // String ss = s.substring(i, i + 1).toUpperCase();
 843       char c = s.charAt(i);
 844       if ('a' <= c && c <= 'z')
 845       {
 846         // TO UPPERCASE !!!
 847         c -= ('a' - 'A');
 848       }
 849
 850       try
 851       {
 852         seq1[i] = charToInt[c]; // set accordingly from setType
 853         if (seq1[i] < 0 || seq1[i] > defInt) // set from setType: 23 for
 854                                              // peptides, or 4 for NA.
 855         {
 856           seq1[i] = defInt;
 857         }
 858
 859       } catch (Exception e)
 860       {
 861         seq1[i] = defInt;
 862       }
 863     }
 864
 865     return seq1;
 866   }
 867
 868   /**
 869    * DOCUMENT ME!
 870    *
 871    * @param g
 872    *          DOCUMENT ME!
 873    * @param mat
 874    *          DOCUMENT ME!
 875    * @param n
 876    *          DOCUMENT ME!
 877    * @param m
 878    *          DOCUMENT ME!
 879    * @param psize
 880    *          DOCUMENT ME!
 881    */
 882   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 883           int psize)
 884   {
 885     int max = -1000;
 886     int min = 1000;
 887
 888     for (int i = 0; i < n; i++)
 889     {
 890       for (int j = 0; j < m; j++)
 891       {
 892         if (mat[i][j] >= max)
 893         {
 894           max = mat[i][j];
 895         }
 896
 897         if (mat[i][j] <= min)
 898         {
 899           min = mat[i][j];
 900         }
 901       }
 902     }
 903
 904     System.out.println(max + " " + min);
 905
 906     for (int i = 0; i < n; i++)
 907     {
 908       for (int j = 0; j < m; j++)
 909       {
 910         int x = psize * i;
 911         int y = psize * j;
 912
 913         // System.out.println(mat[i][j]);
 914         float score = (float) (mat[i][j] - min) / (float) (max - min);
 915         g.setColor(new Color(score, 0, 0));
 916         g.fillRect(x, y, psize, psize);
 917
 918         // System.out.println(x + " " + y + " " + score);
 919       }
 920     }
 921   }
 922
 923   /**
 924    * compute the PID vector used by the redundancy filter.
 925    *
 926    * @param originalSequences
 927    *          - sequences in alignment that are to filtered
 928    * @param omitHidden
 929    *          - null or strings to be analysed (typically, visible portion of
 930    *          each sequence in alignment)
 931    * @param start
 932    *          - first column in window for calculation
 933    * @param end
 934    *          - last column in window for calculation
 935    * @param ungapped
 936    *          - if true then use ungapped sequence to compute PID
 937    * @return vector containing maximum PID for i-th sequence and any sequences
 938    *         longer than that seuqence
 939    */
 940   public static float[] computeRedundancyMatrix(
 941           SequenceI[] originalSequences, String[] omitHidden, int start,
 942           int end, boolean ungapped)
 943   {
 944     int height = originalSequences.length;
 945     float[] redundancy = new float[height];
 946     int[] lngth = new int[height];
 947     for (int i = 0; i < height; i++)
 948     {
 949       redundancy[i] = 0f;
 950       lngth[i] = -1;
 951     }
 952
 953     // long start = System.currentTimeMillis();
 954
 955     float pid;
 956     String seqi, seqj;
 957     for (int i = 0; i < height; i++)
 958     {
 959
 960       for (int j = 0; j < i; j++)
 961       {
 962         if (i == j)
 963         {
 964           continue;
 965         }
 966
 967         if (omitHidden == null)
 968         {
 969           seqi = originalSequences[i].getSequenceAsString(start, end);
 970           seqj = originalSequences[j].getSequenceAsString(start, end);
 971         }
 972         else
 973         {
 974           seqi = omitHidden[i];
 975           seqj = omitHidden[j];
 976         }
 977         if (lngth[i] == -1)
 978         {
 979           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
 980           lngth[i] = ug.length();
 981           if (ungapped)
 982           {
 983             seqi = ug;
 984           }
 985         }
 986         if (lngth[j] == -1)
 987         {
 988           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
 989           lngth[j] = ug.length();
 990           if (ungapped)
 991           {
 992             seqj = ug;
 993           }
 994         }
 995         pid = Comparison.PID(seqi, seqj);
 996
 997         // use real sequence length rather than string length
 998         if (lngth[j] < lngth[i])
 999         {
1000           redundancy[j] = Math.max(pid, redundancy[j]);
1001         }
1002         else
1003         {
1004           redundancy[i] = Math.max(pid, redundancy[i]);
1005         }
1006
1007       }
1008     }
1009     return redundancy;
1010   }
1011 }