src/jalview/analysis/Rna.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   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
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 /* Author: Lauren Michelle Lui
  22  * Methods are based on RALEE methods http://personalpages.manchester.ac.uk/staff/sam.griffiths-jones/software/ralee/
  23  * Additional Author: Jan Engelhart (2011) - Structure consensus and bug fixing
  24  * Additional Author: Anne Menard (2012) - Pseudoknot support and secondary structure consensus
  25  * */
  26
  27 package jalview.analysis;
  28
  29 import jalview.analysis.SecStrConsensus.SimpleBP;
  30 import jalview.datamodel.SequenceFeature;
  31 import jalview.util.MessageManager;
  32
  33 import java.util.ArrayList;
  34 import java.util.HashMap;
  35 import java.util.Hashtable;
  36 import java.util.List;
  37 import java.util.Map;
  38 import java.util.Stack;
  39
  40 public class Rna
  41 {
  42
  43   /**
  44    * Answers true if the character is a valid open pair rna secondary structure
  45    * symbol. Currently accepts A-Z, ([{<
  46    *
  47    * @param c
  48    * @return
  49    */
  50   public static boolean isOpeningParenthesis(char c)
  51   {
  52     return ('A' <= c && c <= 'Z' || c == '(' || c == '[' || c == '{' || c == '<');
  53   }
  54
  55   /**
  56    * Answers true if the string is a valid open pair rna secondary structure
  57    * symbol. Currently accepts A-Z, ([{<
  58    *
  59    * @param s
  60    * @return
  61    */
  62   public static boolean isOpeningParenthesis(String s)
  63   {
  64     return s != null && s.length() == 1
  65             && isOpeningParenthesis(s.charAt(0));
  66   }
  67
  68   /**
  69    * Answers true if the character is a valid close pair rna secondary structure
  70    * symbol. Currently accepts a-z, )]}>
  71    *
  72    * @param c
  73    * @return
  74    */
  75   public static boolean isClosingParenthesis(char c)
  76   {
  77     return ('a' <= c && c <= 'z' || c == ')' || c == ']' || c == '}' || c == '>');
  78   }
  79
  80   /**
  81    * Answers true if the string is a valid close pair rna secondary structure
  82    * symbol. Currently accepts a-z, )]}>
  83    *
  84    * @param s
  85    * @return
  86    */
  87   public static boolean isClosingParenthesis(String s)
  88   {
  89     return s != null && s.length() == 1
  90             && isClosingParenthesis(s.charAt(0));
  91   }
  92
  93   /**
  94    * Returns the matching open pair symbol for the given closing symbol.
  95    * Currently returns A-Z for a-z, or ([{< for )]}>, or the input symbol if it
  96    * is not a valid closing symbol.
  97    *
  98    * @param c
  99    * @return
 100    */
 101   public static char getMatchingOpeningParenthesis(char c)
 102   {
 103     if ('a' <= c && c <= 'z')
 104     {
 105       return (char) (c + 'A' - 'a');
 106     }
 107     switch (c)
 108     {
 109     case ')':
 110       return '(';
 111     case ']':
 112       return '[';
 113     case '}':
 114       return '{';
 115     case '>':
 116       return '<';
 117     default:
 118       return c;
 119     }
 120   }
 121
 122   /**
 123    * Based off of RALEE code ralee-get-base-pairs. Keeps track of open bracket
 124    * positions in "stack" vector. When a close bracket is reached, pair this
 125    * with the last matching element in the "stack" vector and store in "pairs"
 126    * vector. Remove last element in the "stack" vector. Continue in this manner
 127    * until the whole string is processed. Parse errors are thrown as exceptions
 128    * wrapping the error location - position of the first unmatched closing
 129    * bracket, or string length if there is an unmatched opening bracket.
 130    *
 131    * @param line
 132    *          Secondary structure line of an RNA Stockholm file
 133    * @return
 134    * @throw {@link WUSSParseException}
 135    */
 136   protected static List<SimpleBP> getSimpleBPs(CharSequence line)
 137           throws WUSSParseException
 138   {
 139     Hashtable<Character, Stack<Integer>> stacks = new Hashtable<Character, Stack<Integer>>();
 140     List<SimpleBP> pairs = new ArrayList<SimpleBP>();
 141     int i = 0;
 142     while (i < line.length())
 143     {
 144       char base = line.charAt(i);
 145
 146       if (isOpeningParenthesis(base))
 147       {
 148         if (!stacks.containsKey(base))
 149         {
 150           stacks.put(base, new Stack<Integer>());
 151         }
 152         stacks.get(base).push(i);
 153
 154       }
 155       else if (isClosingParenthesis(base))
 156       {
 157
 158         char opening = getMatchingOpeningParenthesis(base);
 159
 160         if (!stacks.containsKey(opening))
 161         {
 162           throw new WUSSParseException(MessageManager.formatMessage(
 163                   "exception.mismatched_unseen_closing_char",
 164                   new String[] { String.valueOf(base) }), i);
 165         }
 166
 167         Stack<Integer> stack = stacks.get(opening);
 168         if (stack.isEmpty())
 169         {
 170           // error whilst parsing i'th position. pass back
 171           throw new WUSSParseException(MessageManager.formatMessage(
 172                   "exception.mismatched_closing_char",
 173                   new String[] { String.valueOf(base) }), i);
 174         }
 175         int temp = stack.pop();
 176
 177         pairs.add(new SimpleBP(temp, i));
 178       }
 179       i++;
 180     }
 181     for (char opening : stacks.keySet())
 182     {
 183       Stack<Integer> stack = stacks.get(opening);
 184       if (!stack.empty())
 185       {
 186         /*
 187          * we have an unmatched opening bracket; report error as at
 188          * i (length of input string)
 189          */
 190         throw new WUSSParseException(MessageManager.formatMessage(
 191                 "exception.mismatched_opening_char",
 192                 new String[] { String.valueOf(opening),
 193                     String.valueOf(stack.pop()) }), i);
 194       }
 195     }
 196     return pairs;
 197   }
 198
 199
 200
 201
 202
 203   /**
 204    * Function to get the end position corresponding to a given start position
 205    *
 206    * @param indice
 207    *          - start position of a base pair
 208    * @return - end position of a base pair
 209    */
 210   /*
 211    * makes no sense at the moment :( public int findEnd(int indice){ //TODO:
 212    * Probably extend this to find the start to a given end? //could be done by
 213    * putting everything twice to the hash ArrayList<Integer> pair = new
 214    * ArrayList<Integer>(); return pairHash.get(indice); }
 215    */
 216
 217   /**
 218    * Figures out which helix each position belongs to and stores the helix
 219    * number in the 'featureGroup' member of a SequenceFeature Based off of RALEE
 220    * code ralee-helix-map.
 221    *
 222    * @param pairs
 223    *          Array of SequenceFeature (output from Rna.GetBasePairs)
 224    */
 225   protected static void computeHelixMap(SequenceFeature[] pairs)
 226   {
 227
 228     int helix = 0; // Number of helices/current helix
 229     int lastopen = 0; // Position of last open bracket reviewed
 230     int lastclose = 9999999; // Position of last close bracket reviewed
 231
 232     int open; // Position of an open bracket under review
 233     int close; // Position of a close bracket under review
 234     int j; // Counter
 235
 236     Map<Integer, Integer> helices = new HashMap<Integer, Integer>();
 237     // Keep track of helix number for each position
 238
 239     // Go through each base pair and assign positions a helix
 240     for (int i = 0; i < pairs.length; i++)
 241     {
 242
 243       open = pairs[i].getBegin();
 244       close = pairs[i].getEnd();
 245
 246       // System.out.println("open " + open + " close " + close);
 247       // System.out.println("lastclose " + lastclose + " lastopen " + lastopen);
 248
 249       // we're moving from right to left based on closing pair
 250       /*
 251        * catch things like <<..>>..<<..>> |
 252        */
 253       if (open > lastclose)
 254       {
 255         helix++;
 256       }
 257
 258       /*
 259        * catch things like <<..<<..>>..<<..>>>> |
 260        */
 261       j = pairs.length - 1;
 262       while (j >= 0)
 263       {
 264         int popen = pairs[j].getBegin();
 265
 266         // System.out.println("j " + j + " popen " + popen + " lastopen "
 267         // +lastopen + " open " + open);
 268         if ((popen < lastopen) && (popen > open))
 269         {
 270           if (helices.containsValue(popen)
 271                   && ((helices.get(popen)) == helix))
 272           {
 273             continue;
 274           }
 275           else
 276           {
 277             helix++;
 278             break;
 279           }
 280         }
 281
 282         j -= 1;
 283       }
 284
 285       // Put positions and helix information into the hashtable
 286       helices.put(open, helix);
 287       helices.put(close, helix);
 288
 289       // Record helix as featuregroup
 290       pairs[i].setFeatureGroup(Integer.toString(helix));
 291
 292       lastopen = open;
 293       lastclose = close;
 294     }
 295   }
 296
 297   /**
 298    * Answers true if the character is a recognised symbol for RNA secondary
 299    * structure. Currently accepts a-z, A-Z, ()[]{}<>.
 300    *
 301    * @param c
 302    * @return
 303    */
 304   public static boolean isRnaSecondaryStructureSymbol(char c)
 305   {
 306     return isOpeningParenthesis(c) || isClosingParenthesis(c);
 307   }
 308
 309   /**
 310    * Answers true if the string is a recognised symbol for RNA secondary
 311    * structure. Currently accepts a-z, A-Z, ()[]{}<>.
 312    *
 313    * @param s
 314    * @return
 315    */
 316   public static boolean isRnaSecondaryStructureSymbol(String s)
 317   {
 318     return isOpeningParenthesis(s) || isClosingParenthesis(s);
 319   }
 320
 321   /**
 322    * Translates a string to RNA secondary structure representation. Returns the
 323    * string with any non-SS characters changed to spaces. Accepted characters
 324    * are a-z, A-Z, and (){}[]<> brackets.
 325    *
 326    * @param ssString
 327    * @return
 328    */
 329   public static String getRNASecStrucState(String ssString)
 330   {
 331     if (ssString == null)
 332     {
 333       return null;
 334     }
 335     StringBuilder result = new StringBuilder(ssString.length());
 336     for (int i = 0; i < ssString.length(); i++)
 337     {
 338       char c = ssString.charAt(i);
 339       result.append(isRnaSecondaryStructureSymbol(c) ? c : " ");
 340     }
 341     return result.toString();
 342   }
 343
 344   /**
 345    * Answers true if the base-pair is either a Watson-Crick (A:T/U, C:G) or a
 346    * wobble (G:T/U) pair (either way round), else false
 347    *
 348    * @param first
 349    * @param second
 350    * @return
 351    */
 352   public static boolean isCanonicalOrWobblePair(char first, char second)
 353   {
 354     if (first > 'Z')
 355     {
 356       first -= 32;
 357     }
 358     if (second > 'Z')
 359     {
 360       second -= 32;
 361     }
 362
 363     switch (first)
 364     {
 365     case 'A':
 366       switch (second)
 367       {
 368       case 'T':
 369       case 'U':
 370         return true;
 371       }
 372       break;
 373     case 'C':
 374       switch (second)
 375       {
 376       case 'G':
 377         return true;
 378       }
 379       break;
 380     case 'T':
 381     case 'U':
 382       switch (second)
 383       {
 384       case 'A':
 385       case 'G':
 386         return true;
 387       }
 388       break;
 389     case 'G':
 390       switch (second)
 391       {
 392       case 'C':
 393       case 'T':
 394       case 'U':
 395         return true;
 396       }
 397       break;
 398     }
 399     return false;
 400   }
 401
 402   /**
 403    * Answers true if the base-pair is Watson-Crick - (A:T/U or C:G, either way
 404    * round), else false
 405    *
 406    * @param first
 407    * @param second
 408    * @return
 409    */
 410   public static boolean isCanonicalPair(char first, char second)
 411   {
 412
 413     if (first > 'Z')
 414     {
 415       first -= 32;
 416     }
 417     if (second > 'Z')
 418     {
 419       second -= 32;
 420     }
 421
 422     switch (first)
 423     {
 424     case 'A':
 425       switch (second)
 426       {
 427       case 'T':
 428       case 'U':
 429         return true;
 430       }
 431       break;
 432     case 'G':
 433       switch (second)
 434       {
 435       case 'C':
 436         return true;
 437       }
 438       break;
 439     case 'C':
 440       switch (second)
 441       {
 442       case 'G':
 443         return true;
 444       }
 445       break;
 446     case 'T':
 447     case 'U':
 448       switch (second)
 449       {
 450       case 'A':
 451         return true;
 452       }
 453       break;
 454     }
 455     return false;
 456   }
 457
 458   /**
 459    * Returns the matching close pair symbol for the given opening symbol.
 460    * Currently returns a-z for A-Z, or )]}> for ([{<, or the input symbol if it
 461    * is not a valid opening symbol.
 462    *
 463    * @param c
 464    * @return
 465    */
 466   public static char getMatchingClosingParenthesis(char c)
 467   {
 468     if ('A' <= c && c <= 'Z')
 469     {
 470       return (char) (c + 'a' - 'A');
 471     }
 472     switch (c)
 473     {
 474     case '(':
 475       return ')';
 476     case '[':
 477       return ']';
 478     case '{':
 479       return '}';
 480     case '<':
 481       return '>';
 482     default:
 483       return c;
 484     }
 485   }
 486
 487   public static SequenceFeature[] getHelixMap(CharSequence rnaAnnotation)
 488           throws WUSSParseException
 489   {
 490     List<SequenceFeature> result = new ArrayList<SequenceFeature>();
 491
 492     int helix = 0; // Number of helices/current helix
 493     int lastopen = 0; // Position of last open bracket reviewed
 494     int lastclose = 9999999; // Position of last close bracket reviewed
 495
 496     Map<Integer, Integer> helices = new HashMap<Integer, Integer>();
 497     // Keep track of helix number for each position
 498
 499     // Go through each base pair and assign positions a helix
 500     List<SimpleBP> bps = getSimpleBPs(rnaAnnotation);
 501     for (SimpleBP basePair : bps)
 502     {
 503       final int open = basePair.getBP5();
 504       final int close = basePair.getBP3();
 505
 506       // System.out.println("open " + open + " close " + close);
 507       // System.out.println("lastclose " + lastclose + " lastopen " + lastopen);
 508
 509       // we're moving from right to left based on closing pair
 510       /*
 511        * catch things like <<..>>..<<..>> |
 512        */
 513       if (open > lastclose)
 514       {
 515         helix++;
 516       }
 517
 518       /*
 519        * catch things like <<..<<..>>..<<..>>>> |
 520        */
 521       int j = bps.size() - 1;
 522       while (j >= 0)
 523       {
 524         int popen = bps.get(j).getBP5();
 525
 526         // System.out.println("j " + j + " popen " + popen + " lastopen "
 527         // +lastopen + " open " + open);
 528         if ((popen < lastopen) && (popen > open))
 529         {
 530           if (helices.containsValue(popen)
 531                   && ((helices.get(popen)) == helix))
 532           {
 533             continue;
 534           }
 535           else
 536           {
 537             helix++;
 538             break;
 539           }
 540         }
 541         j -= 1;
 542       }
 543
 544       // Put positions and helix information into the hashtable
 545       helices.put(open, helix);
 546       helices.put(close, helix);
 547
 548       // Record helix as featuregroup
 549       result.add(new SequenceFeature("RNA helix", "", "", open, close,
 550               String.valueOf(helix)));
 551
 552       lastopen = open;
 553       lastclose = close;
 554     }
 555
 556     return result.toArray(new SequenceFeature[result.size()]);
 557   }
 558 }