src/jalview/analysis/StructureFrequency.java

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