src/jalview/analysis/StructureFrequency.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.1)
   3  * Copyright (C) 2014 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 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  * The Jalview Authors are detailed in the 'AUTHORS' file.
  18  */
  19 package jalview.analysis;
  20
  21 import java.util.*;
  22
  23 import jalview.util.Format;
  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, long nseq)
 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     String fmtstr="%3.1f";
 313     int precision=2;
 314     while (nseq>100) {
 315       precision++;
 316       nseq/=10;
 317     }
 318     if (precision>2)
 319     {
 320       fmtstr = "%"+(2+precision)+"."+precision+"f";
 321     }
 322     Format fmt = new Format(fmtstr);
 323
 324     for (int i = iStart; i < width; i++)
 325     {
 326       Hashtable hci;
 327       if (i >= hconsensus.length || ((hci = hconsensus[i]) == null))
 328       {
 329         // happens if sequences calculated over were shorter than alignment
 330         // width
 331         consensus.annotations[i] = null;
 332         continue;
 333       }
 334       value = 0;
 335       Float fv;
 336       if (ignoreGapsInConsensusCalculation)
 337       {
 338         fv = (Float) hci.get(StructureFrequency.PID_NOGAPS);
 339       }
 340       else
 341       {
 342         fv = (Float) hci.get(StructureFrequency.PID_GAPS);
 343       }
 344       if (fv == null)
 345       {
 346         consensus.annotations[i] = null;
 347         // data has changed below us .. give up and
 348         continue;
 349       }
 350       value = fv.floatValue();
 351       String maxRes = hci.get(StructureFrequency.MAXRESIDUE).toString();
 352       String mouseOver = hci.get(StructureFrequency.MAXRESIDUE) + " ";
 353       if (maxRes.length() > 1)
 354       {
 355         mouseOver = "[" + maxRes + "] ";
 356         maxRes = "+";
 357       }
 358       int[][] profile = (int[][]) hci.get(StructureFrequency.PROFILE);
 359       int[][] pairs = (int[][]) hci.get(StructureFrequency.PAIRPROFILE);
 360
 361       if (pairs != null && includeAllConsSymbols) // Just responsible for the
 362       // tooltip
 363       // TODO Update tooltips for Structure row
 364       {
 365         mouseOver = "";
 366
 367         /*
 368          * TODO It's not sure what is the purpose of the alphabet and wheter it
 369          * is useful for structure?
 370          *
 371          * if (alphabet != null) { for (int c = 0; c < alphabet.length; c++) {
 372          * tval = ((float) profile[0][alphabet[c]]) 100f / (float)
 373          * profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; mouseOver +=
 374          * ((c == 0) ? "" : "; ") + alphabet[c] + " " + ((int) tval) + "%"; } }
 375          * else {
 376          */
 377         Object[] ca = new Object[625];
 378         float[] vl = new float[625];
 379         int x = 0;
 380         for (int c = 65; c < 90; c++)
 381         {
 382           for (int d = 65; d < 90; d++)
 383           {
 384             ca[x] = new int[]
 385             { c, d };
 386             vl[x] = pairs[c][d];
 387             x++;
 388           }
 389         }
 390         jalview.util.QuickSort.sort(vl, ca);
 391         int p = 0;
 392
 393         for (int c = 624; c > 0; c--)
 394         {
 395           if (vl[c] > 0)
 396           {
 397             tval = (vl[c] * 100f / profile[1][ignoreGapsInConsensusCalculation ? 1
 398                     : 0]);
 399             mouseOver += ((p == 0) ? "" : "; ") + (char) ((int[]) ca[c])[0]
 400                     + (char) ((int[]) ca[c])[1] + " " + fmt.form(tval) + "%";
 401             p++;
 402
 403           }
 404         }
 405
 406         // }
 407       }
 408       else
 409       {
 410         mouseOver += (fmt.form(value) + "%");
 411       }
 412       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
 413               value);
 414     }
 415   }
 416
 417   /**
 418    * get the sorted base-pair profile for the given position of the consensus
 419    *
 420    * @param hconsensus
 421    * @return profile of the given column
 422    */
 423   public static int[] extractProfile(Hashtable hconsensus,
 424           boolean ignoreGapsInConsensusCalculation)
 425   {
 426     int[] rtnval = new int[74]; // 2*(5*5)+2
 427     int[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
 428     int[][] pairs = (int[][]) hconsensus
 429             .get(StructureFrequency.PAIRPROFILE);
 430
 431     if (profile == null)
 432       return null;
 433
 434     // TODO fix the object length, also do it in completeConsensus
 435     Object[] ca = new Object[625];
 436     float[] vl = new float[625];
 437     int x = 0;
 438     for (int c = 65; c < 90; c++)
 439     {
 440       for (int d = 65; d < 90; d++)
 441       {
 442         ca[x] = new int[]
 443         { c, d };
 444         vl[x] = pairs[c][d];
 445         x++;
 446       }
 447     }
 448     jalview.util.QuickSort.sort(vl, ca);
 449
 450     rtnval[0] = 2;
 451     rtnval[1] = 0;
 452     for (int c = 624; c > 0; c--)
 453     {
 454       if (vl[c] > 0)
 455       {
 456         rtnval[rtnval[0]++] = ((int[]) ca[c])[0];
 457         rtnval[rtnval[0]++] = ((int[]) ca[c])[1];
 458         rtnval[rtnval[0]] = (int) (vl[c] * 100f / profile[1][ignoreGapsInConsensusCalculation ? 1
 459                 : 0]);
 460         rtnval[1] += rtnval[rtnval[0]++];
 461       }
 462     }
 463
 464     return rtnval;
 465   }
 466
 467   public static void main(String args[])
 468   {
 469     // Short test to see if checkBpType works
 470     ArrayList<String> test = new ArrayList<String>();
 471     test.add("A");
 472     test.add("c");
 473     test.add("g");
 474     test.add("T");
 475     test.add("U");
 476     for (String i : test)
 477     {
 478       for (String j : test)
 479       {
 480         System.out.println(i + "-" + j + ": "
 481                 + StructureFrequency.checkBpType(i.charAt(0), j.charAt(0)));
 482       }
 483     }
 484   }
 485 }