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   /**
 263    * Method to check if a base-pair is a canonical or a wobble bp
 264    *
 265    * @param up
 266    *          5' base
 267    * @param down
 268    *          3' base
 269    * @return True if it is a canonical/wobble bp
 270    */
 271   public static boolean checkBpType(char up, char down)
 272   {
 273     if (up > 'Z')
 274     {
 275       up -= 32;
 276     }
 277     if (down > 'Z')
 278     {
 279       down -= 32;
 280     }
 281
 282     switch (up)
 283     {
 284     case 'A':
 285       switch (down)
 286       {
 287       case 'T':
 288         return true;
 289       case 'U':
 290         return true;
 291       }
 292       break;
 293     case 'C':
 294       switch (down)
 295       {
 296       case 'G':
 297         return true;
 298       }
 299       break;
 300     case 'T':
 301       switch (down)
 302       {
 303       case 'A':
 304         return true;
 305       case 'G':
 306         return true;
 307       }
 308       break;
 309     case 'G':
 310       switch (down)
 311       {
 312       case 'C':
 313         return true;
 314       case 'T':
 315         return true;
 316       case 'U':
 317         return true;
 318       }
 319       break;
 320     case 'U':
 321       switch (down)
 322       {
 323       case 'A':
 324         return true;
 325       case 'G':
 326         return true;
 327       }
 328       break;
 329     }
 330     return false;
 331   }
 332
 333   /**
 334    * Compute all or part of the annotation row from the given consensus
 335    * hashtable
 336    *
 337    * @param consensus
 338    *          - pre-allocated annotation row
 339    * @param hconsensus
 340    * @param iStart
 341    * @param width
 342    * @param ignoreGapsInConsensusCalculation
 343    * @param includeAllConsSymbols
 344    */
 345   public static void completeConsensus(AlignmentAnnotation consensus,
 346           Hashtable[] hconsensus, int iStart, int width,
 347           boolean ignoreGapsInConsensusCalculation,
 348           boolean includeAllConsSymbols)
 349   {
 350     float tval, value;
 351     if (consensus == null || consensus.annotations == null
 352             || consensus.annotations.length < width)
 353     {
 354       // called with a bad alignment annotation row - wait for it to be
 355       // initialised properly
 356       return;
 357     }
 358     for (int i = iStart; i < width; i++)
 359     {
 360       Hashtable hci;
 361       if (i >= hconsensus.length || ((hci = hconsensus[i]) == null))
 362       {
 363         // happens if sequences calculated over were shorter than alignment
 364         // width
 365         consensus.annotations[i] = null;
 366         continue;
 367       }
 368       value = 0;
 369       Float fv;
 370       if (ignoreGapsInConsensusCalculation)
 371       {
 372         fv = (Float) hci.get(StructureFrequency.PID_NOGAPS);
 373       }
 374       else
 375       {
 376         fv = (Float) hci.get(StructureFrequency.PID_GAPS);
 377       }
 378       if (fv == null)
 379       {
 380         consensus.annotations[i] = null;
 381         // data has changed below us .. give up and
 382         continue;
 383       }
 384       value = fv.floatValue();
 385       String maxRes = hci.get(StructureFrequency.MAXRESIDUE).toString();
 386       String mouseOver = hci.get(StructureFrequency.MAXRESIDUE) + " ";
 387       if (maxRes.length() > 1)
 388       {
 389         mouseOver = "[" + maxRes + "] ";
 390         maxRes = "+";
 391       }
 392       int[][] profile = (int[][]) hci.get(StructureFrequency.PROFILE);
 393       int[][] pairs = (int[][]) hci.get(StructureFrequency.PAIRPROFILE);
 394
 395       if (pairs != null && includeAllConsSymbols) // Just responsible for the
 396       // tooltip
 397       // TODO Update tooltips for Structure row
 398       {
 399         mouseOver = "";
 400
 401         /*
 402          * TODO It's not sure what is the purpose of the alphabet and wheter it
 403          * is useful for structure?
 404          *
 405          * if (alphabet != null) { for (int c = 0; c < alphabet.length; c++) {
 406          * tval = ((float) profile[0][alphabet[c]]) 100f / (float)
 407          * profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; mouseOver +=
 408          * ((c == 0) ? "" : "; ") + alphabet[c] + " " + ((int) tval) + "%"; } }
 409          * else {
 410          */
 411         Object[] ca = new Object[625];
 412         float[] vl = new float[625];
 413         int x = 0;
 414         for (int c = 65; c < 90; c++)
 415         {
 416           for (int d = 65; d < 90; d++)
 417           {
 418             ca[x] = new int[]
 419             { c, d };
 420             vl[x] = pairs[c][d];
 421             x++;
 422           }
 423         }
 424         jalview.util.QuickSort.sort(vl, ca);
 425         int p = 0;
 426
 427         for (int c = 624; c > 0; c--)
 428         {
 429           if (vl[c] > 0)
 430           {
 431             tval = (vl[c] * 100f / profile[1][ignoreGapsInConsensusCalculation ? 1
 432                     : 0]);
 433             mouseOver += ((p == 0) ? "" : "; ") + (char) ((int[]) ca[c])[0]
 434                     + (char) ((int[]) ca[c])[1] + " " + ((int) tval) + "%";
 435             p++;
 436
 437           }
 438         }
 439
 440         // }
 441       }
 442       else
 443       {
 444         mouseOver += ((int) value + "%");
 445       }
 446       consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ',
 447               value);
 448     }
 449   }
 450
 451   /**
 452    * get the sorted base-pair profile for the given position of the consensus
 453    *
 454    * @param hconsensus
 455    * @return profile of the given column
 456    */
 457   public static int[] extractProfile(Hashtable hconsensus,
 458           boolean ignoreGapsInConsensusCalculation)
 459   {
 460     int[] rtnval = new int[74]; // 2*(5*5)+2
 461     int[][] profile = (int[][]) hconsensus.get(StructureFrequency.PROFILE);
 462     int[][] pairs = (int[][]) hconsensus
 463             .get(StructureFrequency.PAIRPROFILE);
 464
 465     if (profile == null)
 466       return null;
 467
 468     // TODO fix the object length, also do it in completeConsensus
 469     Object[] ca = new Object[625];
 470     float[] vl = new float[625];
 471     int x = 0;
 472     for (int c = 65; c < 90; c++)
 473     {
 474       for (int d = 65; d < 90; d++)
 475       {
 476         ca[x] = new int[]
 477         { c, d };
 478         vl[x] = pairs[c][d];
 479         x++;
 480       }
 481     }
 482     jalview.util.QuickSort.sort(vl, ca);
 483
 484     rtnval[0] = 2;
 485     rtnval[1] = 0;
 486     for (int c = 624; c > 0; c--)
 487     {
 488       if (vl[c] > 0)
 489       {
 490         rtnval[rtnval[0]++] = ((int[]) ca[c])[0];
 491         rtnval[rtnval[0]++] = ((int[]) ca[c])[1];
 492         rtnval[rtnval[0]] = (int) (vl[c] * 100f / profile[1][ignoreGapsInConsensusCalculation ? 1
 493                 : 0]);
 494         rtnval[1] += rtnval[rtnval[0]++];
 495       }
 496     }
 497
 498     return rtnval;
 499   }
 500
 501   public static void main(String args[])
 502   {
 503     // Short test to see if checkBpType works
 504     ArrayList<String> test = new ArrayList<String>();
 505     test.add("A");
 506     test.add("c");
 507     test.add("g");
 508     test.add("T");
 509     test.add("U");
 510     for (String i : test)
 511     {
 512       for (String j : test)
 513       {
 514         System.out.println(i + "-" + j + ": "
 515                 + StructureFrequency.checkBpType(i.charAt(0), j.charAt(0)));
 516       }
 517     }
 518   }
 519 }