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