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