src/jalview/analysis/AlignSeq.java

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