src/jalview/analysis/StructureFrequency.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.6)
   3  * Copyright (C) 2010 J Procter, AM Waterhouse, G Barton, M Clamp, S Searle
   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 of the License, or (at your option) any later version.
  10  *
  11  * Jalview is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty
  13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14  * PURPOSE.  See the GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  17  */
  18 package jalview.analysis;
  19
  20 import java.util.*;
  21
  22 import jalview.datamodel.*;
  23
  24 /**
  25  * Takes in a vector or array of sequences and column start and column end and
  26  * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied.
  27  * This class is used extensively in calculating alignment colourschemes that
  28  * depend on the amount of conservation in each alignment column.
  29  *
  30  * @author $author$
  31  * @version $Revision$
  32  */
  33 public class StructureFrequency
  34 {
  35   // No need to store 1000s of strings which are not
  36   // visible to the user.
  37   public static final String MAXCOUNT = "C";
  38
  39   public static final String MAXRESIDUE = "R";
  40
  41   public static final String PID_GAPS = "G";
  42
  43   public static final String PID_NOGAPS = "N";
  44
  45   public static final String PROFILE = "P";
  46
  47   public static final String PAIRPROFILE = "B";
  48
  49   /**
  50    * Returns the 3' position of a base pair
  51    *
  52    * @param pairs
  53    *          Secondary structure annotation
  54    * @param indice
  55    *          5' position of a base pair
  56    * @return 3' position of a base pair
  57    */
  58   public static int findPair(SequenceFeature[] pairs, int indice)
  59   {
  60     for (int i = 0; i < pairs.length; i++)
  61     {
  62       if (pairs[i].getBegin() == indice)
  63       {
  64         return pairs[i].getEnd();
  65       }
  66     }
  67     return -1;
  68   }
  69
  70   /**
  71    * Method to calculate a 'base pair consensus row', very similar to nucleotide
  72    * consensus but takes into account a given structure
  73    *
  74    * @param sequences
  75    * @param start
  76    * @param end
  77    * @param result
  78    * @param profile
  79    * @param rnaStruc
  80    */
  81   public static final void calculate(SequenceI[] sequences, int start,
  82           int end, Hashtable[] result, boolean profile,
  83           AlignmentAnnotation rnaStruc)
  84   {
  85     Hashtable residueHash;
  86     String maxResidue;
  87     char[] seq, struc = rnaStruc.getRNAStruc().toCharArray();
  88     SequenceFeature[] rna = rnaStruc._rnasecstr;
  89     char c, s, cEnd;
  90     int count, nonGap = 0, i, bpEnd = -1, j, jSize = sequences.length;
  91     int[] values;
  92     int[][] pairs;
  93     float percentage;
  94
  95     for (i = start; i < end; i++) // foreach column
  96     {
  97       residueHash = new Hashtable();
  98       maxResidue = "-";
  99       values = new int[255];
 100       pairs = new int[255][255];
 101       bpEnd = -1;
 102       if (i < struc.length)
 103       {
 104         s = struc[i];
 105       }
 106       else
 107       {
 108         s = '-';
 109       }
 110       if (s == '.' || s == ' ')
 111       {
 112         s = '-';
 113       }
 114
 115       if (s != '(')
 116       {
 117         if (s == '-')
 118         {
 119           values['-']++;
 120         }
 121       }
 122       else
 123       {
 124         for (j = 0; j < jSize; j++) // foreach row
 125         {
 126           if (sequences[j] == null)
 127           {
 128             System.err
 129                     .println("WARNING: Consensus skipping null sequence - possible race condition.");
 130             continue;
 131           }
 132           seq = sequences[j].getSequence();
 133
 134           if (seq.length > i)
 135           {
 136             c = seq[i];
 137
 138             // standard representation for gaps in sequence and structure
 139             if (c == '.' || c == ' ')
 140             {
 141               c = '-';
 142             }
 143
 144             if (c == '-')
 145             {
 146               values['-']++;
 147               continue;
 148             }
 149             bpEnd = findPair(rna, i);
 150             cEnd = seq[bpEnd];
 151             if (checkBpType(c, cEnd))
 152             {
 153               values['(']++; // H means it's a helix (structured)
 154             }
 155             pairs[c][cEnd]++;
 156
 157             maxResidue = "(";
 158           }
 159         }
 160         // nonGap++;
 161       }
 162       // UPDATE this for new values
 163       if (profile)
 164       {
 165         residueHash.put(PROFILE, new int[][]
 166         { values, new int[]
 167         { jSize, (jSize - values['-']) } });
 168
 169         residueHash.put(PAIRPROFILE, pairs);
 170       }
 171
 172       count = values['('];
 173
 174       residueHash.put(MAXCOUNT, new Integer(count));
 175       residueHash.put(MAXRESIDUE, maxResidue);
 176
 177       percentage = ((float) count * 100) / (float) jSize;
 178       residueHash.put(PID_GAPS, new Float(percentage));
 179
 180       // percentage = ((float) count * 100) / (float) nongap;
 181       // residueHash.put(PID_NOGAPS, new Float(percentage));
 182       if (result[i] == null)
 183       {
 184         result[i] = residueHash;
 185       }
 186       if (bpEnd > 0)
 187       {
 188         values[')'] = values['('];
 189         values['('] = 0;
 190
 191         residueHash = new Hashtable();
 192         maxResidue = ")";
 193
 194         if (profile)
 195         {
 196           residueHash.put(PROFILE, new int[][]
 197           { values, new int[]
 198           { jSize, (jSize - values['-']) } });
 199
 200           residueHash.put(PAIRPROFILE, pairs);
 201         }
 202
 203         residueHash.put(MAXCOUNT, new Integer(count));
 204         residueHash.put(MAXRESIDUE, maxResidue);
 205
 206         percentage = ((float) count * 100) / (float) jSize;
 207         residueHash.put(PID_GAPS, new Float(percentage));
 208
 209         result[bpEnd] = residueHash;
 210       }
 211     }
 212   }
 213
 214   /**
 215    * Method to check if a base-pair is a canonical or a wobble bp
 216    *
 217    * @param up
 218    *          5' base
 219    * @param down
 220    *          3' base
 221    * @return True if it is a canonical/wobble bp
 222    */
 223   public static boolean checkBpType(char up, char down)
 224   {
 225     if (up > 'Z')
 226     {
 227       up -= 32;
 228     }
 229     if (down > 'Z')
 230     {
 231       down -= 32;
 232     }
 233
 234     switch (up)
 235     {
 236     case 'A':
 237       switch (down)
 238       {
 239       case 'T':
 240         return true;
 241       case 'U':
 242         return true;
 243       }
 244       break;
 245     case 'C':
 246       switch (down)
 247       {
 248       case 'G':
 249         return true;
 250       }
 251       break;
 252     case 'T':
 253       switch (down)
 254       {
 255       case 'A':
 256         return true;
 257       case 'G':
 258         return true;
 259       }
 260       break;
 261     case 'G':
 262       switch (down)
 263       {
 264       case 'C':
 265         return true;
 266       case 'T':
 267         return true;
 268       case 'U':
 269         return true;
 270       }
 271       break;
 272     case 'U':
 273       switch (down)
 274       {
 275       case 'A':
 276         return true;
 277       case 'G':
 278         return true;
 279       }
 280       break;
 281     }
 282     return false;
 283   }
 284
 285   /**
 286    * Compute all or part of the annotation row from the given consensus
 287    * hashtable
 288    *
 289    * @param consensus
 290    *          - pre-allocated annotation row
 291    * @param hconsensus
 292    * @param iStart
 293    * @param width
 294    * @param ignoreGapsInConsensusCalculation
 295    * @param includeAllConsSymbols
 296    */
 297   public static void completeConsensus(AlignmentAnnotation consensus,
 298           Hashtable[] hconsensus, int iStart, int width,
 299           boolean ignoreGapsInConsensusCalculation,
 300           boolean includeAllConsSymbols)
 301   {
 302     float tval, value;
 303     if (consensus == null || consensus.annotations == null
 304             || consensus.annotations.length < width)
 305     {
 306       // called with a bad alignment annotation row - wait for it to be
 307       // initialised properly
 308       return;
 309     }
 310     for (int i = iStart; i < width; i++)
 311     {
 312       if (i >= hconsensus.length)
 313       {
 314         // happens if sequences calculated over were shorter than alignment
 315         // width
 316         consensus.annotations[i] = null;
 317         continue;
 318       }
 319       value = 0;
 320       if (ignoreGapsInConsensusCalculation)
 321       {
 322         value = ((Float) hconsensus[i].get(StructureFrequency.PID_NOGAPS))
 323                 .floatValue();
 324       }
 325       else
 326       {
 327         value = ((Float) hconsensus[i].get(StructureFrequency.PID_GAPS))
 328                 .floatValue();
 329       }
 330
 331       String maxRes = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
 332               .toString();
 333       String mouseOver = hconsensus[i].get(StructureFrequency.MAXRESIDUE)
 334               + " ";
 335       if (maxRes.length() > 1)
 336       {
 337         mouseOver = "[" + maxRes + "] ";
 338         maxRes = "+";
 339       }
 340       int[][] profile = (int[][]) hconsensus[i]
 341               .get(StructureFrequency.PROFILE);
 342       int[][] pairs = (int[][]) hconsensus[i]
 343               .get(StructureFrequency.PAIRPROFILE);
 344
 345       if (pairs != null && includeAllConsSymbols) // Just responsible for the
 346       // tooltip
 347       // TODO Update tooltips for Structure row
 348       {
 349         mouseOver = "";
 350
 351         /* TODO It's not sure what is the purpose of the alphabet and wheter it is useful for structure?
 352          *
 353          * if (alphabet != null) { for (int c = 0; c < alphabet.length; c++) {
 354          * tval = ((float) profile[0][alphabet[c]]) 100f / (float)
 355          * profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; mouseOver +=
 356          * ((c == 0) ? "" : "; ") + alphabet[c] + " " + ((int) tval) + "%"; } }
 357          * else {
 358          */
 359         Object[] ca = new Object[625];
 360         float[] vl = new float[625];
 361         int x = 0;
 362         for (int c = 65; c < 90; c++)
 363         {
 364           for (int d = 65; d < 90; d++)
 365           {
 366             ca[x] = new int[]
 367             { c, d };
 368             vl[x] = (float) pairs[c][d];
 369             x++;
 370           }
 371         }
 372         jalview.util.QuickSort.sort(vl, ca);
 373         int p = 0;
 374
 375         for (int c = 624; c > 0; c--)
 376         {
 377           if (vl[c] > 0)
 378           {
 379             tval = ((float) vl[c] * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
 380                     : 0]);
 381             mouseOver += ((p == 0) ? "" : "; ") + (char) ((int[]) ca[c])[0]
 382                     + (char) ((int[]) ca[c])[1] + " " + ((int) tval) + "%";
 383             p++;
 384
 385           }
 386         }
 387
 388         // }
 389       }
 390       else
 391       {
 392         mouseOver += ((int) value + "%");
 393       }
 394       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
 395               value);
 396     }
 397   }
 398
 399   /**
 400    * get the sorted base-pair profile for the given position of the consensus
 401    *
 402    * @param hconsensus
 403    * @return profile of the given column
 404    */
 405   public static int[] extractProfile(Hashtable hconsensus,
 406           boolean ignoreGapsInConsensusCalculation)
 407   {
 408     int[] rtnval = new int[51]; // 2*(5*5)+1
 409     int[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
 410     int[][] pairs = (int[][]) hconsensus
 411             .get(StructureFrequency.PAIRPROFILE);
 412
 413     if (profile == null)
 414       return null;
 415
 416     // TODO fix the object length, also do it in completeConsensus
 417     Object[] ca = new Object[625];
 418     float[] vl = new float[625];
 419     int x = 0;
 420     for (int c = 65; c < 90; c++)
 421     {
 422       for (int d = 65; d < 90; d++)
 423       {
 424         ca[x] = new int[]
 425         { c, d };
 426         vl[x] = (float) pairs[c][d];
 427         x++;
 428       }
 429     }
 430     jalview.util.QuickSort.sort(vl, ca);
 431
 432     rtnval[0] = 1;
 433     for (int c = 624; c > 0; c--)
 434     {
 435       if (vl[c] > 0)
 436       {
 437         rtnval[rtnval[0]++] = ((int[]) ca[c])[0];
 438         rtnval[rtnval[0]++] = ((int[]) ca[c])[1];
 439         rtnval[rtnval[0]++] = (int) ((float) vl[c] * 100f / (float) profile[1][ignoreGapsInConsensusCalculation ? 1
 440                 : 0]);
 441       }
 442     }
 443
 444     return rtnval;
 445   }
 446
 447   public static void main(String args[])
 448   {
 449     // Short test to see if checkBpType works
 450     ArrayList<String> test = new ArrayList<String>();
 451     test.add("A");
 452     test.add("c");
 453     test.add("g");
 454     test.add("T");
 455     test.add("U");
 456     for (String i : test)
 457     {
 458       for (String j : test)
 459       {
 460         System.out.println(i + "-" + j + ": "
 461                 + StructureFrequency.checkBpType(i.charAt(0), j.charAt(0)));
 462       }
 463     }
 464   }
 465 }