src/jalview/analysis/StructureFrequency.java

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