src/jalview/analysis/AlignSeq.java

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