src/jalview/analysis/SecStrConsensus.java

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