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