src/jalview/analysis/StructureFrequency.java

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