3e007c477a31ef3d58aae966999cdc08e052900b
[jalview.git] / SecStrConsensus.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 package jalview.analysis;
22
23 import java.util.ArrayList;
24 import java.util.Hashtable;
25
26 public class SecStrConsensus
27 {
28
29   /**
30    * Internal class to represent a simple base-pair.
31    * 
32    * @author Yawn [JBPNote: ^^ is that Anne Menard or Ya(w)nn Ponty, I wonder !
33    *         ]
34    */
35   public static class SimpleBP
36   {
37     int bp5;
38
39     int bp3;
40
41     public SimpleBP()
42     {
43
44     }
45
46     public SimpleBP(int i5, int i3)
47     {
48       bp5 = i5;
49       bp3 = i3;
50     }
51
52     public void setBP5(int i5)
53     {
54       bp5 = i5;
55     }
56
57     public void setBP3(int i3)
58     {
59       bp3 = i3;
60     }
61
62     public int getBP5()
63     {
64       return bp5;
65     }
66
67     public int getBP3()
68     {
69       return bp3;
70     }
71
72     public String toString()
73     {
74       return "(" + bp5 + "," + bp3 + ")";
75     }
76
77   }
78
79   public static int[] extractConsensus(ArrayList<ArrayList<SimpleBP>> bps)
80   {
81     // We do not currently know the length of the alignment
82     // => Estimate it as the biggest index of a base-pair plus one.
83     int maxlength = 0;
84     for (ArrayList<SimpleBP> strs : bps)
85     {
86       for (SimpleBP bp : strs)
87       {
88
89         maxlength = Math.max(1 + Math.max(bp.bp5, bp.bp3), maxlength);
90
91       }
92     }
93     // Now we have a good estimate for length, allocate and initialize data
94     // to be fed to the dynamic programming procedure.
95     ArrayList<Hashtable<Integer, Double>> seq = new ArrayList<Hashtable<Integer, Double>>();
96     for (int i = 0; i < maxlength; i++)
97     {
98       seq.add(new Hashtable<Integer, Double>());
99     }
100     for (ArrayList<SimpleBP> strs : bps)
101     {
102       for (SimpleBP bp : strs)
103       {
104         int i = bp.bp5;
105         int j = bp.bp3;
106         Hashtable<Integer, Double> h = seq.get(i);
107         if (!h.containsKey(j))
108         {
109           h.put(j, 0.0);
110         }
111         h.put(j, h.get(j) + 1.);
112       }
113     }
114     // At this point, seq contains, at each position i, a hashtable which
115     // associates,
116     // to each possible end j, the number of time a base-pair (i,j) occurs in
117     // the alignment
118
119     // We can now run the dynamic programming procedure on this data
120     double[][] mat = fillMatrix(seq);
121     ArrayList<SimpleBP> res = backtrack(mat, seq);
122
123     // Convert it to an array, ie finalres[i] = j >= 0 iff a base-pair (i,j) is
124     // present
125     // in the consensus, or -1 otherwise
126     int[] finalres = new int[seq.size()];
127     for (int i = 0; i < seq.size(); i++)
128     {
129       finalres[i] = -1;
130     }
131     for (SimpleBP bp : res)
132     {
133       finalres[bp.bp5] = bp.bp3;
134       finalres[bp.bp3] = bp.bp5;
135     }
136
137     return finalres;
138   }
139
140   private static boolean canBasePair(
141           ArrayList<Hashtable<Integer, Double>> seq, int i, int k)
142   {
143     return seq.get(i).containsKey(k);
144   }
145
146   // Returns the score of a potential base-pair, ie the number of structures in
147   // which it is found.
148   private static double basePairScore(
149           ArrayList<Hashtable<Integer, Double>> seq, int i, int k)
150   {
151     return seq.get(i).get(k);
152   }
153
154   private static double[][] fillMatrix(
155           ArrayList<Hashtable<Integer, Double>> seq)
156   {
157     int n = seq.size();
158     double[][] tab = new double[n][n];
159     for (int m = 1; m <= n; m++)
160     {
161       for (int i = 0; i < n - m + 1; i++)
162       {
163         int j = i + m - 1;
164         tab[i][j] = 0;
165         if (i < j)
166         {
167           tab[i][j] = Math.max(tab[i][j], tab[i + 1][j]);
168           for (int k = i + 1; k <= j; k++)
169           {
170             if (canBasePair(seq, i, k))
171             {
172               double fact1 = 0;
173               if (k > i + 1)
174               {
175                 fact1 = tab[i + 1][k - 1];
176               }
177               double fact2 = 0;
178               if (k < j)
179               {
180                 fact2 = tab[k + 1][j];
181               }
182               tab[i][j] = Math.max(tab[i][j], basePairScore(seq, i, k)
183                       + fact1 + fact2);
184             }
185           }
186         }
187       }
188     }
189     return tab;
190   }
191
192   private static ArrayList<SimpleBP> backtrack(double[][] tab,
193           ArrayList<Hashtable<Integer, Double>> seq)
194   {
195     return backtrack(tab, seq, 0, seq.size() - 1);
196   }
197
198   private static ArrayList<SimpleBP> backtrack(double[][] tab,
199           ArrayList<Hashtable<Integer, Double>> seq, int i, int j)
200   {
201     ArrayList<SimpleBP> result = new ArrayList<SimpleBP>();
202     if (i < j)
203     {
204       ArrayList<Integer> indices = new ArrayList<Integer>();
205       indices.add(-1);
206       for (int k = i + 1; k <= j; k++)
207       {
208         indices.add(k);
209       }
210       for (int k : indices)
211       {
212         if (k == -1)
213         {
214           if (tab[i][j] == tab[i + 1][j])
215           {
216             result = backtrack(tab, seq, i + 1, j);
217           }
218         }
219         else
220         {
221           if (canBasePair(seq, i, k))
222           {
223             double fact1 = 0;
224             if (k > i + 1)
225             {
226               fact1 = tab[i + 1][k - 1];
227             }
228             double fact2 = 0;
229             if (k < j)
230             {
231               fact2 = tab[k + 1][j];
232             }
233             if (tab[i][j] == basePairScore(seq, i, k) + fact1 + fact2)
234             {
235               result = backtrack(tab, seq, i + 1, k - 1);
236               result.addAll(backtrack(tab, seq, k + 1, j));
237               result.add(new SimpleBP(i, k));
238             }
239           }
240         }
241       }
242     }
243     else if (i == j)
244     {
245
246     }
247     else
248     {
249
250     }
251     return result;
252   }
253
254 }