From 59a81cfb8e64c68d882182601663d82016f1a866 Mon Sep 17 00:00:00 2001
From: jprocter <Jim Procter>
Date: Fri, 8 Apr 2005 12:47:56 +0000
Subject: [PATCH] Alignment quality with funky coloured histograms.

---
 src/jalview/analysis/Conservation.java |  191 ++++++++++++++++++++++++++++++--
 src/jalview/datamodel/Annotation.java  |    4 +
 src/jalview/gui/AlignViewport.java     |   24 +++-
 src/jalview/gui/AnnotationPanel.java   |    2 +-
 4 files changed, 206 insertions(+), 15 deletions(-)

diff --git a/src/jalview/analysis/Conservation.java b/src/jalview/analysis/Conservation.java
index c5ddd0f..38047f8 100755
--- a/src/jalview/analysis/Conservation.java
+++ b/src/jalview/analysis/Conservation.java
@@ -27,14 +27,52 @@ public class Conservation {
   Vector sequences;
   int    start;
   int    end;
-
-  Vector total = new Vector();
-
-  String consString = "";
-
-  Sequence consSequence;
-  Hashtable        propHash;
-  int              threshold;
+  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 = "";
 
@@ -45,6 +83,8 @@ public class Conservation {
     this.sequences = sequences;
     this.start     = start;
     this.end       = end;
+    seqNums = new Vector(sequences.size());
+    calcSeqNums();
   }
 
 
@@ -189,4 +229,139 @@ public class Conservation {
     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]]++;
+        }
+        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;
+
+    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+1)<lengths[size])
+            ? ((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);
+  }
+
+
 }
diff --git a/src/jalview/datamodel/Annotation.java b/src/jalview/datamodel/Annotation.java
index 6179cfa..16ede50 100755
--- a/src/jalview/datamodel/Annotation.java
+++ b/src/jalview/datamodel/Annotation.java
@@ -18,4 +18,8 @@ public class Annotation
     secondaryStructure = ss;
     value = val;
   }
+  public Annotation(String displayChar, String desc, char ss, float val, Color colour) {
+    this(displayChar, desc, ss, val);
+    this.colour = colour;
+  }
 }
diff --git a/src/jalview/gui/AlignViewport.java b/src/jalview/gui/AlignViewport.java
index 9d63cde..43f6aff 100755
--- a/src/jalview/gui/AlignViewport.java
+++ b/src/jalview/gui/AlignViewport.java
@@ -83,29 +83,41 @@ public class AlignViewport
         alignment.getWidth());
     cons.calculate();
     cons.verdict(false, 100);
-
+    cons.findQuality();
     String sequence = cons.getConsSequence().getSequence();
+    float minR,minG,minB, maxR,maxG,maxB;
+    minR = 0.3f;
+    minG = 0f;
+    minB = 0f;
+    maxR = 1.0f-minR; maxG=0.9f-minG; maxB=0f-minB; // scalable range for colouring
+    float min = cons.qualityRange[0].floatValue();
+    float max = cons.qualityRange[1].floatValue();
+
     for (int i = 0; i < alignment.getWidth(); i++)
     {
-      int value = 0;
+      float value = 0;
       try
       {
-        value = Integer.parseInt(sequence.charAt(i) + "");
+        //value = Integer.parseInt(sequence.charAt(i) + "");
+        value = ((Double) cons.quality.get(i)).floatValue();
       }
       catch (Exception ex)
       {
         if (sequence.charAt(i) == '*') value = 10;
       }
-
+      float vprop = value-min;
+      vprop/=max;
       annotations[i] = new Annotation(sequence.charAt(i) + "",
-                                      "Conservation graph", ' ', value);
+                                      "Conservation graph", ' ', value, new Color(minR+maxR*vprop, minB+maxB*vprop, minB+maxB*vprop));
     }
 
     if(conservation==null)
     {
       conservation = new AlignmentAnnotation("Conservation",
                                              "Conservation of total alignment",
-                                             annotations, 0, 10, 1);
+                                             annotations,
+                                             cons.qualityRange[0].floatValue(),
+                                             cons.qualityRange[1].floatValue(), 1);
       alignment.addAnnotation(conservation);
     }
     else
diff --git a/src/jalview/gui/AnnotationPanel.java b/src/jalview/gui/AnnotationPanel.java
index 48f9c20..76559bd 100755
--- a/src/jalview/gui/AnnotationPanel.java
+++ b/src/jalview/gui/AnnotationPanel.java
@@ -516,7 +516,7 @@ public class AnnotationPanel extends JPanel implements MouseListener, MouseMotio
            height = (int)( (total / row.graphMax) *GRAPH_HEIGHT);
 
          }
-
+         g.setColor(row.annotations[j].colour);
          g.fillRect(x, y-height, av.charWidth, height );
        }
 
-- 
1.7.10.2