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