revised RNA stem/loop processing code that copes with pseudoknots
[jalview.git] / src / jalview / analysis / Rna.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 /* Author: Lauren Michelle Lui 
19  * Methods are based on RALEE methods http://personalpages.manchester.ac.uk/staff/sam.griffiths-jones/software/ralee/
20  * Additional Author: Jan Engelhart (2011) - Structure consensus and bug fixing
21  * Additional Author: Anne Menard (2012) - Pseudoknot support and secondary structure consensus
22  * */
23
24 package jalview.analysis;
25
26
27 import java.util.ArrayList;
28 import java.util.Arrays;
29 import java.util.HashSet;
30 import java.util.Hashtable;
31 import java.util.Stack;
32 import java.util.Vector;
33
34
35 import jalview.analysis.SecStrConsensus.SimpleBP;
36 import jalview.datamodel.SequenceFeature;
37
38 public class Rna
39 {
40         
41         static Hashtable<Integer, Integer> pairHash = new Hashtable();
42         
43   private static final Character[] openingPars = {'(','[','{','<','A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'}; 
44   private static final Character[] closingPars = {')',']','}','>','a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z'}; 
45   
46   private static HashSet<Character> openingParsSet = new HashSet<Character>(Arrays.asList(openingPars)); 
47   private static HashSet<Character> closingParsSet = new HashSet<Character>(Arrays.asList(closingPars));
48   private static Hashtable<Character,Character> closingToOpening = new Hashtable<Character,Character>()
49   // Initializing final data structure
50   {
51         private static final long serialVersionUID = 1L;
52         {
53           for(int i=0;i<openingPars.length;i++)
54           {
55                   System.out.println(closingPars[i]+"->"+openingPars[i]);
56                   put(closingPars[i],openingPars[i]);             
57           }
58   }};
59         
60   private static boolean isOpeningParenthesis(char c)
61   {
62           return openingParsSet.contains(c);
63   }
64         
65   private static boolean isClosingParenthesis(char c)
66   {
67           return closingParsSet.contains(c);
68   }
69
70   private static char matchingOpeningParenthesis(char closingParenthesis) throws WUSSParseException
71   {
72           if (!isClosingParenthesis(closingParenthesis))
73           {
74                   throw new WUSSParseException("Querying matching opening parenthesis for non-closing parenthesis character "+closingParenthesis, -1);
75           }
76         
77           return closingToOpening.get(closingParenthesis);
78   }
79   
80   /**
81    * Based off of RALEE code ralee-get-base-pairs. Keeps track of open bracket
82    * positions in "stack" vector. When a close bracket is reached, pair this
83    * with the last element in the "stack" vector and store in "pairs" vector.
84    * Remove last element in the "stack" vector. Continue in this manner until
85    * the whole string is processed.
86    * 
87    * @param line
88    *          Secondary structure line of an RNA Stockholm file
89    * @return Array of SequenceFeature; type = RNA helix, begin is open base
90    *         pair, end is close base pair
91    */
92   public static Vector<SimpleBP> GetSimpleBPs(CharSequence line) throws WUSSParseException
93   {
94     System.out.println(line);
95         Hashtable<Character,Stack<Integer>> stacks = new Hashtable<Character,Stack<Integer>>();
96     Vector<SimpleBP> pairs = new Vector<SimpleBP>();
97     int i = 0;
98     while (i < line.length())
99     {
100       char base = line.charAt(i);
101      
102       if (isOpeningParenthesis(base))
103       {
104           if (!stacks.containsKey(base)){
105                   stacks.put(base, new Stack<Integer>());
106           }
107         stacks.get(base).push(i);
108        
109       }
110       else if (isClosingParenthesis(base))
111       {
112         
113           char opening = matchingOpeningParenthesis(base);
114           
115           if (!stacks.containsKey(opening)){
116                   throw new WUSSParseException("Mismatched (unseen) closing character "+base, i);
117           }
118           
119          Stack<Integer> stack = stacks.get(opening);
120         if (stack.isEmpty())
121         {
122           // error whilst parsing i'th position. pass back
123                   throw new WUSSParseException("Mismatched closing character "+base, i);
124         }
125         int temp = stack.pop();
126         
127         pairs.add(new SimpleBP(temp,i));
128       }
129       i++;
130     }
131     for(char opening: stacks.keySet())
132     {
133         Stack<Integer> stack = stacks.get(opening);
134         if (!stack.empty())
135         {
136                   throw new WUSSParseException("Mismatched opening character "+opening+" at "+stack.pop(), i);                  
137         }
138     }
139     return pairs;
140   }
141   
142   public static SequenceFeature[] GetBasePairs(CharSequence line) throws WUSSParseException
143   {
144           Vector<SimpleBP> bps = GetSimpleBPs(line);
145           SequenceFeature[] outPairs = new SequenceFeature[bps.size()];
146     for (int p = 0; p < bps.size(); p++)
147     {
148         SimpleBP bp = bps.elementAt(p);
149         outPairs[p] = new SequenceFeature("RNA helix", "", "", bp.getBP5(),bp.getBP3(), "");
150     }
151     return outPairs;
152   }
153   
154   
155   public static  ArrayList<SimpleBP> GetModeleBP(CharSequence line) throws WUSSParseException
156   {
157           Vector<SimpleBP> bps = GetSimpleBPs(line);
158           return new ArrayList<SimpleBP>(bps);
159   }
160   
161   
162   /**
163    * Function to get the end position corresponding to a given start position
164    * @param indice - start position of a base pair
165    * @return - end position of a base pair
166    */
167   /*makes no sense at the moment :(
168   public int findEnd(int indice){
169           //TODO: Probably extend this to find the start to a given end?
170           //could be done by putting everything twice to the hash
171           ArrayList<Integer> pair = new ArrayList<Integer>();
172           return pairHash.get(indice);
173   }*/
174   
175
176   /**
177    * Figures out which helix each position belongs to and stores the helix
178    * number in the 'featureGroup' member of a SequenceFeature Based off of RALEE
179    * code ralee-helix-map.
180    * 
181    * @param pairs
182    *          Array of SequenceFeature (output from Rna.GetBasePairs)
183    */
184   public static void HelixMap(SequenceFeature[] pairs)
185   {
186
187     int helix = 0; // Number of helices/current helix
188     int lastopen = 0; // Position of last open bracket reviewed
189     int lastclose = 9999999; // Position of last close bracket reviewed
190     int i = pairs.length; // Number of pairs
191
192     int open; // Position of an open bracket under review
193     int close; // Position of a close bracket under review
194     int j; // Counter
195
196     Hashtable helices = new Hashtable(); // Keep track of helix number for each
197                                          // position
198
199     // Go through each base pair and assign positions a helix
200     for (i = 0; i < pairs.length; i++)
201     {
202
203       open = pairs[i].getBegin();
204       close = pairs[i].getEnd();
205
206       // System.out.println("open " + open + " close " + close);
207       // System.out.println("lastclose " + lastclose + " lastopen " + lastopen);
208
209       // we're moving from right to left based on closing pair
210       /*
211        * catch things like <<..>>..<<..>> |
212        */
213       if (open > lastclose)
214       {
215         helix++;
216       }
217
218       /*
219        * catch things like <<..<<..>>..<<..>>>> |
220        */
221       j = pairs.length - 1;
222       while (j >= 0)
223       {
224         int popen = pairs[j].getBegin();
225
226         // System.out.println("j " + j + " popen " + popen + " lastopen "
227         // +lastopen + " open " + open);
228         if ((popen < lastopen) && (popen > open))
229         {
230           if (helices.containsValue(popen)
231                   && (((Integer) helices.get(popen)) == helix))
232           {
233             continue;
234           }
235           else
236           {
237             helix++;
238             break;
239           }
240         }
241
242         j -= 1;
243       }
244
245       // Put positions and helix information into the hashtable
246       helices.put(open, helix);
247       helices.put(close, helix);
248
249       // Record helix as featuregroup
250       pairs[i].setFeatureGroup(Integer.toString(helix));
251
252       lastopen = open;
253       lastclose = close;
254
255     }
256   }
257 }
258