Merge branch 'releases/Release_2_11_4_Branch'
[jalview.git] / 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   public static final char GAP_SPACE = ' ';
38
39   public static final char GAP_DOT = '.';
40
41   public 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    * @deprecated use PIDModel.computePID()
139    */
140   @Deprecated
141   public final static float PID(String seq1, String seq2)
142   {
143     return PID(seq1, seq2, 0, seq1.length());
144   }
145
146   static final int caseShift = 'a' - 'A';
147
148   // Another pid with region specification
149   /**
150    * @deprecated use PIDModel.computePID()
151    */
152   @Deprecated
153   public final static float PID(String seq1, String seq2, int start, int end)
154   {
155     return PID(seq1, seq2, start, end, true, false);
156   }
157
158   /**
159    * Calculate percent identity for a pair of sequences over a particular range,
160    * with different options for ignoring gaps.
161    * 
162    * @param seq1
163    * @param seq2
164    * @param start
165    *          - position in seqs
166    * @param end
167    *          - position in seqs
168    * @param wcGaps
169    *          - if true - gaps match any character, if false, do not match
170    *          anything
171    * @param ungappedOnly
172    *          - if true - only count PID over ungapped columns
173    * @return
174    * @deprecated use PIDModel.computePID()
175    */
176   @Deprecated
177   public final static float PID(String seq1, String seq2, int start,
178           int end, boolean wcGaps, boolean ungappedOnly)
179   {
180     int s1len = seq1.length();
181     int s2len = seq2.length();
182
183     int len = Math.min(s1len, s2len);
184
185     if (end < len)
186     {
187       len = end;
188     }
189
190     if (len < start)
191     {
192       start = len - 1; // we just use a single residue for the difference
193     }
194
195     int elen = len - start, bad = 0;
196     char chr1;
197     char chr2;
198     boolean agap;
199     for (int i = start; i < len; i++)
200     {
201       chr1 = seq1.charAt(i);
202
203       chr2 = seq2.charAt(i);
204       agap = isGap(chr1) || isGap(chr2);
205       if ('a' <= chr1 && chr1 <= 'z')
206       {
207         // TO UPPERCASE !!!
208         // Faster than toUpperCase
209         chr1 -= caseShift;
210       }
211       if ('a' <= chr2 && chr2 <= 'z')
212       {
213         // TO UPPERCASE !!!
214         // Faster than toUpperCase
215         chr2 -= caseShift;
216       }
217
218       if (chr1 != chr2)
219       {
220         if (agap)
221         {
222           if (ungappedOnly)
223           {
224             elen--;
225           }
226           else if (!wcGaps)
227           {
228             bad++;
229           }
230         }
231         else
232         {
233           bad++;
234         }
235       }
236
237     }
238     if (elen < 1)
239     {
240       return 0f;
241     }
242     return ((float) 100 * (elen - bad)) / elen;
243   }
244
245   /**
246    * Answers true if the supplied character is a recognised gap character, else
247    * false. Currently hard-coded to recognise '-', '-' or ' ' (hyphen / dot /
248    * space).
249    * 
250    * @param c
251    * 
252    * @return
253    */
254   public static final boolean isGap(char c)
255   {
256     return (c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE) ? true : false;
257   }
258
259   /**
260    * Overloaded method signature to test whether a single sequence is nucleotide
261    * (that is, more than 85% CGTA)
262    * 
263    * @param seq
264    * @return
265    */
266   public static final boolean isNucleotide(SequenceI seq)
267   {
268     return isNucleotide(new SequenceI[] { seq });
269   }
270
271   /**
272    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
273    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
274    * wrong answer (as AGCT are also amino acid codes).
275    * 
276    * @param seqs
277    * @return
278    */
279   public static final boolean isNucleotide(SequenceI[] seqs)
280   {
281     if (seqs == null)
282     {
283       return false;
284     }
285     char[][] letters = new char[seqs.length][];
286     for (int i = 0; i < seqs.length; i++)
287     {
288       if (seqs[i] != null)
289       {
290         char[] sequence = seqs[i].getSequence();
291         if (sequence != null)
292         {
293           letters[i] = sequence;
294         }
295       }
296     }
297
298     return areNucleotide(letters);
299   }
300
301   /**
302    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
303    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
304    * wrong answer (as AGCT are also amino acid codes).
305    * 
306    * @param letters
307    * @return
308    */
309   static final boolean areNucleotide(char[][] letters)
310   {
311     int ntCount = 0;
312     int aaCount = 0;
313     for (char[] seq : letters)
314     {
315       if (seq == null)
316       {
317         continue;
318       }
319       // TODO could possibly make an informed guess just from the first sequence
320       // to save a lengthy calculation
321       for (char c : seq)
322       {
323         if (isNucleotide(c))
324         {
325           ntCount++;
326         }
327         else if (!isGap(c))
328         {
329           aaCount++;
330         }
331       }
332     }
333
334     /*
335      * Check for nucleotide count > 85% of total count (in a form that evades
336      * int / float conversion or divide by zero).
337      */
338     if (ntCount * 100 > EIGHTY_FIVE * (ntCount + aaCount))
339     {
340       return true;
341     }
342     else
343     {
344       return false;
345     }
346
347   }
348
349   /**
350    * Answers true if the character is one of aAcCgGtTuU
351    * 
352    * @param c
353    * @return
354    */
355   public static boolean isNucleotide(char c)
356   {
357     if ('a' <= c && c <= 'z')
358     {
359       c -= TO_UPPER_CASE;
360     }
361
362     switch (c)
363     {
364     case 'A':
365     case 'C':
366     case 'G':
367     case 'T':
368     case 'U':
369       return true;
370     }
371     return false;
372   }
373
374   /**
375    * Answers true if every character in the string is one of aAcCgGtTuU, or
376    * (optionally) a gap character (dot, dash, space), else false
377    * 
378    * @param s
379    * @param allowGaps
380    * @return
381    */
382   public static boolean isNucleotideSequence(String s, boolean allowGaps)
383   {
384     if (s == null)
385     {
386       return false;
387     }
388     for (int i = 0; i < s.length(); i++)
389     {
390       char c = s.charAt(i);
391       if (!isNucleotide(c))
392       {
393         if (!allowGaps || !isGap(c))
394         {
395           return false;
396         }
397       }
398     }
399     return true;
400   }
401
402   /**
403    * Convenience overload of isNucleotide
404    * 
405    * @param seqs
406    * @return
407    */
408   public static boolean isNucleotide(SequenceI[][] seqs)
409   {
410     if (seqs == null)
411     {
412       return false;
413     }
414     List<SequenceI> flattened = new ArrayList<SequenceI>();
415     for (SequenceI[] ss : seqs)
416     {
417       for (SequenceI s : ss)
418       {
419         flattened.add(s);
420       }
421     }
422     final SequenceI[] oneDArray = flattened.toArray(new SequenceI[flattened
423             .size()]);
424     return isNucleotide(oneDArray);
425   }
426
427   /**
428    * Compares two residues either case sensitively or case insensitively
429    * depending on the caseSensitive flag
430    * 
431    * @param c1
432    *          first char
433    * @param c2
434    *          second char to compare with
435    * @param caseSensitive
436    *          if true comparison will be case sensitive otherwise its not
437    * @return
438    */
439   public static boolean isSameResidue(char c1, char c2,
440           boolean caseSensitive)
441   {
442     if (caseSensitive)
443     {
444       return (c1 == c2);
445     }
446     else
447     {
448       return Character.toUpperCase(c1) == Character.toUpperCase(c2);
449     }
450   }
451 }