src/core/net/sf/basedb/util/IntensityCalculatorUtil.java

Code
Comments
Other
Rev Date Author Line
1575 07 Nov 05 nicklas 1 /*
1575 07 Nov 05 nicklas 2   $Id$
1575 07 Nov 05 nicklas 3
3675 16 Aug 07 jari 4   Copyright (C) 2005 Nicklas Nordborg
4889 06 Apr 09 nicklas 5   Copyright (C) 2006 Jari Häkkinen, Nicklas Nordborg, Gregory Vincic
1575 07 Nov 05 nicklas 6   
2304 22 May 06 jari 7   This file is part of BASE - BioArray Software Environment.
2304 22 May 06 jari 8   Available at http://base.thep.lu.se/
1575 07 Nov 05 nicklas 9
1575 07 Nov 05 nicklas 10   BASE is free software; you can redistribute it and/or
1575 07 Nov 05 nicklas 11   modify it under the terms of the GNU General Public License
4479 05 Sep 08 jari 12   as published by the Free Software Foundation; either version 3
1575 07 Nov 05 nicklas 13   of the License, or (at your option) any later version.
1575 07 Nov 05 nicklas 14
1575 07 Nov 05 nicklas 15   BASE is distributed in the hope that it will be useful,
1575 07 Nov 05 nicklas 16   but WITHOUT ANY WARRANTY; without even the implied warranty of
1575 07 Nov 05 nicklas 17   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
1575 07 Nov 05 nicklas 18   GNU General Public License for more details.
1575 07 Nov 05 nicklas 19
1575 07 Nov 05 nicklas 20   You should have received a copy of the GNU General Public License
4515 11 Sep 08 jari 21   along with BASE. If not, see <http://www.gnu.org/licenses/>.
1575 07 Nov 05 nicklas 22 */
1575 07 Nov 05 nicklas 23 package net.sf.basedb.util;
1575 07 Nov 05 nicklas 24
2497 08 Aug 06 nicklas 25 import net.sf.basedb.core.DataQuery;
1575 07 Nov 05 nicklas 26 import net.sf.basedb.core.DbControl;
1575 07 Nov 05 nicklas 27 import net.sf.basedb.core.RawDataType;
1584 10 Nov 05 nicklas 28 import net.sf.basedb.core.ArrayDesign;
1575 07 Nov 05 nicklas 29 import net.sf.basedb.core.Experiment;
1575 07 Nov 05 nicklas 30 import net.sf.basedb.core.RawBioAssay;
1575 07 Nov 05 nicklas 31 import net.sf.basedb.core.Job;
1575 07 Nov 05 nicklas 32 import net.sf.basedb.core.Transformation;
1575 07 Nov 05 nicklas 33 import net.sf.basedb.core.BioAssaySet;
1575 07 Nov 05 nicklas 34 import net.sf.basedb.core.BioAssay;
1575 07 Nov 05 nicklas 35 import net.sf.basedb.core.PositionBatcher;
1575 07 Nov 05 nicklas 36 import net.sf.basedb.core.SpotBatcher;
1575 07 Nov 05 nicklas 37 import net.sf.basedb.core.MappingBatcher;
1577 08 Nov 05 nicklas 38 import net.sf.basedb.core.VirtualColumn;
1575 07 Nov 05 nicklas 39 import net.sf.basedb.core.DataResultIterator;
1576 08 Nov 05 nicklas 40 import net.sf.basedb.core.DynamicRawDataQuery;
1576 08 Nov 05 nicklas 41 import net.sf.basedb.core.DynamicResultIterator;
1873 31 Jan 06 nicklas 42 import net.sf.basedb.core.ProgressReporter;
1575 07 Nov 05 nicklas 43 import net.sf.basedb.core.BaseException;
1575 07 Nov 05 nicklas 44 import net.sf.basedb.core.InvalidDataException;
1575 07 Nov 05 nicklas 45 import net.sf.basedb.core.InvalidUseOfNullException;
4371 03 Jul 08 nicklas 46 import net.sf.basedb.core.data.FeatureData;
1575 07 Nov 05 nicklas 47 import net.sf.basedb.core.data.RawData;
1575 07 Nov 05 nicklas 48 import net.sf.basedb.core.data.ReporterData;
1576 08 Nov 05 nicklas 49 import net.sf.basedb.core.query.Selects;
1577 08 Nov 05 nicklas 50 import net.sf.basedb.core.query.Select;
1577 08 Nov 05 nicklas 51 import net.sf.basedb.core.query.Expression;
1576 08 Nov 05 nicklas 52 import net.sf.basedb.core.query.Dynamic;
1576 08 Nov 05 nicklas 53 import net.sf.basedb.core.query.Aggregations;
1577 08 Nov 05 nicklas 54 import net.sf.basedb.core.query.Expressions;
4118 01 Feb 08 nicklas 55 import net.sf.basedb.core.signal.ThreadSignalHandler;
1581 09 Nov 05 nicklas 56 import net.sf.basedb.util.jep.Jep;
1581 09 Nov 05 nicklas 57 import net.sf.basedb.util.jep.RawFunction;
1915 06 Feb 06 nicklas 58 import net.sf.basedb.util.jep.MeanFunction;
1575 07 Nov 05 nicklas 59
1575 07 Nov 05 nicklas 60 import java.util.Collection;
1575 07 Nov 05 nicklas 61 import java.util.Collections;
1575 07 Nov 05 nicklas 62 import java.util.Map;
1575 07 Nov 05 nicklas 63 import java.util.HashMap;
1576 08 Nov 05 nicklas 64 import java.sql.SQLException;
1575 07 Nov 05 nicklas 65
1575 07 Nov 05 nicklas 66 import org.nfunk.jep.JEP;
1575 07 Nov 05 nicklas 67
1575 07 Nov 05 nicklas 68 /**
1575 07 Nov 05 nicklas 69   This is a utility class for plugins that calculate intensities
1575 07 Nov 05 nicklas 70   from raw spot data.
1575 07 Nov 05 nicklas 71   
1575 07 Nov 05 nicklas 72   @base.modified $Date$
1575 07 Nov 05 nicklas 73   @author Nicklas
1575 07 Nov 05 nicklas 74   @version 2.0
1575 07 Nov 05 nicklas 75 */
1575 07 Nov 05 nicklas 76 public class IntensityCalculatorUtil
1575 07 Nov 05 nicklas 77 {
1575 07 Nov 05 nicklas 78
2726 12 Oct 06 nicklas 79   public static boolean checkSameArrayDesign(Collection<RawBioAssay> rawBioAssays, boolean throwException)
2726 12 Oct 06 nicklas 80   {
2726 12 Oct 06 nicklas 81     ArrayDesign ad = null;
2726 12 Oct 06 nicklas 82     for (RawBioAssay rba : rawBioAssays)
2726 12 Oct 06 nicklas 83     {
2726 12 Oct 06 nicklas 84       if (ad == null)
2726 12 Oct 06 nicklas 85       {
2726 12 Oct 06 nicklas 86         ad = rba.getArrayDesign();
2726 12 Oct 06 nicklas 87         if (ad == null)
2726 12 Oct 06 nicklas 88         {
2726 12 Oct 06 nicklas 89           if (throwException)
2726 12 Oct 06 nicklas 90           {
2726 12 Oct 06 nicklas 91             throw new InvalidDataException("Raw bioassay doesn't have an array design: "+rba);
2726 12 Oct 06 nicklas 92           }
2726 12 Oct 06 nicklas 93           return false;
2726 12 Oct 06 nicklas 94         }
2726 12 Oct 06 nicklas 95         else if (!ad.hasFeatures())
2726 12 Oct 06 nicklas 96         {
2726 12 Oct 06 nicklas 97           if (throwException)
2726 12 Oct 06 nicklas 98           {
2726 12 Oct 06 nicklas 99             throw new InvalidDataException("Array design doesn't have any features: "+ad);
2726 12 Oct 06 nicklas 100           }
2726 12 Oct 06 nicklas 101           return false;
2726 12 Oct 06 nicklas 102         }
2726 12 Oct 06 nicklas 103       }
2726 12 Oct 06 nicklas 104       else
2726 12 Oct 06 nicklas 105       {
2726 12 Oct 06 nicklas 106         if (!ad.equals(rba.getArrayDesign()))
2726 12 Oct 06 nicklas 107         {
2726 12 Oct 06 nicklas 108           if (throwException)
2726 12 Oct 06 nicklas 109           {
2726 12 Oct 06 nicklas 110             throw new InvalidDataException("Raw bioassay uses a different array design: " +rba);
2726 12 Oct 06 nicklas 111           }
2726 12 Oct 06 nicklas 112           return false;
2726 12 Oct 06 nicklas 113         }
2726 12 Oct 06 nicklas 114       }
2726 12 Oct 06 nicklas 115     }  
2726 12 Oct 06 nicklas 116     return true;
2726 12 Oct 06 nicklas 117   }
2726 12 Oct 06 nicklas 118   
1575 07 Nov 05 nicklas 119   /**
1575 07 Nov 05 nicklas 120     Create a root bioassayset for a given experiment and set of raw bioassays.
1584 10 Nov 05 nicklas 121     The raw bioassays must be part of the experiment. The intensities for
1584 10 Nov 05 nicklas 122     the bioassays in the new bioassayset are calculated by adding Query
1584 10 Nov 05 nicklas 123     API {@link Expression} objects to a {@link DynamicRawDataQuery}.
1584 10 Nov 05 nicklas 124     This method is very fast since most of the work is done on the
1584 10 Nov 05 nicklas 125     database only.
1584 10 Nov 05 nicklas 126     <p>
1584 10 Nov 05 nicklas 127     This method can only be used if all raw bioassays are connected
1584 10 Nov 05 nicklas 128     to the same {@link ArrayDesign} (which must have features), since it
1584 10 Nov 05 nicklas 129     doesn't supports remapping of the position numbers.
1575 07 Nov 05 nicklas 130     
1575 07 Nov 05 nicklas 131     @param dc The <code>DbControl</code> object to use for database acecss
1575 07 Nov 05 nicklas 132     @param experiment The experiment where the new bioassayset should be placed
5228 03 Feb 10 nicklas 133     @param rawBioAssays The raw bioassays to create bioassays for, or null
5228 03 Feb 10 nicklas 134       to use all raw bioassays
1584 10 Nov 05 nicklas 135     @param job The job that the transformation should be linked to, or
1575 07 Nov 05 nicklas 136       null if this method isn't called by a job
1584 10 Nov 05 nicklas 137     @param intensities An array containing the expressions for 
1584 10 Nov 05 nicklas 138       calculating the intensities, one for each channel
4020 29 Nov 07 martin 139      @param progress Optional progress reporter if the caller is interested in the progress 
4020 29 Nov 07 martin 140      @return A {@link net.sf.basedb.core.BioAssaySet} object. 
1584 10 Nov 05 nicklas 141     @throws InvalidDataException If any of the required parameters is null
1584 10 Nov 05 nicklas 142       or if the raw bioassays arn't using the same array design or
1584 10 Nov 05 nicklas 143       if the number of intensity expressions doesn't match the number of channels
4020 29 Nov 07 martin 144      @throws BaseException If anything else goes wrong in the process. 
1575 07 Nov 05 nicklas 145   */
1575 07 Nov 05 nicklas 146   public static BioAssaySet createRootBioAssaySet(DbControl dc, Experiment experiment, 
1873 31 Jan 06 nicklas 147     Collection<RawBioAssay> rawBioAssays, Job job, Expression[] intensities, ProgressReporter progress)
1577 08 Nov 05 nicklas 148     throws InvalidDataException, BaseException
1577 08 Nov 05 nicklas 149   {
1577 08 Nov 05 nicklas 150
1577 08 Nov 05 nicklas 151     if (experiment == null) throw new InvalidUseOfNullException("experiment");
1577 08 Nov 05 nicklas 152     if (intensities == null) throw new InvalidUseOfNullException("intensities");
1577 08 Nov 05 nicklas 153   
1585 10 Nov 05 nicklas 154     if (intensities.length != experiment.getRawDataType().getChannels())
1585 10 Nov 05 nicklas 155     {
2790 24 Oct 06 nicklas 156       throw new InvalidDataException("Wrong number of intensity expressions: "+intensities.length+
1585 10 Nov 05 nicklas 157         "; expected "+experiment.getRawDataType().getChannels());
1585 10 Nov 05 nicklas 158     }
1585 10 Nov 05 nicklas 159     
2726 12 Oct 06 nicklas 160     checkSameArrayDesign(rawBioAssays, true);
4118 01 Feb 08 nicklas 161     ThreadSignalHandler.checkInterrupted();
1584 10 Nov 05 nicklas 162     
5228 03 Feb 10 nicklas 163     // Load raw bioassays
5228 03 Feb 10 nicklas 164     if (rawBioAssays == null || rawBioAssays.size() == 0)
5228 03 Feb 10 nicklas 165     {
5228 03 Feb 10 nicklas 166       rawBioAssays = experiment.getRawBioAssays().list(dc);
5228 03 Feb 10 nicklas 167     }
5228 03 Feb 10 nicklas 168
1577 08 Nov 05 nicklas 169     // Create Transformation
1577 08 Nov 05 nicklas 170     Transformation t = experiment.newTransformation(job, rawBioAssays);
2492 08 Aug 06 nicklas 171     dc.saveItem(t);
1577 08 Nov 05 nicklas 172
1577 08 Nov 05 nicklas 173     // Create the root BioAssaySet in a new data cube and layer
2492 08 Aug 06 nicklas 174     BioAssaySet root = t.newProduct("new", "new", false);
1577 08 Nov 05 nicklas 175     root.setName("New root bioassayset");
1577 08 Nov 05 nicklas 176     dc.saveItem(root);
1577 08 Nov 05 nicklas 177
1577 08 Nov 05 nicklas 178     // Create the batchers we need
1577 08 Nov 05 nicklas 179     PositionBatcher posBatcher = root.getPositionBatcher();
1577 08 Nov 05 nicklas 180     SpotBatcher spotBatcher = root.getSpotBatcher();
1577 08 Nov 05 nicklas 181     MappingBatcher rawBatcher = root.getMappingBatcher();
1577 08 Nov 05 nicklas 182
1873 31 Jan 06 nicklas 183     int totalRawBioAssays = rawBioAssays.size();
1873 31 Jan 06 nicklas 184     int rawBioAssaysDone = 0;
2726 12 Oct 06 nicklas 185     // Spot + raw mapping for each bioassay + reporter mapping
2790 24 Oct 06 nicklas 186     float stepSize = 90f / (float)(totalRawBioAssays * 2 + 1);
2726 12 Oct 06 nicklas 187     float percent = 5;
1791 18 Jan 06 gregory 188     
1577 08 Nov 05 nicklas 189     boolean first = true;
1577 08 Nov 05 nicklas 190     for (RawBioAssay rba : rawBioAssays)
1577 08 Nov 05 nicklas 191     {
1577 08 Nov 05 nicklas 192       // Create a new BioAssay for each raw bioassay
1577 08 Nov 05 nicklas 193       BioAssay ba = root.newRootBioAssay(Collections.singleton(rba));
1952 09 Feb 06 gregory 194       ba.setName(rba.getName());
1577 08 Nov 05 nicklas 195       dc.saveItem(ba);
2726 12 Oct 06 nicklas 196
2726 12 Oct 06 nicklas 197       // The query for this raw bioassay
2726 12 Oct 06 nicklas 198       DynamicRawDataQuery query = rba.getDynamicRawData();
1577 08 Nov 05 nicklas 199       Select column = Selects.expression(Expressions.integer(ba.getDataCubeColumnNo()), VirtualColumn.COLUMN.getName());
1577 08 Nov 05 nicklas 200       Select position = Dynamic.selectRawData("position");
2726 12 Oct 06 nicklas 201
2726 12 Oct 06 nicklas 202       // Map reporters to positions
2726 12 Oct 06 nicklas 203       if (first)
2726 12 Oct 06 nicklas 204       {
2726 12 Oct 06 nicklas 205         if (progress != null)
2726 12 Oct 06 nicklas 206         {
3820 12 Oct 07 nicklas 207           progress.display((int)percent, "Mapping reporters for " + rba.getNumDbSpots() + " spots");
2726 12 Oct 06 nicklas 208           percent += stepSize;
2726 12 Oct 06 nicklas 209         }
2726 12 Oct 06 nicklas 210         query.select(position);
2726 12 Oct 06 nicklas 211         query.select(Dynamic.selectRawData("reporter"));
2726 12 Oct 06 nicklas 212         posBatcher.insert(query);
2726 12 Oct 06 nicklas 213         query.reset();
2726 12 Oct 06 nicklas 214         first = false;
4118 01 Feb 08 nicklas 215         ThreadSignalHandler.checkInterrupted();
2726 12 Oct 06 nicklas 216       }
2726 12 Oct 06 nicklas 217
2726 12 Oct 06 nicklas 218       // Calculate spot intensities
2726 12 Oct 06 nicklas 219       if (progress != null)
2726 12 Oct 06 nicklas 220       {
2726 12 Oct 06 nicklas 221         progress.display((int)percent, "Calculating spot intensities; "+rawBioAssaysDone+
2726 12 Oct 06 nicklas 222           " of " + totalRawBioAssays + " raw bioassay(s) done...");
2726 12 Oct 06 nicklas 223         percent += stepSize;
2726 12 Oct 06 nicklas 224       }
1577 08 Nov 05 nicklas 225       query.select(column);
1577 08 Nov 05 nicklas 226       query.select(position);
1577 08 Nov 05 nicklas 227       for (int i = 0; i < intensities.length; ++i)
1577 08 Nov 05 nicklas 228       {
4915 30 Apr 09 nicklas 229         query.select(Selects.expression(intensities[i], "ch" + i));
1577 08 Nov 05 nicklas 230       }
1577 08 Nov 05 nicklas 231       spotBatcher.insert(query);
1577 08 Nov 05 nicklas 232       query.reset();
4118 01 Feb 08 nicklas 233       ThreadSignalHandler.checkInterrupted();
1577 08 Nov 05 nicklas 234       
2726 12 Oct 06 nicklas 235       // Map spot data to raw data
2726 12 Oct 06 nicklas 236       if (progress != null)
2726 12 Oct 06 nicklas 237       {
2726 12 Oct 06 nicklas 238         progress.display((int)percent, "Mapping raw data; "+rawBioAssaysDone+" of " + 
2726 12 Oct 06 nicklas 239           totalRawBioAssays + " raw bioassay(s) done...");
2726 12 Oct 06 nicklas 240         percent += stepSize;
2726 12 Oct 06 nicklas 241       }
1577 08 Nov 05 nicklas 242       query.select(column);
1577 08 Nov 05 nicklas 243       query.select(position);
1577 08 Nov 05 nicklas 244       query.select(Dynamic.selectRawData("id"));
1577 08 Nov 05 nicklas 245       rawBatcher.insert(query);
1577 08 Nov 05 nicklas 246       query.reset();
2726 12 Oct 06 nicklas 247       rawBioAssaysDone++;
4118 01 Feb 08 nicklas 248       ThreadSignalHandler.checkInterrupted();
1577 08 Nov 05 nicklas 249     }
1577 08 Nov 05 nicklas 250     
1577 08 Nov 05 nicklas 251     // Close batchers and clean up
1577 08 Nov 05 nicklas 252     posBatcher.close();
1577 08 Nov 05 nicklas 253     spotBatcher.close();
1577 08 Nov 05 nicklas 254     rawBatcher.close();
1873 31 Jan 06 nicklas 255
1577 08 Nov 05 nicklas 256     // Return  
1577 08 Nov 05 nicklas 257     return root;
1577 08 Nov 05 nicklas 258   }
1577 08 Nov 05 nicklas 259
1577 08 Nov 05 nicklas 260   /**
1577 08 Nov 05 nicklas 261     Create a root bioassayset for a given experiment and set of raw bioassays.
1584 10 Nov 05 nicklas 262     The raw bioassays must be part of the experiment. The intensities for the
1584 10 Nov 05 nicklas 263     bioassays in the new bioassayset are calculated by a {@link IntensityCalculator}
1584 10 Nov 05 nicklas 264     object. This method is not as fast as the other <code>createRootBioAssay</code>
1584 10 Nov 05 nicklas 265     method since all raw data must be loaded into memory and the inserted into
1584 10 Nov 05 nicklas 266     the database one row at a time.
1584 10 Nov 05 nicklas 267     <p>
1584 10 Nov 05 nicklas 268     On the other hand, this method supports remapping of the position numbers
1584 10 Nov 05 nicklas 269     so it is not necessary for all raw bioassays to use the same array design.
1584 10 Nov 05 nicklas 270
1577 08 Nov 05 nicklas 271     @param dc The <code>DbControl</code> object to use for database acecss
1577 08 Nov 05 nicklas 272     @param experiment The experiment where the new bioassayset should be placed
5228 03 Feb 10 nicklas 273     @param rawBioAssays The raw bioassays to create bioassays for, or null
5228 03 Feb 10 nicklas 274       to use all raw bioassays
1584 10 Nov 05 nicklas 275     @param job The job that the transformation should be linked to, or
1577 08 Nov 05 nicklas 276       null if this method isn't called by a job
1584 10 Nov 05 nicklas 277     @param iCalc An object implementing the <code>IntensityCalculator</code>
1577 08 Nov 05 nicklas 278       interface
1873 31 Jan 06 nicklas 279     @param progress A progress reporter object, or null if the caller isn't
1873 31 Jan 06 nicklas 280       interested in the progress
4020 29 Nov 07 martin 281      @return A {@link net.sf.basedb.core.BioAssaySet} object 
4020 29 Nov 07 martin 282      @throws InvalidDataException If any of the required parameters is null. 
4020 29 Nov 07 martin 283      @throws BaseException If there is another error.
1577 08 Nov 05 nicklas 284   */
1577 08 Nov 05 nicklas 285   public static BioAssaySet createRootBioAssaySet(DbControl dc, Experiment experiment, 
1873 31 Jan 06 nicklas 286     Collection<RawBioAssay> rawBioAssays, Job job, IntensityCalculator iCalc, ProgressReporter progress)
1575 07 Nov 05 nicklas 287     throws InvalidDataException, BaseException
1575 07 Nov 05 nicklas 288   {
1575 07 Nov 05 nicklas 289
1575 07 Nov 05 nicklas 290     if (experiment == null) throw new InvalidUseOfNullException("experiment");
1575 07 Nov 05 nicklas 291     if (iCalc == null) throw new InvalidUseOfNullException("iCalc");
1575 07 Nov 05 nicklas 292   
1575 07 Nov 05 nicklas 293     // Create Transformation
1575 07 Nov 05 nicklas 294     Transformation t = experiment.newTransformation(job, rawBioAssays);
2492 08 Aug 06 nicklas 295     dc.saveItem(t);
1873 31 Jan 06 nicklas 296     
1873 31 Jan 06 nicklas 297     if (progress != null) progress.display(5, "Initialising...");
1575 07 Nov 05 nicklas 298
5228 03 Feb 10 nicklas 299     // Load raw bioassays
5228 03 Feb 10 nicklas 300     if (rawBioAssays == null || rawBioAssays.size() == 0)
5228 03 Feb 10 nicklas 301     {
5228 03 Feb 10 nicklas 302       rawBioAssays = experiment.getRawBioAssays().list(dc);
5228 03 Feb 10 nicklas 303     }
5228 03 Feb 10 nicklas 304
1575 07 Nov 05 nicklas 305     // Create the root BioAssaySet in a new data cube and layer
2492 08 Aug 06 nicklas 306     BioAssaySet root = t.newProduct("new", "new", false);
1575 07 Nov 05 nicklas 307     root.setName("New root bioassayset");
1575 07 Nov 05 nicklas 308     dc.saveItem(root);
1575 07 Nov 05 nicklas 309
1575 07 Nov 05 nicklas 310     // Create the batchers we need
1575 07 Nov 05 nicklas 311     PositionBatcher posBatcher = root.getPositionBatcher();
1575 07 Nov 05 nicklas 312     SpotBatcher spotBatcher = root.getSpotBatcher();
1575 07 Nov 05 nicklas 313     MappingBatcher rawBatcher = root.getMappingBatcher();
1575 07 Nov 05 nicklas 314
2492 08 Aug 06 nicklas 315     // A cache holding already seen position --> reporter mappings
4371 03 Jul 08 nicklas 316     Map<Integer, Integer> positionCache = new HashMap<Integer, Integer>();
4371 03 Jul 08 nicklas 317     // A cache holding already seen feature->position mappings
4371 03 Jul 08 nicklas 318     Map<Integer, Integer> featureCache = new HashMap<Integer, Integer>();
1575 07 Nov 05 nicklas 319     
1576 08 Nov 05 nicklas 320     // The max position number in the raw bioassays -- only calculated if needed (see below)
1576 08 Nov 05 nicklas 321     int maxPosition = 0;
1576 08 Nov 05 nicklas 322     
4095 21 Jan 08 enell 323     long spotsDone = 0; // Number of calculated spots
1873 31 Jan 06 nicklas 324     int totalSpots = 0; // Total number of spots
1873 31 Jan 06 nicklas 325     int interval = 0; // Interval between reporter to the progress reporter
1873 31 Jan 06 nicklas 326     if (progress != null) 
1873 31 Jan 06 nicklas 327     {
1873 31 Jan 06 nicklas 328       for (RawBioAssay rba : rawBioAssays)
1873 31 Jan 06 nicklas 329       {
3820 12 Oct 07 nicklas 330         totalSpots += rba.getNumDbSpots();
1873 31 Jan 06 nicklas 331       }
2726 12 Oct 06 nicklas 332       interval = totalSpots / 50;
1873 31 Jan 06 nicklas 333       if (interval < 10) interval = 10;
2726 12 Oct 06 nicklas 334       progress.display(10, "Calculating spot intensities (0 done)...");
1873 31 Jan 06 nicklas 335     }
2497 08 Aug 06 nicklas 336     
1575 07 Nov 05 nicklas 337     for (RawBioAssay rba : rawBioAssays)
1575 07 Nov 05 nicklas 338     {
1575 07 Nov 05 nicklas 339       // Create a new BioAssay for each raw bioassay
1575 07 Nov 05 nicklas 340       BioAssay ba = root.newRootBioAssay(Collections.singleton(rba));
1952 09 Feb 06 gregory 341       ba.setName(rba.getName());
1575 07 Nov 05 nicklas 342       dc.saveItem(ba);
1575 07 Nov 05 nicklas 343     
1575 07 Nov 05 nicklas 344       // This is the column coordinate in the datacube
1575 07 Nov 05 nicklas 345       short column = ba.getDataCubeColumnNo();
1575 07 Nov 05 nicklas 346     
1575 07 Nov 05 nicklas 347       // Load raw data for current raw bioassay
2497 08 Aug 06 nicklas 348       DataQuery<RawData> rawQuery = rba.getRawData();
2497 08 Aug 06 nicklas 349       DataResultIterator<RawData> rawData = rawQuery.iterate(dc);
4118 01 Feb 08 nicklas 350       ThreadSignalHandler.checkInterrupted();
1575 07 Nov 05 nicklas 351       while (rawData.hasNext())
1575 07 Nov 05 nicklas 352       {
4118 01 Feb 08 nicklas 353         ThreadSignalHandler.checkInterrupted();
1575 07 Nov 05 nicklas 354         RawData rd = rawData.next();
1575 07 Nov 05 nicklas 355         int position = rd.getPosition();
1575 07 Nov 05 nicklas 356         ReporterData reporter = rd.getReporter();
4371 03 Jul 08 nicklas 357         FeatureData feature = rd.getFeature();
1575 07 Nov 05 nicklas 358         
1575 07 Nov 05 nicklas 359         // Calculate intensities
1575 07 Nov 05 nicklas 360         float[] intensities = iCalc.calculateIntensities(rd);
2497 08 Aug 06 nicklas 361         
1575 07 Nov 05 nicklas 362         if (intensities != null)
1575 07 Nov 05 nicklas 363         {
4371 03 Jul 08 nicklas 364           // If we have a feature check the feature cache if a position has
4371 03 Jul 08 nicklas 365           // already been assigned to it
4371 03 Jul 08 nicklas 366           Integer cachedFeaturePosition = feature == null ?
4371 03 Jul 08 nicklas 367             null : featureCache.get(feature.getId());
4371 03 Jul 08 nicklas 368           if (cachedFeaturePosition != null)
1575 07 Nov 05 nicklas 369           {
4371 03 Jul 08 nicklas 370             position = cachedFeaturePosition;
4371 03 Jul 08 nicklas 371           }
4371 03 Jul 08 nicklas 372           else
4371 03 Jul 08 nicklas 373           {
5689 11 Aug 11 nicklas 374             // We have no feature or the feature has no position. 
5689 11 Aug 11 nicklas 375             // Check if the position has same reporter
4371 03 Jul 08 nicklas 376             if (positionCache.containsKey(position))
1575 07 Nov 05 nicklas 377             {
4371 03 Jul 08 nicklas 378               Integer cachedReporterId = positionCache.get(position);
4371 03 Jul 08 nicklas 379               Integer reporterId = reporter == null ? null : reporter.getId();
4371 03 Jul 08 nicklas 380               if (reporterId == null && cachedReporterId != null || 
4371 03 Jul 08 nicklas 381                 reporterId != null && !reporterId.equals(cachedReporterId))
1576 08 Nov 05 nicklas 382               {
4371 03 Jul 08 nicklas 383                 // Reporters don't match, we need a new position
4371 03 Jul 08 nicklas 384                 if (maxPosition == 0)
1576 08 Nov 05 nicklas 385                 {
4371 03 Jul 08 nicklas 386                   // We must find the maximum position among all raw bioassays unless we havn't done that already
4371 03 Jul 08 nicklas 387                   try
1576 08 Nov 05 nicklas 388                   {
4371 03 Jul 08 nicklas 389                     for (RawBioAssay rbaMax : rawBioAssays)
4371 03 Jul 08 nicklas 390                     {
4371 03 Jul 08 nicklas 391                       DynamicRawDataQuery query = rbaMax.getDynamicRawData();
4371 03 Jul 08 nicklas 392                       query.select(Selects.expression(Aggregations.max(Dynamic.rawData("position")), "maxPosition"));
4371 03 Jul 08 nicklas 393                       DynamicResultIterator result = query.iterate(dc);
4371 03 Jul 08 nicklas 394                       int thisMaxPosition = result.next().getInt(1);
4371 03 Jul 08 nicklas 395                       if (thisMaxPosition > maxPosition) maxPosition = thisMaxPosition;
4371 03 Jul 08 nicklas 396                       result.close();
4371 03 Jul 08 nicklas 397                     }
1576 08 Nov 05 nicklas 398                   }
4371 03 Jul 08 nicklas 399                   catch (SQLException ex)
4371 03 Jul 08 nicklas 400                   {
4371 03 Jul 08 nicklas 401                     throw new BaseException(ex);
4371 03 Jul 08 nicklas 402                   }
1576 08 Nov 05 nicklas 403                 }
4371 03 Jul 08 nicklas 404                 maxPosition++;
4371 03 Jul 08 nicklas 405                 position = maxPosition;
4371 03 Jul 08 nicklas 406                 posBatcher.insert(position, reporter);
4371 03 Jul 08 nicklas 407                 positionCache.put(position, reporterId);
4371 03 Jul 08 nicklas 408                 if (feature != null) featureCache.put(feature.getId(), position);
1576 08 Nov 05 nicklas 409               }
4371 03 Jul 08 nicklas 410             }
4371 03 Jul 08 nicklas 411             else
4371 03 Jul 08 nicklas 412             {
4371 03 Jul 08 nicklas 413               // This position is not used, map it to the current reporter
1576 08 Nov 05 nicklas 414               posBatcher.insert(position, reporter);
4371 03 Jul 08 nicklas 415               positionCache.put(position, reporter == null ? null : reporter.getId());
4371 03 Jul 08 nicklas 416               if (feature != null) featureCache.put(feature.getId(), position);
1575 07 Nov 05 nicklas 417             }
1575 07 Nov 05 nicklas 418           }
4371 03 Jul 08 nicklas 419           
1575 07 Nov 05 nicklas 420           // Insert mapping to raw data spot
1575 07 Nov 05 nicklas 421           rawBatcher.insert(column, position, rd);
1575 07 Nov 05 nicklas 422           // Insert intensities
1575 07 Nov 05 nicklas 423           spotBatcher.insert(column, position, intensities);
1575 07 Nov 05 nicklas 424         }
1873 31 Jan 06 nicklas 425         if (progress != null)
1873 31 Jan 06 nicklas 426         {
1873 31 Jan 06 nicklas 427           spotsDone++;
1873 31 Jan 06 nicklas 428           if (spotsDone % interval == 0)
1873 31 Jan 06 nicklas 429           {
4095 21 Jan 08 enell 430             int percent = 10 + (int)((90L * spotsDone) / totalSpots);
2726 12 Oct 06 nicklas 431             progress.display(percent, "Calculating spot intensities ("+spotsDone+" done)...");
1873 31 Jan 06 nicklas 432           }
1873 31 Jan 06 nicklas 433         }
1575 07 Nov 05 nicklas 434       }
1575 07 Nov 05 nicklas 435     }
2726 12 Oct 06 nicklas 436     if (progress != null) progress.display(100, "Calculating spot intensities ("+totalSpots+" done)...\n");
5689 11 Aug 11 nicklas 437   
1575 07 Nov 05 nicklas 438     // Close batchers and clean up
1575 07 Nov 05 nicklas 439     positionCache.clear();
1575 07 Nov 05 nicklas 440     posBatcher.close();
1575 07 Nov 05 nicklas 441     spotBatcher.close();
1575 07 Nov 05 nicklas 442     rawBatcher.close();
1575 07 Nov 05 nicklas 443     
1575 07 Nov 05 nicklas 444     // Return  
1575 07 Nov 05 nicklas 445     return root;
1575 07 Nov 05 nicklas 446   }
1575 07 Nov 05 nicklas 447   
1575 07 Nov 05 nicklas 448   /**
1584 10 Nov 05 nicklas 449     Create an intensity calculator from mathematical expressions in strings. The expressions
4020 29 Nov 07 martin 450     may contain most regular mathematical operators, including <code>+, -, / , *, ln</code>.
2322 24 May 06 nicklas 451     Use the name of a raw() function to include the value of that property, or the
2322 24 May 06 nicklas 452     mean() function to include the mean of a property (over the whole data set). Examples
4020 29 Nov 07 martin 453     of valid expressions are (assuming that the properties has been defined for the
2322 24 May 06 nicklas 454     raw data type in question; see {@link RawDataType#getProperties()}):
1575 07 Nov 05 nicklas 455     <pre class="code">
2322 24 May 06 nicklas 456 raw('ch1FgMean') - raw('ch1BgMean')
2322 24 May 06 nicklas 457 raw('ch1FgMedian') - mean('ch1BgMedian')
1575 07 Nov 05 nicklas 458 </pre>
1575 07 Nov 05 nicklas 459     
1590 11 Nov 05 nicklas 460     @param dc The DbControl used for database access
4020 29 Nov 07 martin 461     @param rawDataType The raw data type to create the intensity calculator for. Null is not allowed.
2322 24 May 06 nicklas 462     @param jepFormulas The mathematical expressions in a form supported by JEP,
4020 29 Nov 07 martin 463       one expression is required for each channel. Null is not allowed
1575 07 Nov 05 nicklas 464     @return An intensity calculator object
4020 29 Nov 07 martin 465      @throws InvalidDataException If any of the required arguments is null or
4020 29 Nov 07 martin 466        if the formulas length does not match the number of channels. 
4020 29 Nov 07 martin 467      @throws BaseException If there is another error.
1575 07 Nov 05 nicklas 468     @see RawDataType#getIntensityFormulas()
1584 10 Nov 05 nicklas 469     @see net.sf.basedb.core.IntensityFormula
1584 10 Nov 05 nicklas 470     @see net.sf.basedb.util.jep.Jep
1575 07 Nov 05 nicklas 471   */
1590 11 Nov 05 nicklas 472   public static IntensityCalculator createJepIntensityCalculator(DbControl dc, RawDataType rawDataType, String... jepFormulas)
1584 10 Nov 05 nicklas 473     throws InvalidDataException, BaseException
1575 07 Nov 05 nicklas 474   {
1575 07 Nov 05 nicklas 475     if (rawDataType == null) throw new InvalidUseOfNullException("rawDataType");
1575 07 Nov 05 nicklas 476     if (jepFormulas == null || jepFormulas.length == 0) throw new InvalidUseOfNullException("jepFormulas");
1575 07 Nov 05 nicklas 477   
1575 07 Nov 05 nicklas 478     int channels = rawDataType.getChannels();
1575 07 Nov 05 nicklas 479     if (channels != jepFormulas.length) 
1575 07 Nov 05 nicklas 480     {
1575 07 Nov 05 nicklas 481       throw new InvalidDataException("Invalid number of formulas, expected " + channels + "; found "+jepFormulas.length);
1575 07 Nov 05 nicklas 482     }
1575 07 Nov 05 nicklas 483     JEP[] jeps = new JEP[channels];
2993 01 Dec 06 nicklas 484     RawFunction raw = new RawFunction(dc, rawDataType, true);
1915 06 Feb 06 nicklas 485     MeanFunction mean = new MeanFunction(dc);
1575 07 Nov 05 nicklas 486     for (int i = 0; i < channels; ++i)
1575 07 Nov 05 nicklas 487     {
1575 07 Nov 05 nicklas 488       if (jepFormulas[i] == null)
1575 07 Nov 05 nicklas 489       {
1575 07 Nov 05 nicklas 490         throw new InvalidUseOfNullException("jepFormulas["+i+"]");
1575 07 Nov 05 nicklas 491       }
1915 06 Feb 06 nicklas 492       jeps[i] = Jep.newJep(jepFormulas[i], raw, mean);
1575 07 Nov 05 nicklas 493     }
1915 06 Feb 06 nicklas 494     return new JepIntensityCalculatorImpl(jeps, raw, mean);
1575 07 Nov 05 nicklas 495   }
1575 07 Nov 05 nicklas 496
2726 12 Oct 06 nicklas 497   public static Expression[] createJepExpressions(DbControl dc, RawDataType rawDataType, String[] formulas)
2726 12 Oct 06 nicklas 498   {
2726 12 Oct 06 nicklas 499     Expression[] expressions = new Expression[formulas.length];
2993 01 Dec 06 nicklas 500     RawFunction raw = new RawFunction(dc, rawDataType, true);
2726 12 Oct 06 nicklas 501     MeanFunction mean = new MeanFunction(dc);
2726 12 Oct 06 nicklas 502     for (int i = 0; i < formulas.length; ++i)
2726 12 Oct 06 nicklas 503     {
2726 12 Oct 06 nicklas 504       expressions[i] = Jep.formulaToExpression(formulas[i], raw, mean);
2726 12 Oct 06 nicklas 505     }
2726 12 Oct 06 nicklas 506     return expressions;
2726 12 Oct 06 nicklas 507   }
2726 12 Oct 06 nicklas 508   
1584 10 Nov 05 nicklas 509   /**
1584 10 Nov 05 nicklas 510     An implementation of the {@link IntensityCalculator} interface using
1584 10 Nov 05 nicklas 511     JEP expressions.
1584 10 Nov 05 nicklas 512     @see net.sf.basedb.util.jep.Jep
1584 10 Nov 05 nicklas 513   */
1575 07 Nov 05 nicklas 514   private static class JepIntensityCalculatorImpl
1575 07 Nov 05 nicklas 515     implements IntensityCalculator
1575 07 Nov 05 nicklas 516   {
1575 07 Nov 05 nicklas 517     private final JEP[] jeps;
1581 09 Nov 05 nicklas 518     private final RawFunction raw;
1915 06 Feb 06 nicklas 519     private final MeanFunction mean;
1575 07 Nov 05 nicklas 520     
1915 06 Feb 06 nicklas 521     private JepIntensityCalculatorImpl(JEP[] jeps, RawFunction raw, MeanFunction mean)
1575 07 Nov 05 nicklas 522     {
1575 07 Nov 05 nicklas 523       this.jeps = jeps;
1581 09 Nov 05 nicklas 524       this.raw = raw;
1915 06 Feb 06 nicklas 525       this.mean = mean;
1575 07 Nov 05 nicklas 526     }
1575 07 Nov 05 nicklas 527   
1584 10 Nov 05 nicklas 528     /*
1584 10 Nov 05 nicklas 529       From the IntensityCalculator interface
1584 10 Nov 05 nicklas 530       -------------------------------------------
1584 10 Nov 05 nicklas 531     */
6127 14 Sep 12 nicklas 532     @Override
1575 07 Nov 05 nicklas 533     public float[] calculateIntensities(RawData rawData)
1575 07 Nov 05 nicklas 534     {
1575 07 Nov 05 nicklas 535       float[] intensities = new float[jeps.length];
1575 07 Nov 05 nicklas 536       for (int i = 0; i < jeps.length; ++i)
1575 07 Nov 05 nicklas 537       {
1575 07 Nov 05 nicklas 538         JEP jep = jeps[i];
1581 09 Nov 05 nicklas 539         raw.setRawData(rawData);
1915 06 Feb 06 nicklas 540         mean.setRawData(rawData);
1575 07 Nov 05 nicklas 541         intensities[i] = (float)jep.getValue();
1575 07 Nov 05 nicklas 542       }
1575 07 Nov 05 nicklas 543       return intensities;
1575 07 Nov 05 nicklas 544     }
1584 10 Nov 05 nicklas 545     // -------------------------------------------
1575 07 Nov 05 nicklas 546   }
1575 07 Nov 05 nicklas 547 }

Generated on Sat Oct  7 02:39:49 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)