Changeset 3877

Timestamp:

Oct 25, 2007, 9:57:20 AM (15 years ago)

Author:

Nicklas Nordborg

Message:

References #797: Enhance performance for LOWESS and Medin-ratio plug-ins

LOWESS has now been optimized to use fewer database queries. Initial test results indicates at
least 80% reduced execution time, but I have only tested on a 1-4 bioassays.

File:

• trunk/src/plugins/core/net/sf/basedb/plugins/LowessNormalization.java (modified) (14 diffs)

trunk/src/plugins/core/net/sf/basedb/plugins/LowessNormalization.java

@@ -42 +42 @@
 import net.sf.basedb.core.SpotBatcher;
 import net.sf.basedb.core.Transformation;
+import net.sf.basedb.core.Type;
 import net.sf.basedb.core.VirtualColumn;
 import net.sf.basedb.core.plugin.About;
@@ -51 +52 @@
 import net.sf.basedb.core.plugin.Request;
 import net.sf.basedb.core.plugin.Response;
-import net.sf.basedb.core.query.Aggregations;
 import net.sf.basedb.core.query.Dynamic;
+import net.sf.basedb.core.query.Expression;
 import net.sf.basedb.core.query.Expressions;
-import net.sf.basedb.core.query.JoinType;
 import net.sf.basedb.core.query.Orders;
 import net.sf.basedb.core.query.Restriction;
 import net.sf.basedb.core.query.Restrictions;
-import net.sf.basedb.core.query.Select;
 import net.sf.basedb.core.query.Selects;
 import net.sf.basedb.core.query.SqlResult;
 import net.sf.basedb.util.Values;
 
+import java.sql.SQLException;
 import java.util.ArrayList;
 import java.util.Arrays;
@@ -73 +73 @@
 
 /**
-   @author enell
+   @author enell, Nicklas
    @version 2.0
    @base.modified $Date$
@@ -103 +103 @@
     "method addressing single and multiple slide systematic " +
     "variation. Nucleic Acids Res 2002, 30:e15.",
-    "1.0",
+    "2.5",
     "2006, Base 2 development team",
     null,
@@ -193 +193 @@
     if (command.equals(Request.COMMAND_EXECUTE))
     {
-      SpotBatcher batcher = null;
       DbControl dc = null;
       try
@@ -201 +200 @@
         String childName = Values.getString((String)job.getValue(CHILD_NAME), source.getName());
         String childDescription = (String)job.getValue(CHILD_DESCRIPTION);
-        
-        // Create Transformation
-        Transformation t = source.newTransformation(getCurrentJob(dc));
-        t.setName(about.getName());
-        dc.saveItem(t);
-        
-        // Create the normalized bioassay set
-        BioAssaySet child = t.newProduct(null, "new", true);
-        child.setName(childName);
-        child.setDescription(childDescription);
-        dc.saveItem(child);
-        
-        // Batcher for inserting normalized data
-        batcher = child.getSpotBatcher();
-
         int blockGroupSize = (Integer) job.getValue(blockGroupParameter.getName());
         float delta = (Float) job.getValue(deltaParameter.getName());
         float fitFraction = (Float) job.getValue(fitFractionParameter.getName());
-        int iter = (Integer) job.getValue(iterParameter.getName());
-
-        Select M = Selects.expression
-        (
-          Expressions.log2(Expressions.divide(Dynamic.column(VirtualColumn.channel(1)), Dynamic.column(VirtualColumn.channel(2)))),
-          "m"
-        );
-        Select A = Selects.expression
-        (
-          Expressions.log10(Expressions.sqrt(Expressions.multiply(Dynamic.column(VirtualColumn.channel(1)), Dynamic.column(VirtualColumn.channel(2))))),
-          "a"
-        );
-        Restriction intensityRestriction = Restrictions.and
-        (
-          Restrictions.gt(Dynamic.column(VirtualColumn.channel(1)), Expressions.aFloat(0)), 
-          Restrictions.gt(Dynamic.column(VirtualColumn.channel(2)), Expressions.aFloat(0))
-        );
-
-        DynamicSpotQuery query;
-        DynamicResultIterator resultIter;
-        query = source.getSpotData();
-        query.restrict(intensityRestriction);
-        long numSpots = query.count(dc);
-        int normalizedSpots = 0;
-        if (progress != null) progress.display((int)(normalizedSpots / numSpots * 100), normalizedSpots + " spots normalized");
-        
-        List<BioAssay> assays = source.getBioAssays().list(dc);
-        for (BioAssay assay : assays)
-        {
-          query = assay.getSpotData();
-          query.restrict(intensityRestriction);
-          query.joinRawData(JoinType.LEFT);
-          query.select(Selects.expression(Aggregations.max(Dynamic.rawData("block")), "max"));
-          query.select(Selects.expression(Aggregations.min(Dynamic.rawData("block")), "min"));
-          resultIter = query.iterate(dc);
-          SqlResult result = resultIter.next();
-          int maxBlock = result.getInt(resultIter.getIndex("max"));
-          int minBlock = result.getInt(resultIter.getIndex("min"));
-          resultIter.close();
-
-          for (int i = minBlock; i <= maxBlock; i += blockGroupSize)
-          {
-            Restriction blockRestriction = Restrictions.between
-            ( 
-              Expressions.selected(Dynamic.selectRawData("block")),
-              Expressions.integer(i),
-              Expressions.integer(i+blockGroupSize)
-            );
-            
-            query = assay.getSpotData();
-            query.select(M);
-            query.select(A);
-            query.select(Dynamic.selectRawData("block"));
-            query.joinRawData(JoinType.LEFT);
-            query.restrictPermanent(intensityRestriction);
-            query.restrictPermanent(blockRestriction);
-            query.orderPermanent(Orders.asc(Expressions.selected(A)));
-  
-            int count = (int)query.count(dc);
-            if (count <= 0) continue;
-  
-            List<Double> mValues = new ArrayList<Double>(count);
-            List<Double> aValues = new ArrayList<Double>(count);
-            DynamicResultIterator it = query.iterate(dc);
-            
-            int aIndex = it.getIndex(A.getAlias());
-            int mIndex = it.getIndex(M.getAlias());
-            
-            while (it.hasNext())
-            {
-              SqlResult r = it.next();
-              aValues.add((double) r.getFloat(aIndex));
-              mValues.add((double) r.getFloat(mIndex));
-            }
-            it.close();
-            
-            List<Double> smoothCurve = lowess(aValues, mValues, fitFraction, iter, delta);
-  
-            query.reset();
-            query.select(Dynamic.select(VirtualColumn.COLUMN));
-            query.select(Dynamic.select(VirtualColumn.POSITION));
-            query.select(Dynamic.select(VirtualColumn.channel(1)));
-            query.select(Dynamic.select(VirtualColumn.channel(2)));
-            query.select(A);
-            query.joinRawData(JoinType.LEFT);
-            
-            it = query.iterate(dc);
-            int column = it.getIndex(VirtualColumn.COLUMN.getName());
-            int position = it.getIndex(VirtualColumn.POSITION.getName());
-            int ch1Col = it.getIndex(VirtualColumn.channel(1).getName());
-            int ch2Col = it.getIndex(VirtualColumn.channel(2).getName());
-            
-            for (int j = 0; j < smoothCurve.size() && it.hasNext(); j++)
-            {
-              SqlResult r = it.next();
-  
-              double factor = Math.exp(smoothCurve.get(j) * 0.5);
-              double ch1 = r.getFloat(ch1Col)/factor;
-              double ch2 = r.getFloat(ch2Col)*factor;
-              batcher.insert(r.getShort(column), r.getInt(position), (float) ch1, (float) ch2);
-            }
-            it.close();
-            
-            normalizedSpots += smoothCurve.size();
-            if (progress != null) progress.display((int)((normalizedSpots * 100) / numSpots), normalizedSpots + " spots normalized");
-          }
-        }
-        batcher.close();
+        int iterations = (Integer) job.getValue(iterParameter.getName());
+        Job thisJob = getCurrentJob(dc);
+        
+        BioAssaySet child = normalize(dc, source, thisJob, fitFraction, delta, iterations, blockGroupSize, progress);
+        child.setName(childName);
+        child.setDescription(childDescription);
         dc.commit();
+        int normalizedSpots = child.getNumSpots();
         if (progress != null) progress.display(100, normalizedSpots + " spots normalized\n");
         response.setDone(normalizedSpots + " spots normalized, " + (source.getNumSpots() - normalizedSpots) + " spots removed");
@@ -418 +301 @@
   // -------------------------------------------
 
+  /**
+    Normalise the source bioassay set using LOWESS normalization.
+    @param dc The DbControl to use for database access
+    @param source The source bioassay set that is going to be normalized
+    @param job The job that is doing the normalization, or null
+    @param fitFraction
+    @param delta
+    @param iterations
+    @param blockGroupSize
+    @return The normalized bioassayset
+    @since 2.5
+  */
+  public BioAssaySet normalize(DbControl dc, BioAssaySet source, Job job, float fitFraction, float delta, int iterations, int blockGroupSize, ProgressReporter progress)
+  {
+    if (progress != null) progress.display(0, "Preparing to normalize...");
+    
+    // Create Transformation
+    Transformation t = source.newTransformation(job);
+    t.setName(about.getName());
+    dc.saveItem(t);
+    
+    // Create the normalized bioassay set
+    BioAssaySet child = t.newProduct(null, "new", true);
+    dc.saveItem(child);
+    
+    // Batcher for inserting normalized data
+    SpotBatcher batcher = child.getSpotBatcher();
+
+    // Expressions used to get data
+    Expression block = Dynamic.rawData("block");
+    Expression ch1 = Dynamic.column(VirtualColumn.channel(1));
+    Expression ch2 = Dynamic.column(VirtualColumn.channel(2));
+    // A = log10(sqrt(ch1 * ch2))
+    Expression A = Expressions.log10(Expressions.sqrt(Expressions.multiply(ch1, ch2)));
+
+    // Create restriction: ch1 > 0 and ch2 > 0
+    Restriction intensityRestriction = Restrictions.and(
+        Restrictions.gt(ch1, Expressions.aFloat(0)), 
+        Restrictions.gt(ch2, Expressions.aFloat(0))
+      );
+
+    // Create restriction: column = :bioAssayColumn
+    Restriction bioAssayRestriction = Restrictions.eq(
+        Dynamic.column(VirtualColumn.COLUMN), 
+        Expressions.parameter("bioAssayColumn")
+      );
+      
+    // Count number of spots that is going to be normalized
+    DynamicSpotQuery query = source.getSpotData();
+    query.restrict(intensityRestriction);
+    long numSpots = query.count(dc);
+    int normalizedSpots = 0;
+    if (progress != null) progress.display((int)(normalizedSpots / numSpots * 100), normalizedSpots + " spots normalized");
+    
+    // Create query to retrieve spot data: COLUMN, POSITION, ch1, ch2, block
+    // We use a parameter to restrict the query to return data for one bioassay at a time
+    query.select(Dynamic.select(VirtualColumn.POSITION));
+    query.select(Selects.expression(ch1, "ch1"));
+    query.select(Selects.expression(ch2, "ch2"));
+    query.select(Selects.expression(block, "block"));
+    query.restrict(bioAssayRestriction);
+    query.order(Orders.asc(block));
+    query.order(Orders.asc(A));
+      
+    // Normalize one bioassay at a time
+    List<BioAssay> assays = source.getBioAssays().list(dc);
+  
+    try
+    {
+      for (BioAssay assay : assays)
+      {
+        // Prepare list for holding data
+        int assaySpots = assay.getNumSpots();
+        List<SpotData> data = new ArrayList<SpotData>(assaySpots);
+        
+        // Load spot data for this bioassay
+        short bioassayColumn = assay.getDataCubeColumnNo();
+        query.setParameter("bioAssayColumn", (int)bioassayColumn, Type.INT);
+        
+        DynamicResultIterator it = query.iterate(dc);
+        int positionIndex = it.getIndex(VirtualColumn.POSITION.getName());
+        int ch1Index = it.getIndex("ch1");
+        int ch2Index = it.getIndex("ch2");
+        int blockIndex = it.getIndex("block");
+        
+        // Copy bioassay data to SpotData objects
+        while (it.hasNext())
+        {
+          SqlResult r = it.next();
+          SpotData spot = new SpotData(r.getInt(positionIndex),
+            r.getFloat(ch1Index), r.getFloat(ch2Index), r.getInt(blockIndex));
+          data.add(spot);
+        }
+        it.close();
+        
+        // Continue with next bioassay if there is no data
+        int dataSize = data.size();
+        if (dataSize == 0) continue;
+        
+        // Get range of block numbers - NOTE! query must return spots sorted in block order
+        int minBlock = data.get(0).block;
+        int maxBlock = data.get(data.size()-1).block;
+        
+        int fromIndex = 0;
+        int toIndex = 0;
+        int fromBlock = minBlock;
+        // Normalize each block range independently: fromBlock + blockGroupSize --> toBlock
+        while (fromBlock <= maxBlock)
+        {
+          // Find start and end index for current block range
+          int toBlock = fromBlock + blockGroupSize - 1;
+          if (toBlock > maxBlock) toBlock = maxBlock;
+          fromIndex = toIndex;
+          // Data is sorted by block; find index of last spot with: block <= toBlock
+          // spot given by toIndex should not be included
+          while (toIndex < dataSize && data.get(toIndex).block <= toBlock)
+          {
+            ++toIndex;
+          }
+          
+          if (toIndex > fromIndex)
+          {
+            List<Double> smoothCurve = lowess(data.subList(fromIndex, toIndex), fitFraction, iterations, delta);
+            for (int j = 0; j < smoothCurve.size(); ++j)
+            {
+              SpotData spot = data.get(fromIndex + j);
+              double factor = Math.exp(smoothCurve.get(j) * 0.5);
+              double newCh1 = spot.ch1/factor;
+              double newCh2 = spot.ch2*factor;
+              batcher.insert(bioassayColumn, spot.position, (float) newCh1, (float) newCh2);
+            }
+            normalizedSpots += smoothCurve.size();
+            if (progress != null) progress.display((int)((normalizedSpots * 100) / numSpots), normalizedSpots + " spots normalized");
+          }
+          fromBlock = toBlock + 1;
+        }
+      }
+      batcher.flush();
+      batcher.close();
+    }
+    catch (SQLException e)
+    {
+      throw new BaseException(e);
+    }
+    return child;
+  }
+  
   private RequestInformation getConfigureJobParameters()
   {
@@ -455 +485 @@
     }
     return configureJob;
-  } 
-  
-
-  private static List<Double> lowess(List<Double> x, List<Double> y, double f, int iter, double delta)
-  {
-    Double[] smoothCurve = new Double[x.size()];
-    int windowSize = Math.min(x.size(), (int) (x.size() * f + 0.5));
+  }
+
+  private static List<Double> lowess(List<SpotData> data, double f, int iter, double delta)
+  {
+    int dataSize = data.size();
+    Double[] smoothCurve = new Double[dataSize];
+    int windowSize = Math.min(dataSize, (int) (dataSize * f + 0.5));
 
     List<Double> wFit = new ArrayList<Double>();
-    wFit.addAll(Collections.nCopies(x.size(), 1D));
-    for (int iteration = 0; iteration < iter; iteration++)
+    wFit.addAll(Collections.nCopies(dataSize, 1D));
+    for (int iteration = 0; iteration < iter; ++iteration)
     {
       int windowStart = 0;
       int i = 0;
       int prevI = -1;
-      while (prevI < x.size() - 1)
-      {
-        double xi = x.get(i);
-        // center window around the i:th value in x
+      while (prevI < dataSize - 1)
+      {
+        double Ai = data.get(i).A;
+        // center window around the i:th value in A
         // while distance from windowStart to i is greater then distance from i to windowEnd: move window
-        while (windowStart + windowSize < x.size() && xi - x.get(windowStart) > x.get(windowStart + windowSize) - xi)
+        while ((windowStart + windowSize < dataSize) && 
+            (Ai - data.get(windowStart).A > data.get(windowStart + windowSize).A - Ai))
         {
           windowStart++;
         }
         
-        List<Double> xWindow = x.subList(windowStart, windowStart + windowSize);
-        List<Double> yWindow = y.subList(windowStart, windowStart + windowSize);
+        List<SpotData> window = data.subList(windowStart, windowStart + windowSize);
         List<Double> wFitWindow = wFit.subList(windowStart, windowStart + windowSize);
         
-        List<Double> w = calculateWeights(xWindow, xi, wFitWindow);
-        double[] km = weightedLeastSquaresRegression(xWindow, yWindow, w);
-        smoothCurve[i] = km[0] * xi + km[1];
+        List<Double> w = calculateWeights(window, Ai, wFitWindow);
+        double[] km = weightedLeastSquaresRegression(window, w);
+        smoothCurve[i] = km[0] * Ai + km[1];
 
         // Interpolate skipped points due to delta
         if (prevI + 1 < i)
         {
-          double d = xi - x.get(prevI);
+          double d = Ai - data.get(prevI).A;
           if (d == 0)
           {
@@ -503 +533 @@
             for (int j = prevI + 1; j < i; j++)
             {
-              double a = (x.get(j) - x.get(prevI)) / d;
-              smoothCurve[j] = a * smoothCurve[i] + (1D - a) * smoothCurve[prevI];
+              double t = (data.get(j).A - data.get(prevI).A) / d;
+              smoothCurve[j] = t * smoothCurve[i] + (1D - t) * smoothCurve[prevI];
             }
           }
@@ -511 +541 @@
         // increase i, next x value must be at least delta greater then current
         prevI = i;
-        double cut = xi + delta;
-        while (i < x.size()-1 && x.get(i) <= cut)
+        double cut = Ai + delta;
+        while (i < dataSize-1 && data.get(i).A <= cut)
         {
           i++;
@@ -518 +548 @@
       }
       
-      double invYWRange = 1D/(ywRangeFactor  * medianCorrection(y, smoothCurve));
+      double invYWRange = 1D/(ywRangeFactor  * medianCorrection(data, smoothCurve));
       for (int j = 0; j < wFit.size(); j++)
       {
-        double w = Math.abs((smoothCurve[j] - y.get(j)) * invYWRange);
+        double w = Math.abs((smoothCurve[j] - data.get(j).M) * invYWRange);
         wFit.set(j, w < 1 ? Math.pow(Math.pow(1 - w, 2), 2) : 0);
       }
     }
-
     return Arrays.asList(smoothCurve);
   }
 
-  private static double medianCorrection(List<Double> y, Double[] smoothCurve)
-  {
-    List<Double> temp = new ArrayList<Double>(y.size());
-    for (int i = 0; i < y.size(); i++)
-    {
-      temp.add(Math.abs(y.get(i) - smoothCurve[i]));
+  private static double medianCorrection(List<SpotData> data, Double[] smoothCurve)
+  {
+    List<Double> temp = new ArrayList<Double>(data.size());
+    for (int i = 0; i < data.size(); i++)
+    {
+      temp.add(Math.abs(data.get(i).M - smoothCurve[i]));
     }
     Collections.sort(temp);
@@ -550 +579 @@
 
   
-  private static double[] weightedLeastSquaresRegression(List<Double> x, List<Double> y, List<Double> w)
+  private static double[] weightedLeastSquaresRegression(List<SpotData> data, List<Double> w)
   {
     double k;
     double m;
-    double sumX = 0;
-    double sumY = 0;
-    double sumXX = 0;
-    double sumXY = 0;
+    double sumA = 0;
+    double sumM = 0;
+    double sumAA = 0;
+    double sumAM = 0;
     double sumW = 0;
-    for (int j = 0; j < x.size(); j++)
+    for (int j = 0; j < data.size(); j++)
     {
       double localW = w.get(j);
-      double localX = x.get(j);
-      double localY = y.get(j);
-      sumX += localX * localW;
-      sumY += localY * localW;
-      sumXX += localX * localX * localW;
-      sumXY += localX * localY * localW;
+      double localA = data.get(j).A;
+      double localM = data.get(j).M;
+      sumA += localA * localW;
+      sumM += localM * localW;
+      sumAA += localA * localA * localW;
+      sumAM += localA * localM * localW;
       sumW += localW;
     }
@@ -574 +603 @@
       throw new BaseException("Sum of weigths in line_fit is not positive");
     }
-    double denom = sumW * sumXX - sumX * sumX;
+    double denom = sumW * sumAA - sumA * sumA;
     if (denom != 0.)
     {
-      k = (sumW * sumXY - sumX * sumY) / denom;
-      m = (sumY - k * sumX) / sumW;
+      k = (sumW * sumAM - sumA * sumM) / denom;
+      m = (sumM - k * sumA) / sumW;
     }
     else
     {
       k = 0.;
-      m = sumY / sumW;
+      m = sumM / sumW;
     }
 
@@ -589 +618 @@
   }
 
-  private static List<Double> calculateWeights(List<Double> values, double x1, List<Double> wFit)
-  {
-    List<Double> w = new ArrayList<Double>(values.size());
+  private static List<Double> calculateWeights(List<SpotData> data, double A1, List<Double> wFit)
+  {
+    List<Double> w = new ArrayList<Double>(data.size());
     double invRadius = 0;
-    for (double x2 : values)
-    {
-      double abs = Math.abs(x1 - x2);
+    for (SpotData spot : data)
+    {
+      double abs = Math.abs(A1 - spot.A);
       if (abs > invRadius) invRadius = abs;
     }
     invRadius = 1 / invRadius;
-    for (int i = 0; i < values.size(); i++)
-    {
-      double x2 = values.get(i);
-      double distance = Math.abs(x1 - x2) * invRadius;
+    for (int i = 0; i < data.size(); i++)
+    {
+      double A2 = data.get(i).A;
+      double distance = Math.abs(A1 - A2) * invRadius;
       w.add((distance < 1 ? Math.pow(1D - Math.pow(Math.abs(distance), 3), 3) : 0) * wFit.get(i));
     }
     return w;
   }
-  
+
+  private static class SpotData
+  {
+    final int position;
+    final float ch1;
+    final float ch2;
+    final double M;
+    final double A;
+    final int block;
+    static final double LN2 = Math.log(2);
+
+    public SpotData(int position, float ch1, float ch2, int block)
+    {
+      this.position = position;
+      this.ch1 = ch1;
+      this.ch2 = ch2;
+      this.block = block;
+      this.M = Math.log(ch1 / ch2) / LN2;
+      this.A = Math.log10(Math.sqrt(ch1 * ch2));
+    }
+  }
+
 }