Changeset 3878

Timestamp:

Oct 25, 2007, 11:57:31 AM (15 years ago)

Author:

Nicklas Nordborg

Message:

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

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

File:

• trunk/src/plugins/core/net/sf/basedb/plugins/MedianRatioNormalization.java (modified) (7 diffs)

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

@@ -26 +26 @@
 package net.sf.basedb.plugins;
 
+import java.sql.SQLException;
 import java.util.ArrayList;
 import java.util.Arrays;
@@ -60 +61 @@
 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;
@@ -73 +72 @@
 
 /**
-   @author enell
+   @author enell, Nicklas
    @version 2.0
    @base.modified $Date$
@@ -93 +92 @@
     "The normalizer lets you exclude a percentage of the highest and lowest intensities as wel" +
     "l as all spots with at least one intensity below a fixed value.",
-    "1.0",
+    "2.5",
     "2006, Base 2 development team",
     null,
@@ -189 +188 @@
         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);
+        int blockGroupSize = (Integer) job.getValue(blockGroupParameter.getName());
+        float minIntensity = (Float)job.getValue(minIntensityParameter.getName());
+        float lowExclude = (Float) job.getValue(lowExcludeParameter.getName());
+        float highExclude = (Float) job.getValue(highExcludeParameter.getName());
+        Job thisJob = getCurrentJob(dc);
+        
+        BioAssaySet child = normalize(dc, source, thisJob, minIntensity, lowExclude, highExclude, blockGroupSize, progress);
         child.setName(childName);
         child.setDescription(childDescription);
-        dc.saveItem(child);
-        
-        // Batcher for inserting normalized data
-        batcher = child.getSpotBatcher();
-        
-        Select ratio = Selects.expression(Expressions.divide(Dynamic.column(VirtualColumn.channel(1)), Dynamic.column(VirtualColumn.channel(2))), "ratio");
-        Restriction intensityRestriction = Restrictions.and
-        (
-          Restrictions.gt(Dynamic.column(VirtualColumn.channel(1)), Expressions.aFloat((Float) job.getValue(minIntensityParameter.getName()))), 
-          Restrictions.gt(Dynamic.column(VirtualColumn.channel(2)), Expressions.aFloat((Float) job.getValue(minIntensityParameter.getName())))
-        );
-        int blockGroupSize = (Integer) job.getValue(blockGroupParameter.getName());
-        float lowExclude = (Float) job.getValue(lowExcludeParameter.getName());
-        float highExclude = (Float) job.getValue(highExcludeParameter.getName());
-        
-        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"));
-
-          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.joinRawData(JoinType.LEFT);
-            query.select(Dynamic.select(VirtualColumn.channel(1)));
-            query.select(Dynamic.select(VirtualColumn.channel(2)));
-            query.select(Dynamic.selectRawData("block"));
-            query.select(ratio);
-            query.restrict(intensityRestriction);
-            query.restrict(blockRestriction);
-            query.order(Orders.asc(Expressions.selected(ratio)));
-            
-            int count = (int)query.count(dc);
-            if (count <= 0) continue;
-            int lowCount = (int) Math.floor(count * lowExclude / 100);
-            int highCount = (int) Math.floor(count * highExclude / 100);
-
-            // check count > 0
-            int include = count - lowCount - highCount;
-            if (include <= 0) continue;
-            double median;
-            if (include % 2 == 0)
-            {
-              query.setFirstResult((include / 2) - 1 + lowCount);
-              query.setMaxResults(2);
-              DynamicResultIterator iter = query.iterate(dc);
-              median = Math.sqrt(iter.next().getFloat(iter.getIndex("ratio")) * iter.next().getFloat(iter.getIndex("ratio")));
-            }
-            else
-            {
-              query.setFirstResult(((include - 1) / 2) + lowCount);
-              query.setMaxResults(1);
-              DynamicResultIterator iter = query.iterate(dc);
-              median = iter.next().getFloat(iter.getIndex("ratio"));
-            }
-
-            float sqrtMedRatio = (float) Math.sqrt(median);
-            float invSqrtMedRatio = 1 / sqrtMedRatio;
-            Select newCh1 = Selects.expression(Expressions.multiply(Expressions.selected(Dynamic.select(VirtualColumn.channel(1))),
-              Expressions.aFloat(invSqrtMedRatio)), "newCh1");
-            Select newCh2 = Selects.expression(Expressions.multiply(Expressions.selected(Dynamic.select(VirtualColumn.channel(2))),
-              Expressions.aFloat(sqrtMedRatio)), "newCh2");
-
-            query = assay.getSpotData();
-            query.joinRawData(JoinType.LEFT);
-            query.select(Dynamic.select(VirtualColumn.COLUMN));
-            query.select(Dynamic.select(VirtualColumn.POSITION));
-            query.select(newCh1);
-            query.select(newCh2);
-            query.restrict(intensityRestriction);
-            query.restrict(blockRestriction);
-            query.setParameter("low", new Integer(i), Type.INT);
-            query.setParameter("high", new Integer(i + blockGroupSize), Type.INT);
-
-            DynamicResultIterator iter = query.iterate(dc);
-            int column = iter.getIndex(VirtualColumn.COLUMN.getName());
-            int position = iter.getIndex(VirtualColumn.POSITION.getName());
-            int ch1 = iter.getIndex("newCh1");
-            int ch2 = iter.getIndex("newCh2");
-            while (iter.hasNext())
-            {
-              SqlResult r = iter.next();
-              batcher.insert(r.getShort(column), r.getInt(position), r.getFloat(ch1), r.getFloat(ch2));
-              normalizedSpots++;
-            }
-            if (progress != null) progress.display((int)(normalizedSpots * 100 / numSpots), normalizedSpots + " spots normalized");
-          }
-        }
-        if (progress != null) progress.append("\n");
-        batcher.close();
         dc.commit();
+        int normalizedSpots = child.getNumSpots();
         response.setDone(normalizedSpots + " spots normalized, " + (source.getNumSpots() - normalizedSpots) + " spots removed");
       }
@@ -408 +294 @@
   // -------------------------------------------
   
+  /**
+    Normalise the source bioassay set using MEDIAN-RATIO 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 minIntensity All spots which have a lower intensity value in either channel
+      will be filtered out
+    @param lowExclude A percentage of the spots with lowest ratio that are excuded in the
+      median calculation
+    @param highExclude A percentage of the spots with highest ratio that are excluded in the
+      median calculation
+    @param blockGroupSize
+    @return The normalized bioassayset
+    @since 2.5
+  */
+  public BioAssaySet normalize(DbControl dc, BioAssaySet source, Job job, float minIntensity, float lowExclude, float highExclude, int blockGroupSize, ProgressReporter progress)
+  {
+    if (progress != null) progress.display(0, "Preparing to normalize...");
+
+    // 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);
+    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));
+    // ratio = ch1 / ch2
+    Expression ratio = Expressions.divide(ch1, ch2);
+    
+    //   Create restriction: ch1 > minIntensity and ch2 > minIntensity
+    Restriction intensityRestriction = Restrictions.and
+    (
+      Restrictions.gt(Dynamic.column(VirtualColumn.channel(1)), Expressions.aFloat(minIntensity)), 
+      Restrictions.gt(Dynamic.column(VirtualColumn.channel(2)), Expressions.aFloat(minIntensity))
+    );
+    
+    // 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(ratio));
+
+    // 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)
+          {
+            // Exlude low and high values
+            int count = toIndex - fromIndex;
+            int lowIndex = fromIndex + (int)Math.floor(count * lowExclude / 100);
+            int highIndex = toIndex - (int)Math.floor(count * highExclude / 100);
+            
+            if (highIndex > lowIndex)
+            {
+              // Find median ratio in the current block range
+              double median = median(data.subList(lowIndex, highIndex));
+              
+              // Correct all spot intensities
+              float sqrtMedRatio = (float)Math.sqrt(median);
+              float invSqrtMedRatio = 1 / sqrtMedRatio;
+              for (int j = fromIndex; j < toIndex; ++j)
+              {
+                SpotData spot = data.get(j);
+                double newCh1 = spot.ch1 * invSqrtMedRatio;
+                double newCh2 = spot.ch2 * sqrtMedRatio;
+                batcher.insert(bioassayColumn, spot.position, (float)newCh1, (float)newCh2);
+              }
+            }
+            normalizedSpots += toIndex - fromIndex;
+            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()
   {
@@ -444 +491 @@
     }
     return configureJob;
-  } 
-    
+  }
+  
+  private double median(List<SpotData> data)
+  {
+    int size = data.size();
+    int half = size / 2;
+    if (size % 2 == 0)
+    {
+      return Math.sqrt(data.get(half - 1).ratio * data.get(half).ratio);
+    }
+    else
+    {
+      return data.get(half).ratio;
+    }
+  }
+  
+  
+  private static class SpotData
+  {
+    final int position;
+    final float ch1;
+    final float ch2;
+    final double ratio;
+    final int block;
+
+    SpotData(int position, float ch1, float ch2, int block)
+    {
+      this.position = position;
+      this.ch1 = ch1;
+      this.ch2 = ch2;
+      this.block = block;
+      this.ratio = ch1 / ch2;
+    }
+  }
 }