src/jalview/analysis/SecStrConsensus.java

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