src/jalview/analysis/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 }