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