JAL-2416 order score models by order of addition rather than name
[jalview.git] / src / jalview / util / Comparison.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.util;
22
23 import jalview.datamodel.SequenceI;
24
25 import java.util.ArrayList;
26 import java.util.List;
27
28 /**
29  * Assorted methods for analysing or comparing sequences.
30  */
31 public class Comparison
32 {
33   private static final int EIGHTY_FIVE = 85;
34
35   private static final int TO_UPPER_CASE = 'a' - 'A';
36
37   private static final char GAP_SPACE = ' ';
38
39   private static final char GAP_DOT = '.';
40
41   private static final char GAP_DASH = '-';
42
43   public static final String GapChars = new String(new char[] { GAP_SPACE,
44       GAP_DOT, GAP_DASH });
45
46   /**
47    * DOCUMENT ME!
48    * 
49    * @param ii
50    *          DOCUMENT ME!
51    * @param jj
52    *          DOCUMENT ME!
53    * 
54    * @return DOCUMENT ME!
55    */
56   public static final float compare(SequenceI ii, SequenceI jj)
57   {
58     return Comparison.compare(ii, jj, 0, ii.getLength() - 1);
59   }
60
61   /**
62    * this was supposed to be an ungapped pid calculation
63    * 
64    * @param ii
65    *          SequenceI
66    * @param jj
67    *          SequenceI
68    * @param start
69    *          int
70    * @param end
71    *          int
72    * @return float
73    */
74   public static float compare(SequenceI ii, SequenceI jj, int start, int end)
75   {
76     String si = ii.getSequenceAsString();
77     String sj = jj.getSequenceAsString();
78
79     int ilen = si.length() - 1;
80     int jlen = sj.length() - 1;
81
82     while (Comparison.isGap(si.charAt(start + ilen)))
83     {
84       ilen--;
85     }
86
87     while (Comparison.isGap(sj.charAt(start + jlen)))
88     {
89       jlen--;
90     }
91
92     int count = 0;
93     int match = 0;
94     float pid = -1;
95
96     if (ilen > jlen)
97     {
98       for (int j = 0; j < jlen; j++)
99       {
100         if (si.substring(start + j, start + j + 1).equals(
101                 sj.substring(start + j, start + j + 1)))
102         {
103           match++;
104         }
105
106         count++;
107       }
108
109       pid = (float) match / (float) ilen * 100;
110     }
111     else
112     {
113       for (int j = 0; j < jlen; j++)
114       {
115         if (si.substring(start + j, start + j + 1).equals(
116                 sj.substring(start + j, start + j + 1)))
117         {
118           match++;
119         }
120
121         count++;
122       }
123
124       pid = (float) match / (float) jlen * 100;
125     }
126
127     return pid;
128   }
129
130   /**
131    * this is a gapped PID calculation
132    * 
133    * @param s1
134    *          SequenceI
135    * @param s2
136    *          SequenceI
137    * @return float
138    */
139   public final static float PID(String seq1, String seq2)
140   {
141     return PID(seq1, seq2, 0, seq1.length());
142   }
143
144   static final int caseShift = 'a' - 'A';
145
146   // Another pid with region specification
147   public final static float PID(String seq1, String seq2, int start, int end)
148   {
149     return PID(seq1, seq2, start, end, true, false);
150   }
151
152   /**
153    * Calculate percent identity for a pair of sequences over a particular range,
154    * with different options for ignoring gaps.
155    * 
156    * @param seq1
157    * @param seq2
158    * @param start
159    *          - position in seqs
160    * @param end
161    *          - position in seqs
162    * @param wcGaps
163    *          - if true - gaps match any character, if false, do not match
164    *          anything
165    * @param ungappedOnly
166    *          - if true - only count PID over ungapped columns
167    * @return
168    */
169   public final static float PID(String seq1, String seq2, int start,
170           int end, boolean wcGaps, boolean ungappedOnly)
171   {
172     int s1len = seq1.length();
173     int s2len = seq2.length();
174
175     int len = Math.min(s1len, s2len);
176
177     if (end < len)
178     {
179       len = end;
180     }
181
182     if (len < start)
183     {
184       start = len - 1; // we just use a single residue for the difference
185     }
186
187     int elen = len - start, bad = 0;
188     char chr1;
189     char chr2;
190     boolean agap;
191     for (int i = start; i < len; i++)
192     {
193       chr1 = seq1.charAt(i);
194
195       chr2 = seq2.charAt(i);
196       agap = isGap(chr1) || isGap(chr2);
197       if ('a' <= chr1 && chr1 <= 'z')
198       {
199         // TO UPPERCASE !!!
200         // Faster than toUpperCase
201         chr1 -= caseShift;
202       }
203       if ('a' <= chr2 && chr2 <= 'z')
204       {
205         // TO UPPERCASE !!!
206         // Faster than toUpperCase
207         chr2 -= caseShift;
208       }
209
210       if (chr1 != chr2)
211       {
212         if (agap)
213         {
214           if (ungappedOnly)
215           {
216             elen--;
217           }
218           else if (!wcGaps)
219           {
220             bad++;
221           }
222         }
223         else
224         {
225           bad++;
226         }
227       }
228
229     }
230     if (elen < 1)
231     {
232       return 0f;
233     }
234     return ((float) 100 * (elen - bad)) / elen;
235   }
236
237   /**
238    * Answers true if the supplied character is a recognised gap character, else
239    * false. Currently hard-coded to recognise '-', '-' or ' ' (hyphen / dot /
240    * space).
241    * 
242    * @param c
243    * 
244    * @return
245    */
246   public static final boolean isGap(char c)
247   {
248     return (c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE) ? true : false;
249   }
250
251   /**
252    * Overloaded method signature to test whether a single sequence is nucleotide
253    * (that is, more than 85% CGTA)
254    * 
255    * @param seq
256    * @return
257    */
258   public static final boolean isNucleotide(SequenceI seq)
259   {
260     return isNucleotide(new SequenceI[] { seq });
261   }
262
263   /**
264    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
265    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
266    * wrong answer (as AGCT are also amino acid codes).
267    * 
268    * @param seqs
269    * @return
270    */
271   public static final boolean isNucleotide(SequenceI[] seqs)
272   {
273     if (seqs == null)
274     {
275       return false;
276     }
277     char[][] letters = new char[seqs.length][];
278     for (int i = 0; i < seqs.length; i++)
279     {
280       if (seqs[i] != null)
281       {
282         char[] sequence = seqs[i].getSequence();
283         if (sequence != null)
284         {
285           letters[i] = sequence;
286         }
287       }
288     }
289
290     return areNucleotide(letters);
291   }
292
293   /**
294    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
295    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
296    * wrong answer (as AGCT are also amino acid codes).
297    * 
298    * @param letters
299    * @return
300    */
301   static final boolean areNucleotide(char[][] letters)
302   {
303     int ntCount = 0;
304     int aaCount = 0;
305     for (char[] seq : letters)
306     {
307       if (seq == null)
308       {
309         continue;
310       }
311       // TODO could possibly make an informed guess just from the first sequence
312       // to save a lengthy calculation
313       for (char c : seq)
314       {
315         if (isNucleotide(c))
316         {
317           ntCount++;
318         }
319         else if (!isGap(c))
320         {
321           aaCount++;
322         }
323       }
324     }
325
326     /*
327      * Check for nucleotide count > 85% of total count (in a form that evades
328      * int / float conversion or divide by zero).
329      */
330     if (ntCount * 100 > EIGHTY_FIVE * (ntCount + aaCount))
331     {
332       return true;
333     }
334     else
335     {
336       return false;
337     }
338
339   }
340
341   /**
342    * Answers true if the character is one of aAcCgGtTuU
343    * 
344    * @param c
345    * @return
346    */
347   public static boolean isNucleotide(char c)
348   {
349     if ('a' <= c && c <= 'z')
350     {
351       c -= TO_UPPER_CASE;
352     }
353
354     switch (c)
355     {
356     case 'A':
357     case 'C':
358     case 'G':
359     case 'T':
360     case 'U':
361       return true;
362     }
363     return false;
364   }
365
366   /**
367    * Answers true if every character in the string is one of aAcCgGtTuU, or
368    * (optionally) a gap character (dot, dash, space), else false
369    * 
370    * @param s
371    * @param allowGaps
372    * @return
373    */
374   public static boolean isNucleotideSequence(String s, boolean allowGaps)
375   {
376     if (s == null)
377     {
378       return false;
379     }
380     for (int i = 0; i < s.length(); i++)
381     {
382       char c = s.charAt(i);
383       if (!isNucleotide(c))
384       {
385         if (!allowGaps || !isGap(c))
386         {
387           return false;
388         }
389       }
390     }
391     return true;
392   }
393
394   /**
395    * Convenience overload of isNucleotide
396    * 
397    * @param seqs
398    * @return
399    */
400   public static boolean isNucleotide(SequenceI[][] seqs)
401   {
402     if (seqs == null)
403     {
404       return false;
405     }
406     List<SequenceI> flattened = new ArrayList<SequenceI>();
407     for (SequenceI[] ss : seqs)
408     {
409       for (SequenceI s : ss)
410       {
411         flattened.add(s);
412       }
413     }
414     final SequenceI[] oneDArray = flattened.toArray(new SequenceI[flattened
415             .size()]);
416     return isNucleotide(oneDArray);
417   }
418
419   /**
420    * Compares two residues either case sensitively or case insensitively
421    * depending on the caseSensitive flag
422    * 
423    * @param c1
424    *          first char
425    * @param c2
426    *          second char to compare with
427    * @param caseSensitive
428    *          if true comparison will be case sensitive otherwise its not
429    * @return
430    */
431   public static boolean isSameResidue(char c1, char c2,
432           boolean caseSensitive)
433   {
434     if (caseSensitive)
435     {
436       return (c1 == c2);
437     }
438     else
439     {
440       return Character.toUpperCase(c1) == Character.toUpperCase(c2);
441     }
442   }
443 }