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