X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FConservation.java;h=97863625646b5a901a538aa5787f424972b4afa4;hb=9218eaf932c0847b83e9886f5c06d70a0bc7808e;hp=520036edf8cbaf2566b2f500760a3012c23f3b29;hpb=ff75356cd448ad8bf24733c328054a53e87c242c;p=jalview.git

diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java
index 520036e..9786362 100755
--- a/src/jalview/analysis/Conservation.java
+++ b/src/jalview/analysis/Conservation.java
@@ -27,26 +27,64 @@ public class Conservation {
   Vector sequences;
   int    start;
   int    end;
-
-  Vector total = new Vector();
-
-  String consString = "";
-
-  DrawableSequence consSequence;
-  Hashtable        propHash;
-  int              threshold;
-  Hashtable[]      freqs;
+  Vector seqNums; // vector of int vectors where first is sequence checksum
+  int maxLength=0; //  used by quality calcs
+  boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
+  private void calcSeqNums() {
+    for (int i=0; i<sequences.size(); i++) {
+      calcSeqNum(i);
+    }
+  }
+  private void calcSeqNum(int i) {
+    String sq=null; // for dumb jbuilder not-inited exception warning
+    int[] sqnum = null;
+    if (i>-1 && i<sequences.size()) {
+      sq = ((SequenceI)sequences.elementAt(i)).getSequence();
+      if (seqNums.size()<=i) {
+        seqNums.add(new int[sq.length()+1]);
+      }
+      if (sq.hashCode()!=((int[])seqNums.elementAt(i))[0])
+       {
+         int j, len;
+         seqNumsChanged = true;
+         sq = ((SequenceI) sequences.elementAt(i)).getSequence();
+         len = sq.length();
+         if (maxLength < len)
+           maxLength = len;
+         sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
+         sqnum[0] = sq.hashCode();
+         for (j = 1; j <= len; j++)
+         {
+           sqnum[j] = ( (Integer) jalview.schemes.ResidueProperties.aaHash.get(new
+               String(sq.substring(j - 1, j)))).intValue(); // yuk
+         }
+         seqNums.set(i, sqnum);
+       }
+     } else {
+       // JBPNote INFO level debug
+       System.out.println("calcSeqNum called with out of range sequence index for Alignment\n");
+     }
+   }
+   Vector total = new Vector();
+   public Vector quality;
+   public Double[] qualityRange = new Double[2];
+   String consString = "";
+
+   Sequence consSequence;
+   Hashtable        propHash;
+   int              threshold;
 
   String name = "";
 
-  public Conservation(String name,Hashtable[] freqs,Hashtable propHash, int threshold, Vector sequences, int start, int end) {
+  public Conservation(String name,Hashtable propHash, int threshold, Vector sequences, int start, int end) {
     this.name      = name;
-    this.freqs     = freqs;
     this.propHash  = propHash;
     this.threshold = threshold;
     this.sequences = sequences;
     this.start     = start;
     this.end       = end;
+    seqNums = new Vector(sequences.size());
+    calcSeqNums();
   }
 
 
@@ -121,17 +159,19 @@ public class Conservation {
     }
   }
 
-  public int countGaps(int j) {
+  public int countGaps(int j)
+  {
     int count = 0;
 
-    for (int i = 0; i < sequences.size();i++) {
+    for (int i = 0; i < sequences.size();i++)
+    {
       if( j+1 > ((Sequence)sequences.elementAt(i)).getSequence().length())
       {  count++; continue;}
 
-      String tmp = ((Sequence)sequences.elementAt(i)).getSequence().substring(j,j+1);
-      if (tmp.equals(" ") || tmp.equals(".") || tmp.equals("-")) {
+      char c = ((Sequence)sequences.elementAt(i)).getSequence().charAt(j);
+      if (jalview.util.Comparison.isGap((c)))
         count++;
-      }
+
     }
     return count;
   }
@@ -182,11 +222,149 @@ public class Conservation {
       }
     }
 
-    consSequence = new DrawableSequence(name,consString,start,end);
+    consSequence = new Sequence(name,consString,start,end);
   }
 
-  public jalview.gui.DrawableSequence getConsSequence() {
+  public Sequence getConsSequence() {
     return consSequence;
   }
 
+  // From Alignment.java in jalview118
+
+  public void findQuality() {
+    findQuality(0,maxLength-1);
+  }
+
+  int[][] cons2;
+
+  private void percentIdentity2() {
+    calcSeqNums(); // updates maxLength, too.
+    if (cons2==null || seqNumsChanged) {
+      cons2 = new int[maxLength][24];
+      // Initialize the array
+      for (int j=0;j<24;j++) {
+        for (int i=0; i < maxLength;i++) {
+          cons2[i][j] = 0;
+        }
+      }
+
+      int sqnum[];
+      int j = 0;
+      while(j < sequences.size()) {
+        sqnum=(int[])seqNums.elementAt(j);
+        for (int i = 1; i < sqnum.length; i++) {
+          cons2[i-1][sqnum[i]]++;
+        }
+        for (int i=sqnum.length-1; i<maxLength; i++) {
+          cons2[i][23]++; // gap count
+        }
+        j++;
+      }
+
+      // unnecessary ?
+      /* for (int i=start; i <= end; i++) {
+           int max = -1000;
+    int maxi = -1;
+    int maxj = -1;
+
+    for (int j=0;j<24;j++) {
+        if (cons2[i][j] > max) {
+        max = cons2[i][j];
+        maxi = i;
+        maxj = j;
+      }
+
+    }
+  } */
+    }
+
+}
+
+
+public void findQuality(int start, int end) {
+    quality = new Vector();
+    double max = -10000;
+    String s = "";
+    int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
+    //Loop over columns // JBPNote Profiling info
+    //    long ts = System.currentTimeMillis();
+    //long te = System.currentTimeMillis();
+    percentIdentity2();
+
+    int size = seqNums.size();
+    int[] lengths = new int[size];
+
+    for (int l = 0; l < size; l++)
+      lengths[l] = ((int[]) seqNums.get(l)).length-1;
+
+    for (int j=start; j <= end; j++) {
+      double bigtot = 0;
+
+      // First Xr = depends on column only
+      double x[] = new double[24];
+
+      for (int ii=0; ii < 24; ii++) {
+	x[ii] = 0;
+	try {
+	  for (int i2=0; i2 < 24; i2++) {
+	    x[ii]  += (double)cons2[j][i2] * BLOSUM62[ii][i2]+4;
+	  }
+	} catch (Exception e) {
+	  System.out.println("Exception : "  + e);
+	}
+	//System.out.println("X " + ii + " " + x[ii]);
+	x[ii] /= (size);
+	//System.out.println("X " + ii + " " + x[ii]);
+      }
+      // Now calculate D for each position and sum
+      for (int k=0; k < size; k++) {
+	double tot = 0;
+	double[] xx = new double[24];
+        int seqNum =
+            (j<lengths[k])
+            ? ((int[]) seqNums.get(k))[j+1]
+            : 23; // Sequence, or gap at the end
+
+        // This is a loop over r
+        for (int i=0; i < 23; i++) {
+          double sr = 0;
+          try {
+            sr = (double)BLOSUM62[i][seqNum]+4;
+          } catch (Exception e) {
+            System.out.println("Exception in sr " + e);
+          }
+          //Calculate X with another loop over residues
+
+          //  System.out.println("Xi " + i + " " + x[i] + " " + sr);
+          xx[i] = x[i] - sr;
+
+          tot += xx[i]*xx[i];
+        }
+        bigtot += Math.sqrt(tot);
+      }
+      // This is the quality for one column
+      if (max < bigtot) {max = bigtot;}
+      //      bigtot  = bigtot * (size-cons2[j][23])/size;
+
+      quality.addElement(new Double(bigtot));
+
+      s += "-";
+
+      // Need to normalize by gaps
+    }
+    double newmax=-10000;
+    for (int j=start; j <= end; j++) {
+     double tmp =  ((Double)quality.elementAt(j)).doubleValue();
+     tmp = (max - tmp)*(size-cons2[j][23])/size;
+     //     System.out.println(tmp+ " " + j);
+     quality.setElementAt(new Double(tmp),j);
+     if (tmp>newmax)
+       newmax = tmp;
+    }
+    //    System.out.println("Quality " + s);
+    qualityRange[0] = new Double(0);
+    qualityRange[1] = new Double(newmax);
+  }
+
+
 }