Changeset 1114 for plugins/base2

Timestamp:

Jun 10, 2009, 12:09:36 AM (14 years ago)

Author:

Jari Häkkinen

Message:

Addresses #113. The Pvalue calculation is in place. Need to figure details of std and average calculations.

Location:

plugins/base2/net.sf.basedb.illumina/trunk

Files:

• LICENSE (modified) (1 diff)
• README_PluginDetails (modified) (1 diff)
• src/net/sf/basedb/illumina/plugins/DetectionPValue.java (modified) (5 diffs)
• src/org (added)
• src/org/apache (added)
• src/org/apache/commons (added)
• src/org/apache/commons/math (added)
• src/org/apache/commons/math/ConvergenceException.java (added)
• src/org/apache/commons/math/MathException.java (added)
• src/org/apache/commons/math/MaxIterationsExceededException.java (added)
• src/org/apache/commons/math/special (added)
• src/org/apache/commons/math/special/Erf.java (added)
• src/org/apache/commons/math/special/Gamma.java (added)
• src/org/apache/commons/math/util (added)
• src/org/apache/commons/math/util/ContinuedFraction.java (added)

plugins/base2/net.sf.basedb.illumina/trunk/LICENSE

@@ -1 @@
-Files are copyright by their respective authors. The contributions to
-files where copyright is not explicitly stated can be traced with the
-source code revision system.
+Source code in src/org/apache/commons is unmodified copies of the
+files retrieved from the repository (revision 783095). The files are
+licensed under the Apache license 2.0. The license is available at
+http://www.apache.org/licenses/LICENSE-2.0
+
+
+All other files are copyright by their respective authors. The
+contributions to files where copyright is not explicitly stated can be
+traced with the baseplugins source code revision system.
 
 This file is part of Illumina plug-in package for BASE.

plugins/base2/net.sf.basedb.illumina/trunk/README_PluginDetails

@@ -55 +55 @@
 current bioassay set.
 
+Whole genome beadchips:
+
 For all signals `i` calculate the detection P-value as
-`Pvalue = 1-i/N` where `i` is the rank of the signal relative to
-the negative controls and `N` is the number of negative controls.
+`Pvalue = 1-i/N` where `R` is the rank of the signal i relative to the
+negative controls and `N` is the number of negative controls.
+
+Others (DASL, miRNA, VeraCode DASL, and Focused Arrays):
+
+For all signals `i` calculate the detection P-value as
+`Pvalue = 1/2 - erf( [i-AvgControl]/StdControl/sqrt(2) )` where
+`AvgControl` is the average intensity of the negative controls,
+`StdControl` is the standard deviation of the ... (to be defined and
+implemented accordingly), and `erf` is the error function
+(http://mathworld.wolfram.com/Erf.html). The error function is used
+for arguments witin the range (-4,4). To save CPU cycles, the function
+value for arguments outside this range is set to -1 and 1,
+respectively.
 
 

plugins/base2/net.sf.basedb.illumina/trunk/src/net/sf/basedb/illumina/plugins/DetectionPValue.java

@@ -67 +67 @@
 import net.sf.basedb.util.Values;
 
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.special.Erf;
+
 import java.sql.SQLException;
 import java.util.ArrayList;
@@ -117 +120 @@
     ( "arraytype",
       "Array type",
-      "BeadStudio like detection P-value calculation for array types:\n\n" +
-      "Illumina Whole Genome BeadChips\n" +
-      "Illumina Focused Arrays (not yet implemented)",
+      "Illumina BeadStudio like detection P-value calculation.\n\n" +
+      "Supported array types:\n\n" +
+      " - Whole Genome BeadChips\n" +
+      " - Others (DASL, miRNA, VeraCode DASL, and Focused Arrays)",
       new StringParameterType(255, "Whole Genome BeadChip", true, 1, 0, 0,
         Arrays.asList(
-          new String[] { "Whole Genome BeadChip", "Focused Array (not yet implemented)" }
+          new String[] { "Whole Genome BeadChip",
+                         "Others" }
         )
       )
@@ -352 +357 @@
 
 
-  float pvalue(float[] sortedArray, int startIndex, int stopIndex, float intensity)
+  float pvalueOthers(float intensityA, float intensityS, float intensity)
+    throws MathException
+  {
+    double arg=(intensity-intensityA)/intensityS/1.414213562;
+    // No need to waste iterations in Erf for large or small arguments.
+    if (arg>4)
+      return 0;
+    else if (arg<-4)
+      return 1;
+    return (float)(0.5 - Erf.erf(arg)/2);
+  }
+
+
+  float pvalueWG(float[] sortedArray, int startIndex, int stopIndex, float intensity)
   {
     int i=stopIndex;
@@ -461 +479 @@
             (Float)job.getValue(outlierParameter.getName()));
 
+          // Calculate average control intensity and std only if
+          // needed, i.e., if not calculating detection P-value for
+          // whole genome beadchips.
+          float cIntAvg=0;
+          float cIntStd=0;
+          if (arrayType.equals("Others")) { // must match arrayTypeParameter
+            float x=0;
+            float xx=0;
+            for (int i=controlStartStop[0]; i<=controlStartStop[1]; ++i)
+            {
+              x+=sortedControlIntensity[i];
+              xx+=sortedControlIntensity[i]*sortedControlIntensity[i];
+            }
+            int n=controlStartStop[1]-controlStartStop[0]+1;
+            cIntAvg=x/n;
+            cIntStd=(float)Math.sqrt((xx-x*x/n)/n);
+          }
+
           checkInterrupted();
 
@@ -472 +508 @@
           {
             SqlResult result = dri.next();
-            batcher.insert(result.getShort(columnIndex), result.getInt(positionIndex),
-              pvalue(sortedControlIntensity, controlStartStop[0],
-              controlStartStop[1], result.getFloat(ch1Index)));
+            float pval=( arrayType.equals("Whole Genome BeadChip") ?
+                         pvalueWG(sortedControlIntensity, controlStartStop[0],
+                                  controlStartStop[1], result.getFloat(ch1Index)) :
+                         pvalueOthers(cIntAvg, cIntStd, result.getFloat(ch1Index)) );
+            batcher.insert(result.getShort(columnIndex), result.getInt(positionIndex), pval);
           }
           dri.close();